我有一个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);
并得到了相同的结果。拜托,有人可以向我解释我在做什么错吗?谢谢!
Your passband does not cover any of the three signal frequency components. This can be seen directly on the graph (the second figure, containing the impulse response, has too fast variations compared with the signal). Or it can be seen from the numbers 0.5760
and 0.3655
in
h = (sin(0.5760*pi*(n-M))-sin(0.3665*pi*(n-M)))./pi./(n-M);
Why did you choose those numbers? The normalized frecuencies of the signal are [2 5 8]/60
, that is, 0.0333 0.0833 0.1333
. They are all below the passband [.3665 .5760]
. As a result, the filter greatly attenuates the three components of the input signal.
Say you want to isolate the center frequency component (f = 5000
Hz, or f/fs = 0.08333
normalized frequency). You need a bandpass filter than lets that frequency through, and rejects the others. So you would use for example a normalized passband [.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
A second problem with your code is that it gives two errors because you index h
with 0
. To solve this, change
n = 0:(N-1);
to
n = 1:N;
and
y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j);
to
y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j+1);
So, the modified code to isolate the central frequency component is:
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)];
The result is as follows.
As can be seen, only the central frequency component is present in the output signal.
It's also observed that the envelope of the output signal is not constant. That's because the duration of the input signal is comparable to the filter lenght. That is, you are seeing the transient response of the filter. Note that the envelop rise times is approximately the length of the filter's impulse response h
.
在这里注意到一个有趣的折衷是很有趣的(信号处理充满了折衷)。为了使瞬变更短,可以使用更宽的通带(滤波器长度与通带成反比)。但是,这时滤波器的选择性较低,也就是说,它在不想要的频率处的衰减较小。例如,通过passband查看结果[.04 .12]
:
正如预期的那样,瞬变现在变短了,但是所需的频率却变得不那么纯净了:您可以看到由其他频率的剩余部分引起的“调制”。
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句