小波包谐波干扰去除
这个例子展示了如何使用小波包来消除信号中的谐波干扰。谐波干扰分量是干扰信号的正弦分量。这些伪成分对后续的数据处理和分析有不良影响。
通常,这些谐波干扰成分位于信号的频谱内。因此,需要有技术来消除或减轻这些杂散谐波的影响,而不会对这些谐波附近的初级信号的频率含量产生不利影响。
解决这个问题的常用方法是应用缺口滤波器或梳状滤波器。虽然陷波滤波器是有效的,但不可能使陷波滤波器无限窄,或者等效地,设计一个具有无限Q因子的滤波器。因此,信号的频谱在干扰频率附近不可避免地会有一些失真。此外,陷波滤波器通常是无限脉冲响应(IIR)滤波器,具有非线性相位响应。这在不希望出现相位失真的应用程序中可能会出现问题。
小波包变换与基线移动
在这个例子中,我们使用小波包方法[3.],它有可能是无失真的。看到小波包:细节分解对小波包变换进行了回顾。
谐波干扰小波包技术的基础包括以下步骤:
注意数据采样率和谐波干扰的基频, .
确定小波包变换的最小新采样率和最小正电平J,其中小波包变换的抽取样本对应于的采样 周期的整数倍。
以相应的速率重新采样数据。
使用迭代过程确定每个细节子带中系数的基线水平。
减去细节子带系数的基线,重建数据。
以原始速率重新采样数据。咨询[3]获取技术细节。
上述方法基于小波包变换的详细子带系数为零均值,重采样数据中存在非零基线表明存在谐波干扰的事实。
60hz干扰组件
加载并绘制与MIT-BIH心律失常数据库记录200对应的心电信号[1],[2].采样频率为360hz。
负载mit200ecgsig情节(tm)包含(“秒”) ylabel (“振幅”)轴紧标题(“有60赫兹干扰的心电信号”)
绘制数据的功率谱,以便清楚地看到60hz干扰。注意,这种干扰是不人工注入的数据。它出现在原始录音中。
Pspectrum (ecgsig,360) xlim([45 80]) title(《60赫兹干扰下心电信号功率谱》)
在120hz的数据中也有60hz的谐波可见。然而,这种干扰水平大约比60hz干扰水平低20 dB。因此,我们在这里重点讨论60hz分量。
首先,加载并分析了两种IIR陷波滤波器iirnotch
从DSP系统工具箱™。滤波器设计有两个Q因子,35和100,默认带宽为-3 dB。
负载Q100Filter3负载Q35Filter3
检查带宽水平为-3 dB的两个滤波器的幅度和相位响应。
helperFreqPhasePlot (Q100Filter3.num Q100Filter3.den,…Q35Filter3.num Q35Filter3.den);
采用q因子为100的陷波滤波器对心电数据进行零相位滤波。绘制滤波后数据的功率谱。
y = filtfilt(Q100Filter3.num,Q100Filter3.den,ecgsig);图pspectrum([ecgsig,y],360) xlim([55 65]) title(心电信号与缺口滤波信号的功率谱)传说(心电信号的,“切口过滤”)
陷波滤波器在去除60hz干扰方面做得很好,但在干扰频率附近,信号的频谱显然有一些损坏。
如果缩小到从0到Nyquist的整个功率谱,可以清楚地看到信号频谱在60 Hz附近的失真,而滤波信号的整体频谱与原始信号的频谱匹配得很好。
xlim (180 [0])
辅助函数helperWPHarmonicFilter
实现了中描述的基线移位小波包方法[3].本文采用8个消失矩(16个系数)的Daubechies最小不对称小波。
sigWP = helpwpharmonicfilter (ecgsig,60,“sym8”, 360,“反射”);pspectrum([ecgsig,sigWP],360) xlim([55 65]) title({心电信号功率谱及;“基线移位小波包信号”})(传说心电信号的,“小波包消除”,“位置”,“西南”)
注意,小波包方法去除了谐波干扰没有使信号的频谱在60赫兹左右失真。
缩小以查看小波包方法的功率谱质量与原始信号从0到奈奎斯特频率的比较。
xlim (180 [0])
小波包方法的整体谱与信号谱匹配得很好。奈奎斯特附近有一些小变形。在这一点上,还不清楚这是由于基线移动算法还是重新采样过程。这值得进一步调查。然而,在干扰频率附近没有明显的陷波,因为有陷波滤波器。
小波包方法的优点之一是它在理论上可以同时处理多个谐波。为了说明这一点,加载信号ecgHarmonics
.这是mit200
在60和120 Hz注入强谐波干扰成分的信号。
负载ecgHarmonics图pspectrum(ecgHarmonics,360)标题(“人工谐波心电信号功率谱”)
用小波包方法去除谐波干扰。再次,我们使用Daubechies最小不对称相位小波与8消失矩。
sigWP = helperWPHarmonicFilter(ecgHarmonics,60,“sym8”, 360,“反射”);
绘制原始信号和基线移位小波包重构图。
subplot(2,1,1) plot(tm,ecgHarmonics) ylabel(“振幅”)轴紧标题(“有干扰分量的信号”) subplot(2,1,2) plot(tm,sigWP) xlabel(“秒”) ylabel (“振幅”)标题(“基线移位小波包信号”)轴紧
检查两个信号的功率谱。
图pspectrum([ecgHarmonics,sigWP],360) xlim([30 150]) title({心电信号功率谱及;“小波包基线移位信号”})(传说“有干扰的心电信号”,“基线移位小波包信号”,“位置”,“东北”)
该方法去除了两种谐波,且在干扰谐波附近没有频谱畸变,效果很好。
将其与IIR陷波滤波进行比较。在这里,我们需要对数据进行两次滤波,一次是工作在60 Hz的陷波滤波器,一次是工作在120 Hz的陷波滤波器。或者,您也可以创建一个过滤器级联。这里我们想用filtfilt
对于零相位滤波,我们选择实现两次滤波。
加载120hz陷波滤波器设计使用iirnotch
Q系数为100,默认带宽级别为-3 dB。首先用60hz陷波滤波器对信号进行滤波,然后用120hz陷波滤波器对该操作的输出进行滤波。
负载Q100Filter3_120Hzfilt60Hz = filtfilt(Q100Filter3.num,Q100Filter3.den,ecgHarmonics);filtCascade = filtfilt(Q100Filter3_120Hz.num,Q100Filter3_120Hz.den,filt60Hz);
绘制原始信号和缺口滤波信号。
subplot(2,1,1) plot(tm,ecgHarmonics) ylabel(“振幅”)轴紧标题(“有干扰分量的信号”subplot(2,1,2) plot(tm,filtCascade) xlabel(“秒”) ylabel (“振幅”)标题(“Notch-filtered信号”)轴紧
检查两个信号的功率谱。
图pspectrum([ecgHarmonics,filtCascade],360) xlim([30 150]) title(心电信号与缺口滤波信号的功率谱)传说(心电信号的,“Notch-Filtered信号”)
在这里,我们再次看到在陷波滤波器位置的频谱失真。小波包系数基线移位对去除高振幅正弦干扰分量更为有效。
总结
中描述的基线偏移小波包技术来消除谐波干扰[3].在本文所考虑的心电数据中,小波包技术在真实情况和模拟情况下都能有效去除谐波干扰。小波包技术在干扰频率附近的信号频谱中引入的谐波失真比陷波滤波器少。小波包技术在奈奎斯特频率附近引入了一些失真,值得进一步研究。
参考文献
Goldberger, Ary L., Luis A. N. Amaral, Leon Glass, Jeffrey M. Hausdorff, Plamen Ch. Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody,彭仲康和H. Eugene Stanley。“PhysioBank, PhysioToolkit,和PhysioNet:复杂生理信号新研究资源的组成部分。”循环101年,没有。23(2000年6月13日)https://doi.org/10.1161/01.CIR.101.23.e215.
[2]穆迪,g.b.和R.G.马克。MIT-BIH心律失常数据库的影响IEEE医学与生物工程杂志20日,没有。3(2001年6月):45-50。https://doi.org/10.1109/51.932724.
[3]许丽君。小波包分解系数基线移位消除谐波干扰。IEEE信号处理汇刊53岁的没有。1(2005年1月):222-30。https://doi.org/10.1109/TSP.2004.838954.
函数helpfreqphaseplot (num1,den1,num2,den2) figure [H100,W] = freqz(num1,den1,[],360);Phi100 = phasez(num1,den1);H35 = freqz(num2,den2,[],360);Phi35 = phasez(num2,den2);subplot(2,2,1) plot(W,10*log10(abs(H100)))“幅度- q因子100”) ylabel (“数据库”网格)在轴紧subplot(2,2,2) plot(W,phi100*180/pi)“相位- q因子100”) ylabel (“度”网格)在轴紧subplot(2,2,3) plot(W,10*log10(abs(H35)))“震级- q因子35”)包含(“赫兹”) ylabel (“数据库”网格)在轴紧subplot(2,2,4) plot(W,phi35*180/pi)“相位- q因子35”)包含(“赫兹”) ylabel (“度”网格)在轴紧结束