主要内容

数字滤波实用介绍

此示例演示如何设计、分析数据并对数据应用数字滤波器。它将帮助您回答以下问题:如何补偿滤波器引入的延迟?、如何避免信号失真?、如何从信号中删除不需要的内容?、如何区分信号?以及如何集成信号?

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

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

滤波引起的延迟补偿

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

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

恒定滤波器延迟的补偿

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

考虑一个有噪声的心电信号,你想过滤掉75赫兹以上的高频噪声。你想应用一个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=设计过滤器(“lowpassfir”,“FilterOrder”,70,“截止频率”,Fnorm);

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

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

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

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

在滤波之前,在输入数据向量x的末尾附加D个零。这确保所有有用的样本都从滤波器中清除,并且输入信号和延迟补偿的输出信号具有相同的长度。通过将输出信号移位D个样本,对数据进行滤波并补偿延迟。这是最后一个有效的步骤移除过滤器瞬态。

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

图中包含一个轴对象。带有标题过滤波形的轴对象包含2个line类型的对象。这些对象表示原始噪声信号、滤波信号。

补偿频率相关延迟

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

考虑在前一节中定义的ECG信号,在没有延迟补偿的情况下对该信号进行滤波。设计一个截止频率为75赫兹的第七阶低通IIR椭圆滤波器。

Fnorm=75/(Fs/2);%归一化频率df=设计过滤器(“低通”,...“PassbandFrequency”,Fnorm,...“FilterOrder”7....“通带波纹”,1,...“StopbandAttenuation”,60);

绘制滤波器的群延迟。群延迟随频率变化,表明滤波器延迟与频率相关。

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

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

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

日元=过滤器(df, x);%非线性相位滤波器-无延迟补偿y2=滤波器(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”);

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

[P,F]=pwelch(y,one(8192,1),8192/28192,Fs,“权力”);helperFilterIntroductionPlot1(F,P[60],-9.365-9.365],...{“原始信号功率谱”,“60赫兹音调”})

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表原始信号功率谱,60hz音调。

您可以首先使用低通滤波器删除尽可能多的白噪声频谱内容。滤波器的通带应设置为一个值,该值可在降噪和由于高频内容丢失而导致的音频降级之间提供良好的折衷。在删除60 Hz嗡嗡声之前应用低通滤波器非常方便,因为您可以I’我们能够对带限信号进行降采样。较低的速率信号将允许您设计一个更清晰、更窄的60 Hz带阻滤波器,滤波器阶数更小。

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

Fp=1e3;%通频带频率,单位为HzFst=1.4e3;%阻带频率美联社= 1;%通带纹波(dB)Ast = 95;阻带衰减百分比,单位为dBdf=设计过滤器(“lowpassfir”,“PassbandFrequency”,Fp,...“阻带频率”,Fst,“通带波纹”据美联社,,...“StopbandAttenuation”,Ast,“采样器”Fs);

分析过滤器响应。

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

图形过滤器可视化工具-幅度响应(dB)包含一个轴对象和其他类型的uitoolbar, uimenu对象。标题为“大小响应(dB)”的轴对象包含2个类型为line的对象。

过滤数据并补偿延迟。

D=平均值(grpdelay(df));%过滤器延迟ylp=滤波器(df[y;零(D,1)];ylp=ylp(D+1:结束);

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

[Plp,Flp]=pwelch(ylp,one(8192,1),8192/28192,Fs,“权力”);helperFilterIntroductionPlot1(F、P、Flp、Plp、,...{“原始信号”,“低通滤波信号”})

图中包含一个Axis对象。Axis对象包含2个line类型的对象。这些对象表示原始信号,即低通滤波信号。

从上面的功率谱图可以看出,低通滤波信号的最大不可忽略频率含量为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,...{'在44100 Hz下采样的信号','下采样信号,Fs = 4410hz '})

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象表示在44100hz采样的信号,下采样的信号,Fs = 4410hz。

现在,使用IIR带阻滤波器去除60 Hz的音调。让阻带的宽度为4 Hz,以60 Hz为中心。选择IIR滤波器以获得尖锐的频率陷波、较小的通带纹波和相对较低的阶数。

df=设计过滤器(“bandstopiir”,“PassbandFrequency1”现年55岁的...“阻带频率1”,58,“StopbandFrequency2”, 62,...“PassbandFrequency2”, 65,“PassbandRipple1”,1,...“StopbandAttenuation”,60,“PassbandRipple2”,1,...“采样器”Fs,“设计方法”,“ellip”);

分析震级响应。

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

图形过滤器可视化工具-幅度响应(dB)包含一个轴对象和其他类型的uitoolbar, uimenu对象。标题为“大小响应(dB)”的轴对象包含2个类型为line的对象。

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

ybs=滤波器(df,yds);

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

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

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

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

图中包含一个轴对象。axes对象包含2个line类型的对象。这些对象表示原始信号、最终滤波信号。

如果您的计算机有声卡,您可以在处理前后收听信号。如上所述,最终结果是您有效地衰减了音频文件中的60 Hz嗡嗡声和高频噪声。

%要播放原始信号,取消注释下面两行% hplayer = audioplayer(y, Fs);%播放(hplayer)%要播放降噪信号,请取消对下两行的注释%hplayer=音频播放器(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=设计过滤器(“differentiatorfir”,“FilterOrder”, 50岁,...“PassbandFrequency”,100,“阻带频率”,120,...“采样器”Fs);

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

hfvt=fvtool(df,[1-1],1,“MagnitudeDisplay”,“零”,“Fs”,Fs);图例(hfvt,“50阶FIR微分器”,“diff函数的响应”);

图过滤器可视化工具-零相位响应包含一个轴对象和其他类型的uitoolbar, uimenu对象。标题为“零相位响应”的轴对象包含两个类型为line的对象。这些对象代表50阶FIR微分器,差分函数的响应。

使用差异函数。添加零以补偿由于diff操作而丢失的样本。

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

使用第50阶FIR滤波器进行区分并补偿延迟。

D=平均值(grpdelay(df));%过滤器延迟v2=滤波器(df[drift;zeros(D,1)]);v2=v2(D+1:end);a2=滤波器(df[v2;zeros(D,1)];a2=a2(D+1:end);v2=v2/dt;a2=a2/dt^2;

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

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

图中包含一个axes对象。axes对象包含一个line类型的对象。

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象表示用差分估计速度,用FIR滤波器估计速度。

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

整合一个信号

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

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

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

图形过滤器可视化工具-幅度响应(dB)包含一个轴对象和其他类型的uitoolbar, uimenu对象。标题为“大小响应(dB)”的轴对象包含一个类型为line的对象。

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

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_泄漏、v_原始、v_泄漏)

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

你也可以用cumsum累计梯形积分功能。结果将与泄漏积分器得到的结果相似。

结论

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

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

工具书类

Proakis、J.G.和D.G.Manolakis。数字信号处理:原理、算法和应用. 恩格尔伍德悬崖,新泽西州:普伦蒂斯大厅,1996年。

Orfanidis, s . J。信号处理导论. 恩格尔伍德悬崖,新泽西州:普伦蒂斯大厅,1996年。

附录

本例中使用了以下辅助函数:

另见

||||||||