基于时变滤波器的经验模态分解（TVF-EMD）附matlab代码

% TVF-EMD Algorithm
function [IMFs, Residue] = TVF_EMD(signal, numIMFs, filterOrder)
% 输入参数：
% signal: 输入信号
% numIMFs: 需要提取的 IMF 数量
% filterOrder: 时变滤波器的阶数

% 设置停止条件
epsilon = 0.001;  % 停止条件的阈值
maxIterations = 100;  % 最大迭代次数

% 初始化
IMF = signal;  % 当前的 IMF
IMFs = [];  % 存储提取的 IMF
Residue = signal;  % 存储剩余项

% 迭代提取 IMFs
for k = 1:numIMFs
    iterations = 0;  % 当前迭代次数
    while iterations < maxIterations
        % 计算当前 IMF 的 Hilbert 谱
        hilbSpectrum = abs(hilbert(IMF));

        % 计算滤波器的时变系数
        filterCoeffs = computeFilterCoeffs(hilbSpectrum, filterOrder);

        % 应用时变滤波器
        filteredSignal = applyTimeVaryingFilter(IMF, filterCoeffs);

        % 计算当前 IMF 与滤波后信号的差异
        residue = IMF - filteredSignal;

        % 判断是否满足停止条件
        if norm(residue) / norm(IMF) < epsilon
            break;
        end

        % 更新当前 IMF
        IMF = residue;

        iterations = iterations + 1;
    end

    % 将当前 IMF 存入 IMFs
    IMFs = [IMFs, IMF];

    % 更新剩余项
    Residue = Residue - IMF;
    
    % 检查剩余项是否为全零信号
    if norm(Residue) == 0
        break;
    end

    % 更新当前 IMF 为剩余项
    IMF = Residue;
end

end

% 计算滤波器的时变系数
function filterCoeffs = computeFilterCoeffs(hilbSpectrum, filterOrder)
% 计算时变系数
filterCoeffs = abs(diff(log(hilbSpectrum), filterOrder));
filterCoeffs = filterCoeffs / sum(filterCoeffs); % 归一化
end

% 应用时变滤波器
function filteredSignal = applyTimeVaryingFilter(signal, filterCoeffs)
signalLength = length(signal);
filteredSignal = zeros(1, signalLength);

% 对信号的每个样本应用滤波器
for n = 1:signalLength
    filterLength = min(n, length(filterCoeffs));
    filteredSignal(n) = sum(signal(n-filterLength+1:n) .* filterCoeffs(end-filterLength+1:end));
end

end

% 示例输入信号
t = 0:0.01:10;
f1 = 1;
f2 = 10;
f3 = 20;
signal = sin(2pif1t) + sin(2pif2t) + sin(2pif3*t);

% 调用 TVF-EMD 函数进行信号分解
numIMFs = 3;
filterOrder = 10;
[IMFs, Residue] = TVF_EMD(signal, numIMFs, filterOrder);

% 绘制分解后的 IMFs 和剩余项
figure;
subplot(numIMFs+1, 1, 1);
plot(t, signal);
title(‘Original Signal’);
xlabel(‘Time’);
ylabel(‘Amplitude’);
for k = 1:numIMFs
subplot(numIMFs+1, 1, k+1);
plot(t, IMFs(:, k));
title(['IMF ', num2str(k)]);
xlabel(‘Time’);
ylabel(‘Amplitude’);
end