matlab - 查找带通滤波器的频率响应

Question

我需要执行此MATLAB任务，到目前为止，我已经生成了一个混合信号，该信号由三个频率值的余弦函数组成:86Hz，159Hz和392Hz。我使用了1Khz的采样频率，并生成了一秒钟的数据。然后，我创建了此代码的离散傅立叶变换，该变换显示了(我认为FT的目的是)复合信号中最常见的频率。下面显示了此代码和生成的图:

fs = 1000;
t = (0 : 1/fs : 1)';
f = [ 86, 159, 392 ];
x = sum( cos(2*pi * t * f), 2 );
y = fft(x);
subplot(2,1,1), plot(t,x);
subplot(2,1,2), stem(y,t);

(为了方便起见，我将省略轴和图标题的代码)。这是我到目前为止为报告生成的两个数字。

报告的下一部分向我提出以下要求:

“如果仅考虑感兴趣的中频，则可以找到使用傅立叶变换过滤掉较高和较低频率的LTI系统的频率响应”

我需要知道是否有某种方法可以在MATLAB中自动完成。我不确定如何创建这个带通滤波器，然后手动找到它的频率响应，所以我想知道MATLAB是否会以某种方式自动完成它，因为我可能对此太过思索了。

Answer 1

通过“中频”，我假设您想滤除信号，以便仅存在159 Hz的分量。这意味着输出结果应仅包含一个159 Hz的正弦波。首先，我想提一下傅立叶变换显示了信号的频率分解。您可以将信号视为许多不同频率的正弦波的总和，并且它们可能因不同的相移而异相。对于每个频率的正弦波，都有一个幅度分量和一个相位分量。幅度告诉您在该频率下正弦波有多大，相位告诉您正弦波正经历多少延迟。

现在，在计算FFT时，它将执行Cooley-Tukey算法，这是一种非常有效的傅立叶变换计算方法。它的计算方式是信号的上半部分显示0 <= f < fn的频谱，信号的下半部分显示-fn <= f < 0的频谱。 fn是所谓的Nyquist frequency，它是采样频率的一半。请注意，频谱的前半部分与后半部分相比，其排他性。我们在上半部分不包括fn，但在下半部分包括-fn。

这仅适用于实际信号，这是将三个正弦波加在一起的情况。为了使这一点更有意义，可以使用 fftshift 以便重新组织频谱，使0 Hz / DC频率出现在中间，这将允许您绘制频谱，使其位于-fn <= f < fn之间。

同样，绘制信号的方式是归一化，这意味着您看到的频率实际上在-1 <= f <= 1之间。要将其转换为实际频率，请使用以下关系:

freq = i * fs / N

i是您想要的FFT上的bin号或点，在您的情况下为0到500。请记住，信号的前半部分代表从0到fn的频率分布，从0到500的频率分布。信号就是您所需要的。只需将i替换为-i即可找到负频率。 fs是采样频率，N是FFT的大小，在您的情况下，它是信号的总长度1001。要生成与正确点相对应的正确频率，可以使用 linspace 并在N+1之间生成-fn点和fn以确保间距正确，但是由于不在的正端包含fn，因此我们将其从范围的末尾删除。

另外，要绘制信号的幅度和相位，请分别使用 abs 和 angle 。

因此，尝试以这种方式进行绘制，并且也要关注幅度和相位。仅绘制光谱图是不确定的，因为通常存在实部和虚部。

%// Your code
fs = 1000;
t = (0 : 1/fs : 1)';
f = [ 86, 159, 392 ];
x = sum( cos(2*pi * t * f), 2 );
y = fft(x);

%// New code
%// Shift spectrum
ys = fftshift(y);
N = numel(x); %// Determine number of points
mag = abs(ys); %// Magnitude
pha = angle(ys); %// Phase

%// Generate frequencies
freq = linspace(-fs/2, fs/2, N+1);
freq(end) = [];

%// Draw stuff
figure;
subplot(2,1,1);
plot(freq, mag);
xlabel('Frequency (Hz)');
ylabel('Magnitude');

subplot(2,1,2);
plot(freq, pha);
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');

这是我得到的:



如您所见，存在三个与三个正弦波分量相对应的尖峰。您需要中间一个，即159 Hz。要创建带通滤波器，您希望滤除除+/- 159 Hz以外的所有分量。如果要自动执行此操作，则会找到最接近+/- 159 Hz的分档位置，然后展开围绕这两个点的邻域，并确保在将其余组件清零时保持不变。

由于您具有精确的正弦曲线，因此在这方面使用带通滤波器是完全可以接受的。通常，由于振铃和混叠效应而不会执行此操作，因为以这种方式从带通滤波器截取的截止频率会在时域中引入不必要的混叠效应。参见Wikipedia article on aliasing for more details。

因此，要弄清楚我们需要滤除的位置，请尝试使用 min 并找到159 Hz以上在上面生成的频率之间的绝对差-特别是找到位置。一旦找到+159 Hz和-159 Hz的频率，就在这些点周围扩展邻域，确保在频谱中其余点都设置为0的情况下保持它们不变。

[~,min_pt_pos] = min(abs(freq - f(2))); %// Find location where +159 Hz is located
[~,min_pt_neg] = min(abs(freq + f(2))); %// Find location of where -159 Hz is

%// Neighbourhood size
ns = 100; %// Play with this parameter

%// Filtered signal
yfilt = zeros(1,numel(y));

%// Extract out the positive and negative frequencies centered at
%// 159 Hz
yfilt(min_pt_pos-ns/2 : min_pt_pos+ns/2) = ys(min_pt_pos-ns/2 : min_pt_pos+ns/2);
yfilt(min_pt_neg-ns/2 : min_pt_neg+ns/2) = ys(min_pt_neg-ns/2 : min_pt_neg+ns/2);

yfilt现在包含经过滤波的信号，该信号除去了159 Hz之外的所有分量。如果要显示此滤波信号的幅度和相位，我们可以这样做:

mag2 = abs(yfilt);
pha2 = phase(yfilt);
figure;
subplot(2,1,1);
plot(freq, mag2);
xlabel('Frequency (Hz)');
ylabel('Magnitude');

subplot(2,1,2);
plot(freq, pha2);
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');

这是我们得到的:



正如我们所期望的，只有一个很强的频率可以分解该信号，即159 Hz。现在要重构此信号，您​​必须撤消我们所做的居中，并且必须在此滤波后的结果上使用 ifftshift ，然后通过 ifft 进行逆。您可能还会得到小的残留虚数分量，因此在输出结果上使用 real 是一个好主意。

out = real(ifft(ifftshift(yfilt)));

如果绘制此图，我们得到:

plot(t, out);
xlabel('Time (seconds)');
ylabel('Height');

如您所见，有一个具有一个频率的正弦波，频率为159 Hz。不过不要介意振幅。这完全是由于您选择绘制信号的点数不精确，因此某些时间点可能与信号的真实峰值不完全一致。请记住，如果存在一个以上的正弦曲线，那么所有正弦曲线的峰都会相遇，并且您会得到一个更高的峰，而不是单个正弦曲线提供的峰。因为振幅在-1到1之间徘徊，所以可以确定只有一个正弦波存在。如果要选择更精细的粒度步长，可以避免这种情况，因此可以在FFT中选择更多的点。

希望这足以让您入门。祝好运!