Matlab：与带通滤波器的卷积不会消除不需要的频率

Question

我有一个 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)]);andyyy = conv(h,x);并得到了相同的结果。请问有人可以解释一下我做错了什么吗？谢谢你！

Answer 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 = 5000Hz，或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]：

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