主要内容

数字滤波实用导论

这个例子展示了如何设计、分析和应用数字滤波器到您的数据。它将帮助您回答以下问题:我如何补偿滤波器引入的延迟?、如何避免信号失真?,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,并通过在输入信号上添加d0并通过D采样在时间上移动输出信号来补偿这个延迟。

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

Fs = 500;%采样频率,单位为HzN = 500;%信号采样个数rng默认的;x =心电图(N) + 0.25 * randn (N, 1);%的波形t = (0: n - 1) / Fs;%的时间向量

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

Fnorm = 75 / (Fs / 2);%归一化频率df = designfilt (“lowpassfir”“FilterOrder”, 70,“CutoffFrequency”, Fnorm);

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

grpdelay (df、2048 Fs)图群延迟

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

D =意味着(grpdelay (df))采样滤波器延迟%
D = 35

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

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

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

补偿频率相关延迟

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

考虑上一节定义的心电信号。滤波器这个信号与没有延迟补偿。设计一个截止频率为75hz的七阶低通IIR椭圆滤波器。

Fnorm = 75 / (Fs / 2);%归一化频率df = designfilt (“lowpassiir”...“PassbandFrequency”Fnorm,...“FilterOrder”7...“PassbandRipple”, 1...“StopbandAttenuation”、60);

绘制滤波器的组延迟。组延迟随频率变化,说明滤波器延迟与频率有关。

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

图形过滤可视化工具-组延迟包含一个轴对象和其他类型的uitoolbar, uimenu对象。具有标题组延迟的axis对象包含一个类型为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的语气”})

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

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

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

《外交政策》= 1 e3;%通频带频率,单位为Hz置= 1.4 e3;停止频带频率,单位为Hz美联社= 1;通频带纹波百分比,单位为dBAst = 95;阻带衰减百分比,单位为dBdf = designfilt (“lowpassfir”“PassbandFrequency”《外交政策》,...“StopbandFrequency”置,“PassbandRipple”据美联社,,...“StopbandAttenuation”Ast,“SampleRate”Fs);

分析过滤器响应。

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

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

过滤数据并补偿延迟。

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

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

Flp (Plp) = pwelch (ylp (8192 1), 8192/2, 8192 Fs,“权力”);Flp helperFilterIntroductionPlot1 (F P、Plp...原始信号的“低通滤过的信号”})

图中包含一个轴对象。轴对象包含两个类型为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,...“信号以44100hz采样”'下采样信号,Fs = 4410hz '})

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

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

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

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

yb = filtfilt (df、码);

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

yf =插值函数(yb, 10);Fs = f * 10;

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

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

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表原始信号,最终滤波信号。

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

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

区分一个信号

MATLABdiff函数区分信号的缺点是,您可能会在输出时增加噪声水平。一个更好的选择是使用微分器滤波器,在感兴趣的频带中作为微分器,在所有其他频率上作为衰减器,有效地去除高频噪声。

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

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

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

设计一个50阶微分器滤波器,其通频带频率为100hz,这是信号能量的主要来源。设置滤波器的阻带频率为120hz。

df = designfilt (“differentiatorfir”“FilterOrder”, 50岁,...“PassbandFrequency”, 100,“StopbandFrequency”, 120,...“SampleRate”Fs);

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

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

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

区分使用diff函数。添加零来补偿由于diff操作而丢失的样本。

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滤波器计算速度和加速度的一些数据点。注意噪声是如何在速度估计中被略微放大而在加速度估计中被大大放大的diff

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

图中包含一个轴对象。axis对象包含一个类型为line的对象。

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

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

整合一个信号

漏式积分器滤波器是一种具有传递函数的全极滤波器 H Z 1 / 1 - c Z - 1 在哪里 c 为必须小于1的常数,以确保滤波器的稳定性。这是不足为奇的 c 趋近于1,漏积器趋近于diff传递函数。将漏积器应用于前一节得到的加速度和速度估计,分别得到速度和漂移。使用得到的估计diff因为它们噪音更大。

使用漏泄积分器 一个 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_leakyint、v_original v_leakyint)

图中包含2个轴对象。坐标轴对象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。

附录

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

另请参阅

||||||||