数字滤波实用介绍
这个例子展示了如何设计、分析和对数据应用数字滤波器。它将帮助您回答以下问题:如何补偿过滤器带来的延迟?如何避免信号失真?,How do I remove unwanted content from my signal?, How do I differentiate my signal?, and How do I integrate my signal?
滤波器可用于以所需的方式塑造信号频谱或执行微分和积分等数学运算。在接下来的内容中,您将学习一些实用的概念,以便在需要时轻松使用过滤器。
这个例子关注的是数字滤波器的应用而不是它们的设计。如果您想了解有关如何设计数字滤波器的更多信息,请参阅数字滤波器设计实用介绍的例子。
滤波引入的延迟补偿
数字滤波器在信号中引入延迟。根据滤波器的特性,延迟可以在所有频率上保持恒定,也可以随频率变化。延迟的类型决定了您必须采取的补偿措施。的grpdelay
函数允许您将滤波器延迟视为频率的函数。查看这个函数的输出可以让您确定滤波器的延迟是恒定的,还是随频率变化(换句话说,它是否依赖于频率)。
滤波器延迟在所有频率上都是恒定的,可以通过及时改变信号来很容易地补偿。FIR滤波器通常具有恒定的延迟。另一方面,随频率变化的延迟会引起相位失真,并能显著改变信号波形。对频率相关的延迟进行补偿并不像对常数延迟情况那样微不足道。IIR滤波器引入频率相关的延迟。
补偿常量滤波器延迟
如前所述,您可以测量滤波器的延迟组,以验证它是频率的常数函数。您可以使用grpdelay
函数来测量滤波器延迟D,并通过对输入信号附加D个零并及时移动输出信号D个样本来补偿这种延迟。
考虑一个有噪声的心电图信号,您希望对其进行滤波以去除75 Hz以上的高频噪声。您希望应用FIR低通滤波器并补偿滤波器延迟,以便噪声信号和滤波信号正确对齐,并可以相互绘制以进行比较。
Fs = 500;%采样率(Hz)N = 500;%信号样本数rng默认的;x = ecg(N)'+0.25*randn(N,1);噪声波形t = (0:N-1)/Fs;%时间向量
设计一个70阶低通FIR滤波器,截止频率为75hz。
Fnorm = 75/(Fs/2);归一化频率%Df = designfilt(“lowpassfir”,“FilterOrder”, 70,“CutoffFrequency”, Fnorm);
绘制滤波器的组延迟,以验证它在所有频率上都是常数,表明滤波器是线性相位。使用组延迟来测量滤波器的延迟。
grpdelay (df、2048 Fs)%绘图组延迟
D = mean(grpdelay(df))样本中的过滤延迟
D = 35
在滤波之前,在输入数据向量x的末尾附加D个零。这确保了所有有用的样本都从滤波器中被清除,并且输入信号和延迟补偿的输出信号具有相同的长度。对数据进行滤波,通过对输出信号进行D个采样的移位来补偿延迟。最后一步有效地删除了暂态过滤器。
Y = filter(df,[x;0 (D, 1)]);在输入数据后附加D个零y = y(D+1:结束);%移位数据以补偿延迟情节(t t, x,, y,“线宽”, 1.5)标题(过滤后的波形的)包含(“时间(s)”)传说(“原始噪声信号”,“过滤信号”网格)在轴紧
补偿频率相关的延迟
频率相关的延迟导致信号的相位失真。对这种类型的延迟进行补偿并不像对恒定延迟情况那样微不足道。方法实现零相位滤波,如果应用程序允许脱机处理,则可以消除频率相关的延迟filtfilt
函数。filtfilt
通过正向和反向处理输入数据来执行零相位滤波。主要效果是获得零相位失真,即使用等效滤波器过滤数据,该滤波器具有恒定的0个样本延迟。其他的影响是,你得到的滤波器传递函数等于原始滤波器传递函数的模的平方,并且滤波器的阶数是原始滤波器阶数的两倍。
考虑上一节中定义的ECG信号。用或不带延迟补偿对该信号进行滤波。设计了一个截止频率为75 Hz的7阶低通IIR椭圆滤波器。
Fnorm = 75/(Fs/2);归一化频率%Df = designfilt(“lowpassiir”,...“PassbandFrequency”Fnorm,...“FilterOrder”7...“PassbandRipple”, 1...“StopbandAttenuation”、60);
绘制过滤器的组延迟。群时延随频率变化,表明滤波器时延与频率相关。
grpdelay (df, 2048,“一半”Fs)
过滤数据并查看每个过滤器实现对时间信号的影响。零相滤波有效地消除了滤波延时。
Y1 = filter(df,x);%非线性相位滤波器-无延迟补偿Y2 = filtfilt(df,x);%零阶段实现-延迟补偿情节(t, x)在情节(t, y1,“r”,“线宽”1.5)情节(t, y2,“线宽”, 1.5)标题(过滤后的波形的)包含(“时间(s)”)传说(原始信号的,“非线性相位IIR输出”,...“零相IIR输出”) xlim([0.25 0.55]在
如果您的应用程序允许非因果正向/反向过滤操作,并且允许将过滤器响应更改为原始响应的平方,那么零相位过滤是一个很好的工具。
引入常数延迟的滤波器是线性相位滤波器。引入频率相关延迟的滤波器是非线性相位滤波器。
从信号中移除不需要的频谱内容
滤波器通常用于从信号中去除不需要的光谱内容。你可以从各种各样的过滤器中进行选择。当您想要删除高频内容时,您可以选择低通滤波器,当您想要删除低频内容时,您可以选择高通滤波器。您还可以选择带通滤波器来去除低频和高频内容,同时保留中间频带的频率。当您想要去除给定频带上的频率时,可以选择带阻滤波器。
考虑一个有电线嗡嗡声和白噪声的音频信号。电源线嗡嗡声是由60hz的音引起的。白噪声是一种存在于所有音频带宽中的信号。
加载音频信号。指定44.1 kHz的采样率。
Fs = 44100;Y = audieread (“noisymusic.wav”);
绘制信号的功率谱。红色三角形标记表示干扰音频信号的60hz强音调。
[P,F] = pwelch(y,ones(8192,1),8192/2,8192,F,“权力”);helperFilterIntroductionPlot1(F,P,[60 60],[-9.365 -9.365],...{“原始信号功率谱”,“60赫兹音”})
您可以首先使用低通滤波器去除尽可能多的白噪声频谱内容。滤波器的通频带应设置为在降低噪声和由于高频内容丢失而导致的音频退化之间提供良好平衡的值。在去除60hz嗡嗡声之前应用低通滤波器非常方便,因为您将能够对带限信号进行降采样。低速率信号将允许您设计一个更清晰和更窄的60 Hz带阻滤波器与更小的滤波器阶数。
设计一个低通滤波器,通带频率为1khz,阻带频率为1.4 kHz。选择最小订货量的设计。
Fp = 1e3;%通频带频率,单位为HzFst = 1.4e3;%停止频带频率(Hz)Ap = 1;%通带纹波,单位为dBAst = 95;%阻带衰减,以dB为单位Df = designfilt(“lowpassfir”,“PassbandFrequency”《外交政策》,...“StopbandFrequency”置,“PassbandRipple”据美联社,,...“StopbandAttenuation”Ast,“SampleRate”Fs);
分析过滤器响应。
fvtool (df,“Fs”Fs,“FrequencyScale”,“日志”,...“FrequencyRange”,'指定频率矢量',“FrequencyVector”F)
过滤数据并补偿延迟。
D = mean(grpdelay(df));%滤波器延迟Ylp = filter(df,[y;0 (D, 1)]);ylp = ylp(D+1:end);
看看低通滤波信号的频谱。1400hz以上的频率内容已被删除。
[Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192, f,“权力”);Flp helperFilterIntroductionPlot1 (F P、Plp...{原始信号的,低通滤波信号})
从上面的功率谱图可以看出,低通滤波信号的最大不可忽略频率含量在1400hz。由抽样定理,抽样率为 Hz足以正确地表示信号,然而,您使用的是44100 Hz的采样率,这是一种浪费,因为您将需要处理更多的样本。通过减少需要处理的样本数量,您可以对信号进行低采样以降低采样率并减少计算负载。较低的采样率也将允许您设计一个更清晰和更窄的带阻滤波器,以较小的滤波器阶数去除60 Hz噪声。
对低通滤波信号进行10倍的下采样,得到Fs/10 = 4.41 kHz的采样率。绘制下采样前后信号的频谱图。
Fs = Fs/10;Yds = downsample(ylp,10);[Pds,Fds] = pwelch(yds,ones(8192,1),8192/2,8192,Fs,“权力”);helperFilterIntroductionPlot1 (F P Fds, Pds,...{“信号以44100赫兹采样”,下采样信号,Fs = 4410 Hz})
现在使用IIR带阻滤波器去除60赫兹的音调。让止带的宽度为4hz,居中于60hz。选择IIR滤波器,以实现尖锐的频率缺口,小的通带纹波,和相对较低的阶数。
Df = designfilt(“bandstopiir”,“PassbandFrequency1”现年55岁的...“StopbandFrequency1”58岁的“StopbandFrequency2”, 62,...“PassbandFrequency2”, 65,“PassbandRipple1”, 1...“StopbandAttenuation”现年60岁的“PassbandRipple2”, 1...“SampleRate”Fs,“DesignMethod”,“ellip”);
分析震级响应。
fvtool (df,“Fs”Fs,“FrequencyScale”,“日志”,...“FrequencyRange”,'指定频率矢量',“FrequencyVector”Fds (Fds > F (2)))
进行零相位滤波,避免相位失真。使用filtfilt
函数来处理数据。
Ybs = filtfilt(df,yds);
最后,对信号进行上采样,使其恢复到与音频声卡兼容的原始音频采样率44.1 kHz。
Yf = interp(ybs,10);Fs = Fs*10;
最后看一下原始信号和经过处理的信号的频谱。高频噪声地板和60赫兹的音调已被滤波器衰减。
[Pfinal, final] = pwelch(yf,ones(8192,1),8192/2,8192, f,“权力”);helperFilterIntroductionPlot1 (F P Ffinal Pfinal,...{原始信号的,“最终过滤信号”})
如果你的电脑有声卡,你可以听处理前后的信号。如上所述,最终结果是您已经有效地减弱了音频文件中的60 Hz嗡嗡声和高频噪声。
要播放原始信号,取消接下来两行的注释% hplayer = audioplayer(y, Fs);% (hplayer)玩要播放降噪信号,取消下面两行注释% hplayer = audioplayer(yf, Fs);% (hplayer);
对信号求导
MATLABdiff
函数区分信号的缺点是可以潜在地增加输出的噪声级。更好的选择是使用微分器滤波器,在感兴趣的波段作为微分器,并在所有其他频率作为衰减器,有效地去除高频噪声。
作为一个例子,分析地震时建筑物地板的位移速度。地震条件下,位移或漂移测量被记录在一个三层测试结构的一楼,并保存在地震流中。垫文件。数据向量的长度为10e3,采样率为1 kHz,测量单位为cm。
微分位移数据,以获得地震期间建筑物地板的速度和加速度的估计。比较使用diff和FIR微分器滤波器的结果。
负载quakedrift.matFs = 1000;%抽样率dt = 1/Fs;%时间差T =(0:长度(漂移)-1)*dt;%时间向量
设计一个50阶微分器滤波器,其通频带频率为100 Hz,这是发现大部分信号能量的带宽。设置滤波器的阻带频率为120hz。
Df = designfilt(“differentiatorfir”,“FilterOrder”, 50岁,...“PassbandFrequency”, 100,“StopbandFrequency”, 120,...“SampleRate”Fs);
的diff
函数可以看作是一个有响应的一阶FIR滤波器
.利用FVTool比较了50阶微分器FIR滤波器的幅值响应和50阶微分器FIR滤波器的幅值响应diff
函数。显然,这两种响应在通带区域(从0到100 Hz)是等效的。然而,在阻带区域,50阶滤波器衰减分量,而差分响应放大分量。这有效地提高了高频噪声的水平。
HFVT = fvtool(df,[1 -1],1,“MagnitudeDisplay”,“零”,“Fs”Fs);传奇(hfvt“50阶FIR微分器”,diff函数的响应);
区分使用diff
函数。加零来补偿由于差异操作导致的缺失样本。
V1 = diff(漂移)/dt;A1 = diff(v1)/dt;V1 = [0;v1);A1 = [0;0;a1];
使用50阶FIR滤波器进行微分,并补偿延迟。
D = mean(grpdelay(df));%滤波器延迟V2 = filter(df,[漂移;0 (D, 1)]);v2 = v2(D+1:结束);A2 = filter(df,[v2;0 (D, 1)]);a2 = a2(D+1:结束);V2 = V2 /dt;A2 = A2 /dt^2;
绘制楼层位移的几个数据点。还绘制了用diff和50阶FIR滤波器计算的速度和加速度的一些数据点。注意噪声是如何在速度估计中被略微放大的,而在加速度估计中被大大放大的diff
.
helperFilterIntroductionPlot2 (t,漂移,v1、v2, a1, a2)
对信号进行积分
漏积分滤波器是一种具有传递函数的全极滤波器
在哪里
是一个常数,必须小于1,以确保过滤器的稳定性。这并不奇怪
趋近于1,漏积分器趋近于的逆diff
传递函数。将泄漏积分器应用于前一节中获得的加速度和速度估计,分别得到速度和漂移。使用得到的估计的
diff
因为它们噪音更大。
使用漏积分器 .绘制泄漏积分器滤波器的幅值响应。该滤波器作为低通滤波器有效地消除高频噪声。
fvtool (1 [1 -.999]“Fs”Fs)
用漏积分器对速度和加速度进行滤波。乘以时间差。
V_original = v1;A_original = a1;D_leakyint = filter(1,[1 -0.999],v_original);V_leakyint = filter(1,[1 -0.999],a_original);D_leakyint = D_leakyint * dt;V_leakyint = V_leakyint * dt;
绘制位移和速度估计值,并与原始信号进行比较。
helperFilterIntroductionPlot3 (t,漂移,d_leakyint、v_original v_leakyint)
函数也可以对信号进行积分cumsum
而且cumtrapz
功能。结果将与使用漏积分器得到的结果相似。
结论
在本例中,您学习了线性和非线性相位滤波器,并学习了如何补偿每种滤波器类型引入的相位延迟。您还学习了如何应用滤波器从信号中去除不需要的频率成分,以及如何在使用低通滤波器限制信号带宽后对其进行低采样。最后,您学习了如何使用数字滤波器设计来区分和集成信号。在整个示例中,您还学习了如何使用分析工具查看过滤器的响应和分组延迟。
有关滤波器应用程序的更多信息,请参阅信号处理工具箱™文档。有关如何设计数字滤波器的更多信息,请参阅数字滤波器设计实用介绍的例子。
参考文献
普罗基斯,J. G.和D. G.马诺拉基斯。数字信号处理:原理,算法和应用.恩格尔伍德悬崖,新泽西州:Prentice-Hall, 1996年。
奥法尼迪斯,s.j.。信号处理导论.恩格尔伍德悬崖,新泽西州:Prentice-Hall, 1996年。
附录
本例中使用了以下helper函数:
另请参阅
带通
|bandstop
|designfilt
|fftfilt
|过滤器
|filtfilt
|grpdelay
|高通滤波
|低通滤波器