基于小波包分析的一维信号谐波干扰去除方法(MATLAB R2021B)

小波包是由Coifman、Meyer及Wickhauser提出的。他们在研究正交小波基的基础上创立了正交小波包的概念，后来又发展到半正交小波包和广义小波包。

由于正交小波变换采用Mallat算法实现，只对信号的低频部分做进一步分解，而对高频部分不再继续分解，所以小波变换仅能很好的表征以低频信息为主要成分的信号。但是在实际应用中非平稳机械振动信号、遥感图像、地震信号和生物医学信号等包含大量细节信息，正交小波变换不能很好地分解和表示细小边缘或纹理的信号。小波包变换可以对高频部分提供更精细的分解，而且这种分解既无冗余，也无疏漏，所以对包含大量中、高频信息的信号能够进行更好的时频局部化分析。小波包变换在信号去噪、滤波、压缩、非平稳机械振动信号的分析与故障诊断、非平稳信号的特征提取及多载波调制技术等方面具有重要应用。

信号小波包分析的基本实现步骤如下:

(1)选择适当的小波滤波器，对给定的采样信号进行小波包变换，获得树形结构的小波包系数。常用的正交滤波器包括Daubechies小波滤波器、Symlets小波滤波器等，小波滤波器的性能通常与所分析的信号类型相关，可通过实验比较选取合适的滤波器。

(2)选择信息代价函数，并利用最佳小波包基选取算法选取最佳基。我们曾给出几个常用的信息代价函数的选取方法。需要说明的是，信息代价函数的选取与小波包的应用是相关的，而且人们也在不断提出新的选取方法。通常需要通过实验比较选取合适的代价函数。当代价函数选定后，可采用自底向顶的方法搜索最佳小波包基，该基能够提供所分析信号的高效表示。

(3)对最佳正交小波包基对应的小波包系数进行处理。处理方法与具体应用有关，如在去噪时，可对系数进行阈值化处理。

鉴于此，采用小波包分析对一维时间序列信号进行谐波干扰去除，运行环境为MATLAB R2021B。

function [P,F,PWR] = computePersistence(opts, plotFlag)
% Compute persistence spectrum

% Cache probability min threshold and set it to zero in opts so that
% spectrogram is computed with no threshold.
minThreshold = opts.MinThreshold;
opts.MinThreshold = 0;

[Pspectrogram,F,~,FRES,TRES,zMin,zMax] = computeSpectrogram(opts, false);
Npoints = length(F);
numSpectralWindows = size(Pspectrogram,2);

f1 = opts.FrequencyLimits(1);
f2 = opts.FrequencyLimits(2);

% Add 5% cushion so that there is some empty space above and below the
% persistence spectrum image.
pwrCushion = 0.05*abs(zMax-zMin);
zMin = zMin-pwrCushion;
zMax = zMax+pwrCushion;

if coder.target('MATLAB')
    persistenceObj = signalanalyzer.internal.PersistenceSpectrum;
else
    persistenceObj = signal.internal.codegenable.pspectrum.PersistenceSpectrum;
end

persistenceObj.setup(f1,f2,Npoints,zMin,zMax,opts.NumPowerBins);
persistenceObj.computeSpectrum(Pspectrogram(:));
P = persistenceObj.fetchPersistenceSpectrum2D();
P = 100*(P/numSpectralWindows); % convert to probability in percentage

F = persistenceObj.fetchFrequencyVector();
PWR = persistenceObj.fetchMagnitudeVector(); % this is in dB

if minThreshold > 0
    P(P<minThreshold) = 0;
end

if plotFlag
    if opts.IsSingle
        P = single(P);
        F = single(F);
        PWR = single(PWR);
    end
    signal.internal.pspectrum.displayPersistence(PWR,F,P,...
        opts.IsNormalizedFreq,opts.MinThreshold,FRES,TRES,numSpectralWindows);
end

% Convert PWR to linear before returning
PWR = 10.^(PWR/10);
%完整代码可通过知乎学术咨询获得：https://www.zhihu.com/consult/people/792359672131756032?isMe=1
end

工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。