找出带通滤波器的频率响应
Finding the frequency response of a bandpass filter
我有这个 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);
（为了方便起见，我将省略轴和绘图标题的代码）。这是我迄今为止为报告生成的两个数字。
报告的下一部分问我这个：
"If we consider only the middle frequency of interest find the frequency response of a LTI system that filters out the higher and lower frequencies using the fourier transform"
我想知道是否有某种方法可以在 MATLAB 中自动完成。我不确定如何创建这个带通滤波器，然后手动找到它的频率响应，所以我想知道 MATLAB 是否有某种方法可以自动完成它，因为我可能想多了。
通过 "middle" 频率，我假设您想过滤掉信号，以便只存在 159 Hz 的分量。这意味着输出结果应该只包含一个 159 Hz 的正弦波。首先，我想提一下，傅立叶变换显示了信号的频率分解。您可以将信号视为许多不同频率的正弦波的总和，它们可能因不同的相移而异相。对于一个频率的每个正弦波，都有一个 幅度 分量和一个 相位 分量。幅度告诉您正弦波在此频率下有多大，相位告诉您正弦波经历了多少延迟。
现在，在计算 FFT 时，这会执行 Cooley-Tukey algorithm, which is a very efficient way of computing the Fourier Transform. It's computed in such a way where the first half of the signal shows the spectrum from 0 <= f < fn and the second half of the signal shows the spectrum from -fn <= f < 0. fn is what is known as the Nyquist frequency，这是采样频率的一半。请注意与后半部分相比，前半部分光谱范围内的排他性。我们不在上半场包括fn，但我们在下半场包括-fn。
这仅适用于 真实 信号，这就是您将三个正弦信号加在一起的情况。为了使这更有意义，您可以使用 fftshift 重新组织频谱，使 0 Hz / DC 频率出现在中间，这将允许您绘制频谱，使其在 -fn <= f < fn.
此外，信号的绘制方式是 标准化 ，这意味着您看到的频率实际上在 -1 <= f <= 1 之间。要将其转换为实际频率，请使用以下关系：
freq = i * fs / N
i 是您想要的二进制数或 FFT 上的点，在您的情况下是从 0 到 500。请记住，信号的前半部分代表从 0 到fn 并且信号中从 0 到 500 是您所需要的。只需将 i 替换为 -i 即可找到负频率。 fs 是采样频率，N 是 FFT 的大小，在您的情况下是信号的总长度，1001。要生成与正确点对应的正确频率，您可以使用 linspace 并在 -fn 和 fn 之间生成 N+1 点以确保间距正确，但因为我们 不 包括fn 在正端，我们将其从范围末尾删除。
此外，要绘制信号的幅度和相位，请分别使用 abs and 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 的 bin 位置，然后扩展围绕这两个点的邻域并确保它们未被触及，同时将其余组件归零。
因为你有精确的正弦波，所以在这方面使用带通滤波器是完全可以接受的。通常，由于振铃和混叠效应，您不会这样做，因为带通滤波器截止的锐度会在时域中引入不需要的混叠效应。见 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 on this filtered result, then do the inverse via ifft. You may also get small residual imaginary components, so it's a good idea to use real。
out = real(ifft(ifftshift(yfilt)));
如果我们绘制这个，我们得到：
plot(t, out);
xlabel('Time (seconds)');
ylabel('Height');
如您所见，有一个频率为 159 Hz 的正弦波。不过不要介意振幅。这仅仅是由于您选择绘制信号的点数不精确，因此某些时间点可能与信号的真实峰值不完全一致。请记住，如果存在多个正弦波，那么所有正弦波的峰值都会在一个点相遇，您将获得更高的峰值，而不是单个正弦波提供的峰值 1。由于振幅在 -1 到 1 之间徘徊，您可以确定只存在一个正弦波。如果您要选择更细粒度的步长，并因此在 FFT 中选择更多的点，则可以避免这种情况。
希望这足以让您入门。祝你好运！
我有这个 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);
（为了方便起见，我将省略轴和绘图标题的代码）。这是我迄今为止为报告生成的两个数字。
报告的下一部分问我这个：
"If we consider only the middle frequency of interest find the frequency response of a LTI system that filters out the higher and lower frequencies using the fourier transform"
我想知道是否有某种方法可以在 MATLAB 中自动完成。我不确定如何创建这个带通滤波器，然后手动找到它的频率响应，所以我想知道 MATLAB 是否有某种方法可以自动完成它，因为我可能想多了。
通过 "middle" 频率，我假设您想过滤掉信号，以便只存在 159 Hz 的分量。这意味着输出结果应该只包含一个 159 Hz 的正弦波。首先，我想提一下，傅立叶变换显示了信号的频率分解。您可以将信号视为许多不同频率的正弦波的总和，它们可能因不同的相移而异相。对于一个频率的每个正弦波，都有一个 幅度 分量和一个 相位 分量。幅度告诉您正弦波在此频率下有多大，相位告诉您正弦波经历了多少延迟。
现在，在计算 FFT 时，这会执行 Cooley-Tukey algorithm, which is a very efficient way of computing the Fourier Transform. It's computed in such a way where the first half of the signal shows the spectrum from 0 <= f < fn and the second half of the signal shows the spectrum from -fn <= f < 0. fn is what is known as the Nyquist frequency，这是采样频率的一半。请注意与后半部分相比，前半部分光谱范围内的排他性。我们不在上半场包括fn，但我们在下半场包括-fn。
这仅适用于 真实 信号，这就是您将三个正弦信号加在一起的情况。为了使这更有意义，您可以使用 fftshift 重新组织频谱，使 0 Hz / DC 频率出现在中间，这将允许您绘制频谱，使其在 -fn <= f < fn.
此外，信号的绘制方式是 标准化 ，这意味着您看到的频率实际上在 -1 <= f <= 1 之间。要将其转换为实际频率，请使用以下关系：
freq = i * fs / N
i 是您想要的二进制数或 FFT 上的点，在您的情况下是从 0 到 500。请记住，信号的前半部分代表从 0 到fn 并且信号中从 0 到 500 是您所需要的。只需将 i 替换为 -i 即可找到负频率。 fs 是采样频率，N 是 FFT 的大小，在您的情况下是信号的总长度，1001。要生成与正确点对应的正确频率，您可以使用 linspace 并在 -fn 和 fn 之间生成 N+1 点以确保间距正确，但因为我们 不 包括fn 在正端，我们将其从范围末尾删除。
此外，要绘制信号的幅度和相位，请分别使用 abs and 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 的 bin 位置，然后扩展围绕这两个点的邻域并确保它们未被触及，同时将其余组件归零。
因为你有精确的正弦波，所以在这方面使用带通滤波器是完全可以接受的。通常，由于振铃和混叠效应，您不会这样做，因为带通滤波器截止的锐度会在时域中引入不需要的混叠效应。见 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 on this filtered result, then do the inverse via ifft. You may also get small residual imaginary components, so it's a good idea to use real。
out = real(ifft(ifftshift(yfilt)));
如果我们绘制这个，我们得到：
plot(t, out);
xlabel('Time (seconds)');
ylabel('Height');
如您所见，有一个频率为 159 Hz 的正弦波。不过不要介意振幅。这仅仅是由于您选择绘制信号的点数不精确，因此某些时间点可能与信号的真实峰值不完全一致。请记住，如果存在多个正弦波，那么所有正弦波的峰值都会在一个点相遇，您将获得更高的峰值，而不是单个正弦波提供的峰值 1。由于振幅在 -1 到 1 之间徘徊，您可以确定只存在一个正弦波。如果您要选择更细粒度的步长，并因此在 FFT 中选择更多的点，则可以避免这种情况。
希望这足以让您入门。祝你好运！