主要内容

数字滤波实用导论

此示例显示如何为数据设计,分析和应用数字过滤器。它将帮助您回答问题:如何弥补过滤器引入的延迟?,如何避免扭曲我的信号?,如何从​​我的信号中删除不需要的内容?,如何区分我的信号?,我如何整合我的信号?

滤波器可用于以期望的方式塑造信号频谱或执行数学运算,如微分和积分。在接下来的内容中,您将学习一些实用的概念,当您需要时,这些概念将简化过滤器的使用。

这个例子着重于数字滤波器的应用而不是其设计。如果你想了解更多关于如何设计数字滤波器,请参阅数字滤波器设计的实用介绍的例子。

滤波引起的延迟补偿

数字滤波器会在信号中引入延迟。根据滤波器特性,延迟可以在所有频率上保持不变,也可以随频率变化。延迟的类型决定您必须采取的补偿措施grpdelay函数可以让你看到滤波器延迟作为频率的函数。通过查看这个函数的输出,您可以确定滤波器的延迟是常量还是随频率变化(换句话说,它是否与频率有关)。

滤波器延迟在所有频率上都是恒定的,可以很容易地通过在时间上移动信号来补偿。FIR滤波器通常有恒定的延迟。另一方面,随频率变化的延迟会引起相位失真,并显著改变信号波形。对频率相关延迟的补偿不像对常量延迟的补偿那么简单。IIR滤波器引入频率相关延迟。

恒定滤波器延迟的补偿

如前所述,可以测量滤波器的延迟组,以验证它是频率的常数函数。你可以使用grpdelay函数来测量滤波器延迟D,并通过在输入信号上添加d0并通过D采样在时间上移动输出信号来补偿这个延迟。

考虑一个要过滤的嘈杂的心电图信号,以删除75 Hz以上的高频噪声。您想应用FIR低通滤波器并补偿过滤器延迟,以便噪声和滤波的信号正确对齐,可以绘制在彼此的顶部以进行比较。

Fs = 500;Hz中的样品率%N = 500;%信号样本数rng默认的;X = ECG(n)'+ 0.25 * RANDN(N,1);%的波形t = (0: n - 1) / Fs;%的时间向量

设计一个截止频率为75hz的70阶低通FIR滤波器。

fnorm = 75 /(fs / 2);%归一化频率df = designfilt (“lowpassfir”“FilterOrder”,70,'cutfffrequency', Fnorm);

绘制滤波器的组延迟,以验证它是常数跨越所有频率,表明滤波器是线性相位。使用组延迟来测量滤波器的延迟。

grpdelay (df、2048 Fs)%绘图组延迟

图形过滤器可视化工具-组延迟包含uitoolbar、uimenu类型的轴和其他对象。标题组延迟的轴包含类型为line的对象。

D=平均值(GRP延迟(df))采样滤波器延迟%
D = 35

在过滤之前,请在输入数据向量的末尾附加D Zeros,x。这确保了所有有用的样本被冲出过滤器,并且输入信号和延迟补偿输出信号具有相同的长度。通过D样品转换输出信号来筛选数据并补偿延迟。最后一步有效地去除过滤器瞬态。

y =滤波器(df,[x; zeros(d,1)]);%将D个零附加到输入数据y = y (D + 1:结束);%移动数据以补偿延迟绘图(t,x,t,y,“线宽”,1.5)标题(过滤后的波形的)包含(“时间(s)”)传说('原始嘈杂的信号'“过滤信号”网格)

图中包含一个坐标轴。标题为“过滤波形”的轴包含两个类型为line的对象。这些对象分别代表原始噪声信号和滤波信号。

补偿频率相关延迟

频率相关的延迟导致信号的相位失真。这种类型的延迟的补偿不像常量延迟的补偿那么简单。如果您的应用程序允许脱机处理,您可以通过使用filtfilt.函数。filtfilt.通过对输入数据进行正向和反向处理,实现零相位滤波。主要的效果是获得零相位失真,也就是说,你用一个等效的滤波器过滤数据,它有一个恒定的0采样延迟。其他的效果是你得到一个滤波器传递函数它等于原滤波器传递函数的平方大小,一个滤波器的阶数是原滤波器阶数的两倍。

考虑上一节中定义的ECG信号。使用且不延迟补偿过滤此信号。设计具有75Hz的截止频率的7阶低通IIR椭圆滤波器。

fnorm = 75 /(fs / 2);%归一化频率df = designfilt (“低通”...“PassbandFrequency”,Fnorm,...“FilterOrder”7,...'passbandropple',1,...“StopbandAttenuation”、60);

绘制过滤器的组延迟。组延迟随频率而变化,表示过滤器延迟是频率依赖性的。

grpdelay (df, 2048,“一半”Fs)

图形过滤器可视化工具-组延迟包含uitoolbar、uimenu类型的轴和其他对象。标题组延迟的轴包含类型为line的对象。

过滤数据,并查看每个滤波器实现对时间信号的影响。零相位滤波有效地消除了滤波器延迟。

日元=过滤器(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])网格

图中包含一个轴。带标题过滤波形的轴包含3个line类型的对象。这些对象表示原始信号、非线性相位IIR输出、零相位IIR输出。

如果您的应用程序允许非因果的前向/后向过滤操作,并且允许将过滤器响应更改为原始响应的平方,那么零相位过滤是一个很好的工具。

引入恒定延迟的滤波器是线性相位滤波器。引入频率相关延迟的滤波器是非线性相位滤波器。

去除信号中不需要的光谱内容

滤波器通常用于去除信号中不需要的光谱内容。您可以从各种过滤器中进行选择。当你想要去除高频内容时,你可以选择低通滤波器,或者当你想要去除低频内容时,你可以选择高通滤波器。您也可以选择带通滤波器,以去除低频和高频内容,而留下一个中间频带的频率完好无损。当你想要去除给定频带上的频率时,你可以选择带阻滤波器。

考虑一个带有电力线嗡嗡声和白噪声的音频信号。电力线的嗡嗡声是由60hz的音调引起的。白噪声是存在于所有音频带宽的信号。

加载音频信号。指定采样率为44.1 kHz。

Fs = 44100;y = audioread (“noisymusic.wav”);

绘制信号的功率谱。红色三角形标记显示强烈的60hz音调干扰音频信号。

[P F] = pwelch (y)的(8192 1),8192/2,8192 Fs,“权力”);helperFilterIntroductionPlot1 (F P [60 60], [-9.365 - -9.365],...{'原始信号功率谱'“60 Hz的语气”})

图中包含一个坐标轴。轴线包含2个线型对象。这些对象代表原始信号功率谱,60hz音调。

首先可以使用低通滤波器去除尽可能多的白噪声频谱内容。滤波器的通频带应该设置为一个值,在噪声降低和由于高频内容丢失而导致的音频退化之间提供一个良好的平衡。在去除60hz的嗡嗡声之前应用低通滤波器是非常方便的,因为你将能够向下采样带限信号。较低的速率信号将允许您设计一个更尖锐和更窄的60赫兹带阻滤波器较小的滤波器订单。

设计一个通频带频率为1 kHz、阻频带频率为1.4 kHz的低通滤波器。选择最小订单设计。

《外交政策》= 1 e3;%通频带频率,单位为Hz置= 1.4 e3;%阻带频率(单位:Hz)美联社= 1;通频带纹波百分比,单位为dBAst = 95;阻带衰减百分比,单位为dBdf = designfilt (“lowpassfir”“PassbandFrequency”,fp,...'stopband职业',fst,'passbandropple'据美联社,,...“StopbandAttenuation”Ast,“SampleRate”,fs);

分析过滤器响应。

fvtool (df,“Fs”Fs,“FrequencyScale”'日志'...“FrequencyRange”“指定freq.向量”“FrequencyVector”F)

图筛选筛选器可视化工具 - 幅度响应(DB)包含UIMEnu类型uitoolbar类型的轴和其他对象。具有标题幅度响应(DB)的轴包含2个类型的线。

过滤数据并补偿延迟。

D =意味着(grpdelay (df));%过滤器延迟ylp =滤波器(df,[y; zeros(d,1)]);YLP = YLP(D + 1:结束);

看看低通滤波信号的频谱。去掉1400 Hz以上的频率含量。

[PLP,FLP] = PWELCH(YLP,ylp(8192,1),8192 / 2,8192,FS,“权力”);HelperfilterIntroductionPlot1(F,P,FLP,PLP,...{原始信号的“低通滤过的信号”})

图中包含一个坐标轴。轴线包含2个线型对象。这些对象代表原始信号,低通滤波信号。

从上面的功率谱图,您可以看到低通滤波信号的最大不可忽略的频率内容为1400 Hz。通过采样定理,采样率 2 × 1400 = 2800 Hz将足以代表信号的正确,然而,您使用的采样率为44100hz,这是一种浪费,因为您将需要处理更多的样本比那些必要的。您可以对信号进行下采样,以减少采样率,并通过减少需要处理的采样数量来减少计算负载。更低的采样率也将允许您设计一个更尖锐和更窄的带阻滤波器,需要去除60 Hz的噪声,与更小的滤波器顺序。

对低通滤波信号进行10倍的下采样,以获得Fs/10 = 4.41 kHz的采样率。绘制下采样前后的信号频谱。

Fs = f / 10;码= downsample (ylp 10);(Pds, Fds) = pwelch(码,(8192 1),8192/2,8192 Fs,“权力”);helperFilterIntroductionPlot1 (F P Fds, Pds,...{“信号以44100hz采样”'下采样信号,Fs = 4410hz '})

图中包含一个坐标轴。轴线包含2个线型对象。这些对象表示在44100hz采样的信号,下采样的信号,Fs = 4410hz。

现在,使用IIR带阻滤波器去除60 Hz的音调。让阻带的宽度为4 Hz,以60 Hz为中心。选择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”“指定freq.向量”“FrequencyVector”Fds (Fds > F (2)))

图筛选筛选器可视化工具 - 幅度响应(DB)包含UIMEnu类型uitoolbar类型的轴和其他对象。具有标题幅度响应(DB)的轴包含2个类型的线。

执行零相位滤波以避免相位失真。使用filtfilt.函数来处理数据。

yb = filtfilt (df、码);

最后,使信号使其恢复到44.1 kHz的原始音频采样率,这与音频声卡兼容。

yf=interp(ybs,10);Fs=Fs*10;

最后看一下原始和处理过的信号的频谱。高频底噪声和60赫兹的音调已被滤波器衰减。

[pfinal,ffinal] = pwelch(yf,yf,8192,1),8192 / 2,8192,fs,“权力”);helperFilterIntroductionPlot1 (F P Ffinal Pfinal,...{原始信号的“最终滤波信号”})

图中包含一个坐标轴。轴线包含2个线型对象。这些对象代表原始信号,最终滤波信号。

如果你的电脑有声卡,你可以听信号处理前后。如上所述,最终结果是有效地衰减了音频文件中的60hz嗡嗡声和高频噪声。

%要播放原始信号,取消注释下面两行% hplayer = audioplayer(y, Fs);%播放(HPlayer)%播放降噪信号,取消注释接下来的两行%hplayer = audioplayer(yf,fs);% (hplayer);

区分一个信号

MATLAB功能将信号与缺点区分开来表示您可以增加输出处的噪声水平。更好的选择是使用差分器过滤器,该滤光器充当感兴趣的乐队中的差分器,并且作为所有其他频率的衰减器,有效地去除高频噪声。

以某建筑物为例,分析了地震作用下建筑物楼板的位移速度。在地震条件下,对三层试验结构的一层进行了位移或漂移测量,并保存在地震裂缝中。垫文件。数据向量的长度为10e3,采样率为1khz,测量单位为cm。

对位移数据进行微分,得到地震时建筑物楼层的速度和加速度估计值。比较使用diff和FIR微分器滤波器的结果。

加载quakedrift.matFs = 1000;%采样率dt = 1 / fs;%的时间差异t =(0:长度(漂移)1)* dt;%的时间向量

设计一个50阶微分器滤波器,通带频率为100 Hz,这是发现大部分信号能量的带宽。将滤波器的阻带频率设置为120 Hz。

df = designfilt (“differentiatorfir”“FilterOrder”, 50岁,...“PassbandFrequency”,100,'stopband职业',120,...“SampleRate”,fs);

函数可以看作是一阶FIR滤波器的响应 H Z. = 1 - Z. - 1 .使用FVTool比较50阶微分器FIR滤波器的幅值响应与该滤波器的响应函数。显然,这两个响应在通带区域(从0到100hz)是等效的。然而,在阻带区域,50阶滤波器衰减分量,而差分响应放大分量。这有效地增加了高频噪声的水平。

HFVT = fvtool(df,[1 -1],1, 0);“MagnitudeDisplay”“零”“Fs”,fs);传奇(HFVT,“50阶FIR微分器”“diff函数的响应”);

图筛选筛选器可视化工具 - 零阶段响应包含UIMenu型UIToolbar类型的轴和其他对象。具有标题零阶段响应的轴包含2个类型的2个对象。这些对象代表50阶FIR差异化器,差异响应。

使用函数。添加零以补偿由于差异操作引起的缺失样本。

v1 = diff / dt(漂移);a1 = diff / dt (v1);v1 = [0;v1);a1 = [0;0;a1];

使用50阶FIR滤波器进行微分和延时补偿。

D =意味着(grpdelay (df));%过滤器延迟v2 =过滤器(df,漂移;0 (D, 1)]);v2 = v2 (D + 1:结束);a2 =过滤器(df, v2;0 (D, 1)]);a2 = a2 (D + 1:结束);v2 = v2 / dt;a2 = a2 / dt ^ 2;

绘制一些楼层位移的数据点。用diff和50阶FIR滤波器计算速度和加速度的一些数据点。注意噪声是如何在速度估计中被略微放大而在加速度估计中被大大放大的

helperFilterIntroductionPlot2 (t,漂移,v1、v2, a1, a2)

图中包含一个坐标轴。轴包含一个线型对象。

图中包含一个坐标轴。轴线包含2个线型对象。这些对象表示用差分估计速度,用FIR滤波器估计速度。

图中包含一个坐标轴。轴线包含2个线型对象。这些对象表示估计加速度使用diff,估计加速度使用FIR滤波器。

整合一个信号

漏式积分器滤波器是一种具有传递函数的全极滤波器 H Z. = 1 / [ 1 - C Z. - 1 ] 在哪里 C 是一个必须小于1的常数,以确保过滤器的稳定性。这并不奇怪 C 趋近于1,漏积器趋近于转换功能。将泄漏积分器应用于前一节中获得的加速度和速度估计,以分别恢复速度和漂移。使用获得的估计因为它们噪音更大。

使用泄漏的集成器 一种 = 0. 9. 9. 9. .绘制漏出的积分器滤波器的幅度响应。该滤波器作为低通滤波器有效地消除高频噪声。

fvtool (1 [1 -.999]“Fs”Fs)

图筛选筛选器可视化工具 - 幅度响应(DB)包含UIMEnu类型uitoolbar类型的轴和其他对象。具有标题幅度响应(DB)的轴包含类型线的对象。

用泄漏的积分器过滤速度和加速度。乘以时间微分。

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,Drive,D_Leakyint,V_Original,V_Leakyint)

图中包含2个轴。axis 1包含2个类型为line的对象。这些对象代表Leaky积分估计,原始位移。axis 2包含2个类型为line的对象。这些对象表示漏积分估计,原始速度。

你也可以用cumsumCumtrapz.功能。结果将与泄漏积分器得到的结果相似。

结论

在这个例子中,您学习了线性和非线性相位滤波器,并学习了如何补偿由每种滤波器类型引入的相位延迟。您还学习了如何应用滤波器从信号中去除不需要的频率成分,以及如何在用低通滤波器限制其带宽后对信号进行下采样。最后,您学习了如何使用数字滤波器设计来区分和集成信号。在整个示例中,您还学习了如何使用分析工具查看过滤器的响应和组延迟。

有关过滤器应用程序的更多信息,请参阅信号处理工具箱™文档。有关如何设计数字滤波器的更多信息,请参阅数字滤波器设计的实用介绍的例子。

工具书类

Proakis,J.G和D. G.Manolakis。数字信号处理:原理,算法和应用.Englewood Cliffs, NJ: Prentice-Hall, 1996。

Orfanidis, s . J。信号处理概论.Englewood Cliffs, NJ: Prentice-Hall, 1996。

附录

在此示例中使用以下辅助功能:

也可以看看

||||||||