Matlab中的傅里叶变换和LTI滤波器以及频率响应

Question

我是 Matlab 的 LTI 信号处理新手，想知道是否有人可以帮助我确定一些基本的东西。我花了好几个小时研究和获取背景信息，但仍然无法找到解决这些问题的明确途径。到目前为止，从头开始，我已经生成了一个所需的信号，并设法使用 fft 函数来生成信号的 DFT：

function x = fourier_rikki(A,t,O)

Fs = 1000;
t = 0:(1/Fs):1;

A = [0.5,0,0.5];
N = (length(A) - 1)/2;
x = zeros(size(t));
f1 = 85;
O1 = 2*pi*f1;

for k = 1:length(A)
x1 = x + A(k)*exp(1i*O1*t*(k-N-1));
end

f2 = 150;
O2 = 2*pi*f2;

for k = 1:length(A);
x2 = x + A(k)*exp(1i*O2*t*(k-N-1));
end

f3 = 330;
O3 = 2*pi*f3;

for k = 1:length(A);
x3 = x + A(k)*exp(1i*O3*t*(k-N-1));
end

signal = x1 + x2 + x3;

figure(1);
subplot(3,1,1);
plot(t, signal);
title('Signal x(t) in the Time Domain');
xlabel('Time (Seconds)');
ylabel('x(t)');

X = fft(signal);  %DFT of the signal

subplot(3,1,2);
plot(t, X);
title('Power Spectrum of Discrete Fourier Transform of x(t)');
xlabel('Time (Seconds)');
ylabel('Power');

f = linspace(0, 1000, length(X)); %?

subplot(3,1,3);
plot(f, abs(X));  %Only want the positive values
title('Spectral Frequency');
xlabel('Frequency (Hz)'); ylabel('Power');

end

在这个阶段，我假设这是正确的：

“使用 1000Hz 的采样频率生成频率为 85,150,330Hz 的信号 - 绘制 1 秒的信号及其离散傅立叶变换。”

下一步是“查找使用傅里叶变换滤除较高和较低频率的 LTI 系统的频率响应”。我一直在尝试创建一个可以做到这一点的 LTI 系统！我必须留下 150Hz 信号，我猜我在 FFT 上执行滤波，可能使用 conv。

我的课程不是编程课程 - 我们没有评估我们的编程技能，而且我的 Matlab 经验极少 - 基本上我们只能靠自己的设备来挣扎，因此我们将不胜感激任何帮助！我正在筛选大量不同的示例并使用“帮助”等搜索 Matlab 函数，但是由于每个示例都不同，并且没有分解所使用的变量，因此解释了为什么选择某些参数/值等。它只是添加混乱。

在我看过的许多（许多）其他人中： http://www.ee.columbia.edu/~ronw/adst-spring2010/lectures/matlab/lecture1.html http://gribblelab.org/scicomp/09_Signals_and_sampling.html 尤其是第 10.4 节。 以及 Matlab Geeks 示例和 Mathworks Matlab 函数说明。 我想可能发生的最糟糕的情况是没有人回答，我继续燃烧我的眼球，直到我设法想出一些东西:) 提前谢谢。

我发现这个带通滤波器代码是一个 Mathworks 示例，这正是需要应用于我的 fft 信号的代码，但我不了解衰减值 Ast 或纹波量 Ap。

n = 0:159;
x = cos(pi/8*n)+cos(pi/2*n)+sin(3*pi/4*n);

d = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2',1/4,3/8,5/8,6/8,60,1,60);
Hd = design(d,'equiripple');

y = filter(Hd,x);
freq = 0:(2*pi)/length(x):pi;
xdft = fft(x);
ydft = fft(y);

plot(freq,abs(xdft(1:length(x)/2+1)));
hold on;
plot(freq,abs(ydft(1:length(x)/2+1)),'r','linewidth',2);
legend('Original Signal','Bandpass Signal');

Answer 1

这里有一些你可以用作参考的东西。我想我明白了你想要做什么的要点。如果您有任何问题，请告诉我。

clear all
close all

Fs = 1000;
t = 0:(1/Fs):1;
N = length(t);

% 85, 150, and 330 Hz converted to radian frequency
w1 = 2*pi*85;
w2 = 2*pi*150;
w3 = 2*pi*330;

% amplitudes
a1 = 1;
a2 = 1.5;
a3 = .75;

% construct time-domain signals
x1 = a1*cos(w1*t);
x2 = a2*cos(w2*t);
x3 = a3*cos(w3*t);

% superposition of 85, 150, and 330 Hz component signals
x = x1 + x2 + x3; 

figure
plot(t(1:100), x(1:100));
title('unfiltered time-domain signal, amplitude vs. time');
ylabel('amplitude');
xlabel('time (seconds)');

% compute discrete Fourier transform of time-domain signal
X = fft(x); 
Xmag = 20*log10(abs(X)); % magnitude spectrum
Xphase = 180*unwrap(angle(X))./pi; % phase spectrum (degrees)
w = 2*pi*(0:N-1)./N; % normalized radian frequency
f = w./(2*pi)*Fs; % radian frequency to Hz
k = 1:N; % bin indices 

% plot magnitude spectrum
figure
plot(f, Xmag)
title('frequency-domain signal, magnitude vs. frequency');
xlabel('frequency (Hz)');
ylabel('magnitude (dB)');

% frequency vector of the filter. attenuates undesired frequency components
% and keeps desired components. 
H = 1e-3*ones(1, length(k));
H(97:223) = 1;
H((end-223):(end-97)) = 1;

% plot magnitude spectrum of signal and filter
figure
plot(k, Xmag)
hold on
plot(k, 20*log10(H), 'r')
title('frequency-domain signal (blue) and filter (red), magnitude vs. bin index');
xlabel('bin index');
ylabel('magnitude (dB)');

% filtering in frequency domain is just multiplication
Y = X.*H;

% plot magnitude spectrum of filtered signal
figure
plot(f, 20*log10(abs(Y)))
title('filtered frequency-domain signal, magnitude vs. frequency');
xlabel('frequency (Hz)');
ylabel('magnitude (dB)');

% use inverse discrete Fourier transform to obtain the filtered time-domain
% signal. This signal is complex due to imperfect symmetry in the
% frequency-domain, however the imaginary components are nearly zero.
y = ifft(Y);

% plot overlay of filtered signal and desired signal
figure
plot(t(1:100), x(1:100), 'r')
hold on
plot(t(1:100), x2(1:100), 'linewidth', 2)
plot(t(1:100), real(y(1:100)), 'g')
title('input signal (red), desired signal (blue), signal extracted via filtering (green)');
ylabel('amplitude');
xlabel('time (seconds)');

这是最终结果...