主要内容

数字滤波实用介绍“,

这个例子展示了如何设计、分析和应用数字滤波器到您的数据。它将帮助您回答以下问题:我如何补偿由过滤器带来的延迟?如何避免信号失真?,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赫兹以上的高频噪声。你想要应用一个FIR低通滤波器,并补偿滤波器延迟,以便噪声和滤波信号正确对齐,并可以绘制在彼此的顶部进行比较。

Fs = 500;%采样率,单位为HzN = 500;%信号样本数rng默认的;x = ecg(N)'+0.25*randn(N,1);%噪声波形t = (0:N-1)/Fs;%时间矢量

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

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

绘制滤波器的组延迟,以验证它在所有频率上是恒定的,表明滤波器是线性相位。使用组时延来度量过滤器的时延。

grpdelay (df、2048 Fs)% Plot group延迟

图过滤器可视化工具-组延迟包含一个axis对象和其他类型为uitoolbar、uimenu的对象。标题为Group delay的axes对象包含一个类型为line的对象。

D = mean(grpdelay(df))样本中的滤波延迟
D = 35

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

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

图中包含一个axes对象。标题为Filtered Waveforms的axis对象包含两个类型为line的对象。这些对象分别代表原始噪声信号、滤波信号。

补偿与频率相关的延迟

频率依赖性延迟导致信号相位失真。对这种延迟的补偿并不像对恒定延迟的补偿那么简单。方法实现零相位滤波,如果应用程序允许离线处理,则可以消除与频率相关的延迟filtfilt函数。filtfilt通过正向和反向处理输入数据来执行零相位滤波。其主要效果是获得零相位失真,即使用具有恒定0个样本延迟的等效滤波器对数据进行滤波。其他的效果是,你得到一个滤波器传递函数等于原始滤波器传递函数的大小的平方,一个滤波器阶数是原始滤波器阶数的两倍。

考虑上一节中定义的心电信号。有或无延迟补偿对该信号进行滤波。设计一个截止频率为75 Hz的七阶低通IIR椭圆滤波器。

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

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

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

图过滤器可视化工具-组延迟包含一个axis对象和其他类型为uitoolbar、uimenu的对象。标题为Group delay的axes对象包含一个类型为line的对象。

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

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])网格

图中包含一个axes对象。标题为Filtered Waveforms的axis对象包含3个类型为line的对象。这些对象代表原始信号,非线性相位IIR输出,零相位IIR输出。

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

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

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

滤波器通常用于去除信号中不需要的光谱内容。您可以从各种各样的过滤器中进行选择。如果要删除高频内容,则选择低通滤波器;如果要删除低频内容,则选择高通滤波器。您还可以选择带通滤波器来去除低频和高频内容,同时保留中间频带的频率不变。当您想要去除给定频带上的频率时,可以选择带阻滤波器。

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

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

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

画出信号的功率谱。红色三角形标记表示干扰音频信号的60hz强音调。

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

图中包含一个axes对象。坐标轴对象包含两个line类型的对象。这些物体代表原始信号功率谱,60赫兹的音调。

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

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

Fp = 1e3;%通带频率,单位为HzFst = 1.4e3;%阻带频率,单位为HzAp = 1;%通带波纹在dBAst = 95;%阻带衰减单位为dBDf = designfilt(“lowpassfir”“PassbandFrequency”《外交政策》,...“StopbandFrequency”置,“PassbandRipple”据美联社,,...“StopbandAttenuation”Ast,“SampleRate”Fs);

分析过滤器响应。

fvtool (df,“Fs”Fs,“FrequencyScale”“日志”...“FrequencyRange”'指定频率。vector'“FrequencyVector”F)

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

筛选数据并补偿延迟。

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,Fs,“权力”);Flp helperFilterIntroductionPlot1 (F P、Plp...原始信号的“低通滤波信号”})

图中包含一个axes对象。坐标轴对象包含两个line类型的对象。这些对象分别代表原始信号、低通滤波信号。

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

对低通滤波后的信号下采样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”})

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

现在使用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”'指定频率。vector'“FrequencyVector”Fds (Fds > F (2)))

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

进行零相位滤波,避免相位失真。使用filtfilt函数处理数据。

Ybs = filtfilt(df,yds);

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

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

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

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

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

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

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

信号的微分

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

作为一个例子,分析了建筑物在地震时的位移速度。位移或漂移测量记录在地震条件下三层试验结构的一层,并保存在地震裂谷中。垫文件。数据向量的长度为10e3,采样率为1khz,测量单位为cm。

微分位移数据,以获得地震期间建筑物地板的速度和加速度的估计。比较使用差分和FIR微分器滤波器的结果。

负载quakedrift.matFs = 1000;抽样率%dt = 1/Fs;%时间差异T =(0:长度(漂移)-1)*dt;%时间矢量

设计一个通频带频率为100hz的50阶微分器滤波器,这是大部分信号能量被发现的带宽。将滤波器的阻带频率设置为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,“MagnitudeDisplay”“零”“Fs”Fs);传奇(hfvt“50阶FIR微分器”“diff函数的响应”);

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

使用diff函数。添加0来补偿由于差值操作而丢失的样本。

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:end);V2 = V2 /dt;A2 = A2 /dt^2;

画出地板位移的几个数据点。还绘制了一些用diff和用50阶FIR滤波器计算的速度和加速度的数据点。注意噪声是如何在速度估计中被轻微放大,而在加速度估计中被大大放大的diff

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

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

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

图中包含一个axes对象。坐标轴对象包含两个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)的axis对象包含一个类型为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 = * dt;V_leakyint = V_leakyint * dt;

绘制位移和速度估计图,并与原始信号进行比较。

helperFilterIntroductionPlot3 (t,漂移,d_leakyint、v_original v_leakyint)

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

函数也可以对信号进行积分cumsum而且cumtrapz功能。结果将与使用泄漏积分器得到的结果相似。

结论

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

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

参考文献

普罗基斯,J. G.和D. G.马诺拉基斯。数字信号处理:原理,算法和应用.恩格尔伍德悬崖,新泽西州:Prentice-Hall, 1996。

奥法尼迪斯,s.j。信号处理概论.恩格尔伍德悬崖,新泽西州:Prentice-Hall, 1996。

附录

本例中使用了以下helper函数:

另请参阅

||||||||