Matlab：带通滤波器的卷积不会减少不需要的频率答案

【问题标题】：Matlab: convolution with bandpass filter does not cut the unwanted frequenciesMatlab：带通滤波器的卷积不会减少不需要的频率
【发布时间】：2015-12-11 11:40:51
【问题描述】：

我有一个 6ms 长的信号，其中三个频率分量以 60kHz 采样：

fs = 60000;
T = 0.006;
t = 0:1/fs:T;
x = 0.3*sin(2*pi*2000*t) + sin(2*pi*5000*t) + 0.4*sin(2*pi*8000*t);

我有一个带通滤波器，其脉冲响应是两个 sinc 函数之间的差异：

M = 151;
N = 303;
n = 0:(N-1);
h = (sin(0.5760*pi*(n-M))-sin(0.3665*pi*(n-M)))./pi./(n-M);
h(n==M) = 0.2094;

我设计了一个将输入与过滤器进行卷积的函数：

function y = fir_filter(h,x)
y = zeros(1,length(x)+length(h)-1);
for i = 1:length(x)
for j = 1:length(h)
    y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j);
end
end

然后应用过滤器：

y = fir_filter(h,x);

这产生了奇怪的结果：

figure(21)
ax1 = subplot(311);
plot(x);
title('Input Signal');
ax2 = subplot(312);
plot(h);
title('FIR');
ax3 = subplot(313);
plot(y);
title('Output Signal');
linkaxes([ax1,ax2,ax3],'x')
ax2.XLim = [0,length(y)];

由于滤波器是带通的，因此预计只有一个频率分量可以保留。

我尝试使用yy = filter(h,1,[x,zeros(1,length(h)-1)]); 和yyy = conv(h,x); 并得到相同的结果。拜托，谁能解释我做错了什么？谢谢！

【问题讨论】：

剧情是什么？频率响应？你的问题是fft镜像？阅读dsp.stackexchange.com/questions/4825/why-is-the-fft-mirrored或phys.nsu.ru/cherk/fft.pdf
问题是我的输入与 FIR 滤波器的卷积产生了错误的结果。结果应该接近纯正弦曲线
可以添加代码来生成绘图吗？
好的，我编辑了帖子
您的代码对我不起作用，因为fir_filter 中的`h(length(h)-j)` 在j=length(h) 时访问h(0)。只是出于好奇，为什么不conv？

标签： matlab filtering signal-processing

【解决方案1】：

您的通带不涵盖三个信号频率分量中的任何一个。这可以直接在图表上看到（第二个图包含脉冲响应，与信号相比变化太快）。或者从数字0.5760和0.3655中可以看出

h = (sin(0.5760*pi*(n-M))-sin(0.3665*pi*(n-M)))./pi./(n-M);

您为什么选择这些数字？信号的归一化频率为[2 5 8]/60，即0.0333 0.0833 0.1333。它们都低于通带[.3665 .5760]。结果，滤波器大大衰减了输入信号的三个分量。

假设您要隔离中心频率分量（f = 5000 Hz 或 f/fs = 0.08333 归一化频率）。您需要一个带通滤波器，而不是让该频率通过并拒绝其他频率。因此，您将使用例如标准化通带[.06 .1]：

h = (sin(0.06*pi*(n-M))-sin(0.1*pi*(n-M)))./pi./(n-M);
h(n==M) = (h(n==M+1)+h(n==M-1))/2; %// auto adjustment to avoid the 0/0 sample

您的代码的第二个问题是它给出了两个错误，因为您将h 索引为0。为了解决这个问题，改变

n = 0:(N-1);

到

n = 1:N;

和

y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j);

到

y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j+1);

所以，隔离中心频率分量的修改代码是：

fs = 60000;
T = 0.006;
t = 0:1/fs:T;
x = 0.3*sin(2*pi*2000*t) + sin(2*pi*5000*t) + 0.4*sin(2*pi*8000*t); 

M = 151;
N = 303;
n = 1:N; %// modified
h = (sin(0.06*pi*(n-M))-sin(0.1*pi*(n-M)))./pi./(n-M); %// modified
 h(n==M) = (h(n==M+1)+h(n==M-1))/2; %// modified


y = zeros(1,length(x)+length(h)-1);
for i = 1:length(x)
for j = 1:length(h)
    y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j+1); %// modified
end
end

figure(21)
ax1 = subplot(311);
plot(x);
title('Input Signal');
ax2 = subplot(312);
plot(h);
title('FIR');
ax3 = subplot(313);
plot(y);
title('Output Signal');
linkaxes([ax1,ax2,ax3],'x')
ax2.XLim = [0,length(y)];

结果如下。

可以看出，输出信号中只有中心频率分量。

还观察到输出信号的包络不是恒定的。这是因为输入信号的持续时间与滤波器长度相当。也就是说，您看到的是过滤器的瞬态响应。请注意，包络上升时间大约是滤波器脉冲响应h 的长度。

在这里注意到一个有趣的权衡很有趣（信号处理充满了权衡）。为了使瞬态更短，您可以使用更宽的通带（滤波器长度与通带成反比）。但是滤波器的选择性会降低，也就是说，它在不需要的频率处的衰减会更小。例如，使用 passband [.04 .12] 查看结果：

正如预期的那样，瞬态现在更短了，但所需的频率不太纯净：您可以看到由其他频率的剩余部分引起的一些“调制”。

【讨论】：