当信号通过FIR滤波器的输出,我们能发现输出数据的前部经常为0值,这可能是由于滤波器的延迟造成的,那么对滤波器的延迟量能否修正?一旦修正后输出数据就变短了,又如何处理呢?
这是在通过FIR滤波器时经常遇到的问题。实际上通过FIR滤波器输出的延迟量能修正,而且输出数据也不会变短,这就是下面要讨论的问题。
由于主要是使用第1类FIR滤波器,可知第1类滤波器的相位角为
可看出延迟量τg是(N-1)/2(单位:样点),所以在时间域上消除延迟的(N-1)/2个样点,就
能把FIR滤波器输出数据校正。
案例、 有一个信号由两个正弦信号组成,这两正弦信号的频率分别为5Hz和20Hz,信号幅值都为1,初始相位分别为45°和60°,对信号采样频率为100Hz,信号长为100。在该信号中要消除20Hz,需要设计低通滤波器,通带和阻带为fp=10Hz,fs=15Hz,通带波纹和阻带衰减为Rp=3dB,As=60dB。对滤波器的输出要消除延迟的影响。
由于阻带衰减为60dB,故选择布莱克曼窗函数。程序如下:
- clear all; clc; close all;
-
- % 信号构成
- N=100; % 数据长
- Fs=100; % 采样频率
- Fs2=Fs/2; % 奈奎斯特频率
- f1=5; f2=20; % 两正弦信号频率
- phy1=pi/4; phy2=pi/3; % 两正弦信号初始相位角
- t=(0:N-1)/Fs; % 时间刻度
- x=cos(2*pi*f1*t+phy1)+cos(2*pi*f2*t+phy2);% 构成输入信号
- % FIR低通滤波器设计
- fp=10; fs=15; % 通带和阻带频率
- Rp=3; As=60; % 通带波纹和阻带衰减
- wp=fp*pi/Fs2; ws=fs*pi/Fs2; % 归一化角频率
- deltaw= ws - wp; % 过渡带宽Δω的计算
- M = ceil(11*pi/ deltaw); % 按布莱克曼窗计算滤波器阶数M
- M = M + mod(M,2); % 保证滤波器系数长M+1为奇数
- wind = (blackman(M+1))'; % 布莱克曼窗函数
- Wn=(fp+fs)/Fs; % 计算截止频率
- b=fir1(M,Wn,wind); % 用fir1函数设计FIR第1类滤波器
- [db,mag,phs,gdy,w]=freqz_m(b,1); % 计算滤波器响应
- % 信号滤波
- y1=filter(b,1,x); % 用filter函数进行滤波
- y2=conv(b,x); % 用conv函数进行滤波
- y21=y2(M/2+1:M/2+N); % 取y2中的有效部分
- N2=length(y2); % 求y2长度
- t2=(0:N2-1)/Fs; % y2的时间刻度
- % 作图
- figure(1)
- pos = get(gcf,'Position');
- set(gcf,'Position',[pos(1), pos(2)-100,pos(3),(pos(4)-160)]);
- plot(w/pi*Fs2,db,'k')
- title('低通滤波器的幅值响应');
- grid; axis([0 Fs2 -150 10]);
- xlabel('频率/Hz'); ylabel('幅值/dB')
- set(gca,'XTickMode','manual','XTick',[0,10,15,20,50])
- set(gca,'YTickMode','manual','YTick',[-100,-60,0])
- set(gcf,'color','w');
- figure(2)
- subplot 311; plot(t,y1,'k');
- ylim([-1.1 1.1]);
- title('(a)通过filter函数FIR滤波器的输出');
- xlabel('时间/s'); ylabel('幅值')
- subplot 312; plot(t2,y2,'k');
- ylim([-1.1 1.1]);
- title('(b)通过conv函数FIR滤波器的输出');
- xlabel('时间/s'); ylabel('幅值')
- subplot 313; plot(t,y21,'k');
- ylim([-1.1 1.1]);
- title('(c)对(b)的输出进行修正');
- xlabel('时间/s'); ylabel('幅值')
- set(gcf,'color','w');
运行结果如下:
分析:
对FIR滤波器的滤波过程可以用filter函数也可用conv函数,但这两个函数滤波的结果是完全不同的。从基础理论部分,我们知道第1类FIR滤波器的输出将会产生延迟(M-1)/2,M是滤波器的长度(或为N/2,N是第1类FIR滤波器的阶数)。当用函数conv来滤波时滤波器输出:y=conv(x,b),y序列长度是L+M-1,其中L为数据x的长度,M是滤波器的长度。由于x只有L长,在输出y中只有L长的数据是滤波器输出的有效部分,而我们已知滤波器的输出有(M-1)/2个样点的延迟,所以输出y中从(M-1)/2+1开始取数值,取L长的样点,构成y的有效部分。
filter函数输出只有L长,它把(M-1)/2个由延迟产生的无效样点也包含在内。从图(a)中看到,filter函数输出y1的波形中一大半是接近0值,这是因为FIR滤波器是N=112阶,延迟量N/2=56个样点,在采样频率100Hz时占大于0.5s的时间,所以y1的波形中大于一半是接近0值;如果把前面的时延部分删除,则输出数据变短了。
通过conv函数的输出y2,波形图显示在图(b)中,由于y2长L+M-1,所以时间上大于2s。取了y21=y2((M-1)/2+1:(M-1)/2+L)后得y21,显示在图(c)中,这就是我们所需要的滤波后输出的正弦信号。
但是得到的滤波输出y21也不是完全没有问题,在初始部分由于瞬态现象的存在,与原始输入信号之间还存在一些偏差。
参考文献:MATLAB数字信号处理85个实用案例精讲——入门到进阶;宋知用(编著)