基于ASSA-WFEE的新型VMD优化策略及其在滚动轴承故障诊断中的应用

项目简介

❝

用于支持ASSA算法的先进性及其在轴承故障诊断应用中的合理性的代码和数据，其中ASSA存储了基于测试函数的ASSA算法与WOA、GWO、PSO、HHO、SSA、MFO的比较测试。PARPOR存储了项目中使用的轴承的振动信号数据以及用于故障特征提取的不同算法。比较代码的效果，迭代次数都是50次，使用WFEE作为适应度函数。

其他系列文章集合详见：

• 公众号：[「滚动轴承故障诊断与寿命预测」]

前言

考虑经验模态分解(EMD)与集合经验模态分解(EEMD)存在缺乏自适应性、端点效应和模态混叠等问题, Dragomiretskiy等提出了一种变尺度非平稳信号处理方法—变分模态分解(VMD)。VMD方法基于无冗余方式对信号进行正交分解, 可反映原信号在不同频带上的分布特征。但VMD方法需要人为设置模态分解数 K 与惩罚因子c,二者对信号分解精度影响显著。Huang 等基于改进尺度空间预设模态分解数 K,虽然有效提高了 VMD 的自适应性,但是惩罚因子c仍靠经验得到。唐贵基等采用遗传算法(GA)对 VMD 参数进行优化,但未考虑变异率及迭代次数的影响。张俊等利用粒子群算法(PSO)优化 VMD 参数,选取最佳参数组合,虽然避免了人为经验影响,但是仍易陷入局部最优,导致误差较大。

VMD理论分析基础

优化变分模态分解

其他创新模型结果

clear
clc
close all
addpath(genpath(pwd)) % 设置访问路径


%     % 齿轮振动信号
load 276.mat  % inner
%         load 288.mat  % ball
%                 load 272.mat  % inner 21
%        load 209.mat  %
X_Gear = X276_FE_time(20001:30000); % 10000
%         X_Gear = X209_FE_time(20001:30000); % 10000
X_Gear = X_Gear';
%
num_point=10000;
fs = 12000; % 假设采样频率为1000 Hz
f = (0:num_point-1)*(fs/num_point); % f是一个频率向量
%% figure 1 原始信号图及其频域图

figure(1);
x2 = zeros(1, num_point); % 预分配空间
t2 = zeros(1, num_point); % 预分配空间——在Matlab中，对于循环内部需要动态增加大小的变量，
% 频繁的重新分配内存会导致程序的性能下降。因此，为了避免这种情况，我们可以在循环之
% 前先预分配足够大的空间。
for i=1:num_point
    x2(i)=X_Gear(i);
    t2(i)=i/fs;
end
subplot(1,2,1)
plot(t2,x2);xlabel('Time/s','FontSize', 12, 'FontName', 'Arial','LineWidth', 0.6) % 画出输入信号
ylabel('Amplitude','FontSize', 12, 'FontName', 'Arial'); % 画出输入信号
% 设置纵坐标轴的上下限
% yLowerLimit = -0.5;  % 设置下限
% yUpperLimit = 0.5; % 设置上限
% ylim([yLowerLimit,yUpperLimit]);
% 设置横坐标轴的上下限
xLowerLimit = 0;  % 设置下限
xUpperLimit = 0.834; % 设置上限
xlim([xLowerLimit,xUpperLimit]);

% 计算信号的频域表示
N_hat = length(X_Gear);    % 信号的样本数
freq_hat = (0:N_hat-1)*(fs/N_hat); % 频率向量
fft_hat = abs(fft(X_Gear)); % 取频谱的幅值并进行归一化
% 只绘制频谱的一半，并处理幅值
half_N_Gear = floor(N_hat/2);   % 取频谱的一半
fft_signal_half_Gear = abs(fft_hat(1:half_N_Gear))/N_hat; % 取频谱的幅值并进行归一化
% 绘制频域图
subplot(1,2,2)
plot(freq_hat(1:half_N_Gear), fft_signal_half_Gear);
% title('信号的频域图');
xlabel('Frequency(Hz)','FontSize', 12, 'FontName', 'Arial');
ylabel('Amplitude','FontSize', 12, 'FontName', 'Arial');
% 设置纵坐标轴的上下限
% yLowerLimit = 0;  % 设置下限
% yUpperLimit = 1600; % 设置上限
% ylim([yLowerLimit,yUpperLimit]);
% 设置横坐标轴的上下限
% xLowerLimit = 0;  % 设置下限
% xUpperLimit = 50; % 设置上限
% xlim([xLowerLimit,xUpperLimit]);

% 设置图像的宽度和高度（以厘米为单位）
width_inches = 26;
height_inches = 4;
set(gcf, 'Units', 'centimeters', 'Position', [0, 0, width_inches, height_inches]);
%% 绘制子模态图和频谱图
tau=0;
DC=0;
init=1;
tol=1e-7;
X=X_Gear;

K_PSO=7;
alpha_PSO=381;
K_GWO=4;
alpha_GWO=3574;
K_ABC=6;
alpha_ABC=871;
K_AOA=2;
alpha_AOA=3836;
K_ASSA=2;
alpha_ASSA=685;

[u_PSO, u_hat_PSO, omega_PSO] = VMD(X, alpha_PSO, tau, K_PSO, DC, init, tol);
[u_GWO, u_hat_GWO, omega_GWO] = VMD(X, alpha_GWO, tau, K_GWO, DC, init, tol);
[u_ABC, u_hat_ABC, omega_ABC] = VMD(X, alpha_ABC, tau, K_ABC, DC, init, tol);
[u_AOA, u_hat_AOA, omega_AOA] = VMD(X, alpha_AOA, tau, K_AOA, DC, init, tol);
[u_ASSA, u_hat_ASSA, omega_ASSA] = VMD(X, alpha_ASSA, tau, K_ASSA, DC, init, tol);
u_hat_PSO=u_hat_PSO';
u_hat_GWO=u_hat_GWO';
u_hat_ABC=u_hat_ABC';
u_hat_AOA=u_hat_AOA';
u_hat_ASSA=u_hat_ASSA';

x2 = zeros(1, num_point); % 预分配空间
t2 = zeros(1, num_point); % 预分配空间——在Matlab中，对于循环内部需要动态增加大小的变量，
x1 = zeros(1, num_point); % 预分配空间
t1 = zeros(1, num_point); % 预分配空间——在Matlab中，对于循环内部需要动态增加大小的变量，
% 频繁的重新分配内存会导致程序的性能下降。因此，为了避免这种情况，我们可以在循环之
% 前先预分配足够大的空间。
%     figure(2);
%     time_fre_fig(K_ASSA,num_point,fs,u_ASSA,u_hat_ASSA,t1);
%     figure(3);
%     time_fre_fig(K_AOA,num_point,fs,u_AOA,u_hat_AOA,t1);
%     figure(4);
%     time_fre_fig(K_PSO,num_point,fs,u_PSO,u_hat_PSO,t1);
%     figure(5);
%     time_fre_fig(K_GWO,num_point,fs,u_GWO,u_hat_GWO,t1);
%     figure(6);
%     time_fre_fig(K_ABC,num_point,fs,u_ABC,u_hat_ABC,t1);

%% 绘制包络谱图
figure(7); % 绘制原始信号的包络谱
h = hilbert(X_Gear); % h是一个复数矩阵
env = abs(h); % env是一个实数矩阵，每一列是一个IMF的包络信号
env_hat = fft(env); % env_hat是一个复数矩阵，每一列是一个包络信号的傅里叶变换
%     subplot(K_SSA,1,1); % 在第k个子图中绘制第k个IMF的包络谱
plot(f,abs(env_hat)); % 绘制包络谱的幅度和频率
xlabel('Frequency (Hz)'); % 设置x轴标签为频率
ylabel('A '); % 设置y轴标签为模式编号
f_hat_X=abs(env_hat(1:1000)); % 包络谱数据
xLowerLimit = 1;  % 设置下限
xUpperLimit =1000 ; % 设置上限
xlim([xLowerLimit,xUpperLimit]);

% 绘制子模态的包络谱图
f_fault_list=[89.05,116.30,144.28];
%     f_fault_list=[66,100,130];
figure(8);
[env_hat_return_ASSA,f_rerurn_ASSA]=es_fig(K_ASSA,u_ASSA);
fcc_ASSA=FCC_rank(env_hat_return_ASSA,f_rerurn_ASSA,f_hat_X,f_fault_list);
xlabel('ASSA','FontSize', 12, 'FontName', 'Arial') % 画出输入信号
figure(10);
[env_hat_return_PSO,f_rerurn_PSO]=es_fig(K_PSO,u_PSO);
fcc_PSO=FCC_rank(env_hat_return_PSO,f_rerurn_PSO,f_hat_X,f_fault_list);
xlabel('PSO','FontSize', 12, 'FontName', 'Arial') % 画出输入信号
figure(11);
[env_hat_return_GWO,f_rerurn_GWO]=es_fig(K_GWO,u_GWO);
fcc_GWO=FCC_rank(env_hat_return_GWO,f_rerurn_GWO,f_hat_X,f_fault_list);
xlabel('GWO','FontSize', 12, 'FontName', 'Arial') % 画出输入信号
figure(12);
[env_hat_return_ABC,f_rerurn_ABC]=es_fig(K_ABC,u_ABC);
fcc_ABC=FCC_rank(env_hat_return_ABC,f_rerurn_ABC,f_hat_X,f_fault_list);
xlabel('ABC','FontSize', 12, 'FontName', 'Arial') % 画出输入信号
figure(13);
[env_hat_return_AOA,f_rerurn_AOA]=es_fig(K_AOA,u_AOA);
fcc_AOA=FCC_rank(env_hat_return_AOA,f_rerurn_AOA,f_hat_X,f_fault_list);
xlabel('AOA','FontSize', 12, 'FontName', 'Arial') % 画出输入信号

% 开始排名
%     rank_list=[fcc_ASSA,fcc_PSO,fcc_GWO,fcc_ABC,fcc_AOA];
rank_list=[fcc_AOA',fcc_ABC',fcc_GWO',fcc_PSO',fcc_ASSA'];
figure(9);
% 设置柱子的颜色

% 定义算法名称和对应的颜色
% algorithm_names = {'f1', 'f2', 'f3'};
algorithm_colors = struct();
% algorithm_colors.f1 = [0.9, 0.7, 0.1];
% algorithm_colors.f2 = [0.3, 0.6, 0.2];
% algorithm_colors.f3 = [0.1, 0.7, 0.9];
algorithm_colors.f1 = [184, 168, 207];
algorithm_colors.f2 = [78, 101, 155];
algorithm_colors.f3 = [253, 207, 158];


h=bar(rank_list', 'grouped');
set(h(1), 'FaceColor', algorithm_colors.f1/255)
set(h(2), 'FaceColor', algorithm_colors.f2/255)
set(h(3), 'FaceColor', algorithm_colors.f3/255)
% 设置柱子的颜色 (使用RGB 255表示)
colors = [
    255, 0, 0;  % 红色
    0, 255, 0;  % 绿色
    0, 0, 255   % 蓝色
    ];

colormap(colors/255);  % 注意要将0-255的RGB值转换为0-1的范围
% for i = 1:5
%     bar(i, rank_list(i), 'FaceColor', algorithm_colors.(algorithm_names{i}));
%     hold on;
% end
% xlabel('Method');
ylabel('FFC');
% title('每个算法找到最小值的次数');
% 设置柱子颜色


% 设置x轴标签
algorithm_names = {'VMD-AOA', 'VMD-ABC', 'VMD-GWO', 'VMD-PSO', 'VMD-ASSA'};
set(gca, 'XTick', 1:5, 'XTickLabel', algorithm_names,'FontSize', 11, 'FontName', 'Times New Roman');
xtickangle(45); % 旋转刻度标签，使其在45度角显示
% 添加图例
legend('Outer ring', 'Rolling element', 'Inner ring', 'FontSize', 10, 'FontName', 'Times New Roman','Location','northwest','LineWidth', 0.7);

%     for i = 1:K_ASSA
%         for j=1:num_point
%         t21(j)=j/fs;
%         end
%             subplot(K_ASSA,2,2*i-1)
%     plot(t21,u_ASSA(i,:));
%         xlabel('Time/s','FontSize', 12, 'FontName', 'Arial') % 画出输入信号
%     ylabel('Amplitude','FontSize', 12, 'FontName', 'Arial'); % 画出输入信号
yLowerLimit = 0;  % 设置下限
yUpperLimit = 0.018; % 设置上限
ylim([yLowerLimit,yUpperLimit]);
%     % 设置图像的宽度和高度（以厘米为单位）
% width_inches = 26;
% height_inches = 4*K_ASSA;
%     set(gcf, 'Units', 'centimeters', 'Position', [0, 0, width_inches, height_inches]);
%     end
%
%     for i = 1:K_SSA
%         N_Gear_u = length(u_hat_SSA(:,i));    % 信号的样本数
%         freq_Gear_u = (0:N_Gear_u-1)*(fs2/N_Gear_u); % 频率向量
%         fft_Gear_u = abs(fft(u_hat_SSA(:,i))); % 取频谱的幅值并进行归一化
%         % 只绘制频谱的一半，并处理幅值
%         half_N_Gear_u = floor(N_Gear_u/2);   % 取频谱的一半
%         fft_signal_half_Gear_u = abs(fft_Gear_u(1:half_N_Gear_u))/N_Gear_u; % 取频谱的幅值并进行归一化
%         % 绘制频域图
%         subplot(K_SSA,1,i)
%         plot(freq_Gear_u(1:half_N_Gear_u), fft_signal_half_Gear_u);
%         % title('信号的频域图');
%                 xlabel('Frequency(Hz)','FontSize', 12, 'FontName', 'Arial');
%         ylabel('Amplitude','FontSize', 12, 'FontName', 'Arial');
%            xLowerLimit = 0;  % 设置下限
%     xUpperLimit = 1000; % 设置上限
%     xlim([xLowerLimit,xUpperLimit]);
%     end

% % 去除直流分量
% u_hat_SSA = detrend(u_hat_SSA);
% u_SSA = detrend(u_SSA);
%% 计算FFC并画柱状图
function fcc_rank_value = FCC_rank(env_hat_return,f_return,f_hat_X,f_fault_list)
% env_hat_return  子模态的包络谱
% f_return  % 横轴频率点，间隔1.2
% f_hat_X   % 原始信号的包络谱
% list_env=zeros(2,1000);
% list_env(1:2,:)=env_hat_return;
% list_env(3,:)=f_return;
fs=1200;
n=1000;
% 计算每个故障频率对应的FCC值，并累加
freq_offset=2;% 允许频率偏移的范围
% 计算频率对应的索引

%     fs * length(u_hat,2)
index_1 = round(f_fault_list(1)*n/fs);
%       index_3=[121,241,361,482,601];
index_2 = round(f_fault_list(2)*n/fs);
index_3 = round(f_fault_list(3)*n/fs);
% index_3 = 121;
% 倍频最高倍数
nf=5;
all_list_1=zeros(nf,size(env_hat_return,1));
all_list_2=zeros(nf,size(env_hat_return,1));
all_list_3=zeros(nf,size(env_hat_return,1));
for j = 1:nf
    index1=j*index_1; % 设置倍频
    index2=j*index_2;
    index3=j*index_3;%
    for i = 1:size(env_hat_return,1)

        % 前后3Hz范围内求和能量
        energy_sum_1=0;
        energy_sum_2=0;
        energy_sum_3=0;
        for f = (index1 - freq_offset) : (index1 + freq_offset)
            energy_sum_1 = energy_sum_1 + abs(env_hat_return(i,f))^2;
        end
        for f = (index2 - freq_offset) : (index2 + freq_offset)
            energy_sum_2 = energy_sum_2 + abs(env_hat_return(i,f))^2;
        end
        for f = (index3 - freq_offset) : (index3 + freq_offset)
            energy_sum_3 = energy_sum_3 + abs(env_hat_return(i,f))^2;
        end
        all_list_1(j,i)=energy_sum_1;
        all_list_2(j,i)=energy_sum_2;
        all_list_3(j,i)=energy_sum_3;
    end

end
% 计算所有模态的故障频率能量和
% all_list
% [total_energy,max_index]=max(all_list);
% max_index
sum_of_max_two_values_1 = 0; % 用于存储最大两个值相加的结果
sum_of_max_two_values_2 = 0; % 用于存储最大两个值相加的结果
sum_of_max_two_values_3 = 0; % 用于存储最大两个值相加的结果
for r=1:size(all_list_1,1)
    sorted_r_1=sort(all_list_1(r,:),'descend');
    sum_of_max_two_values_1 = sum_of_max_two_values_1 + sum(sorted_r_1(1));
end
for r=1:size(all_list_2,1)
    sorted_r_2=sort(all_list_2(r,:),'descend');
    sum_of_max_two_values_2 = sum_of_max_two_values_2 + sum(sorted_r_2(1));
end
for r=1:size(all_list_3,1)
    sorted_r_3=sort(all_list_3(r,:),'descend');
    sum_of_max_two_values_3 = sum_of_max_two_values_3 + sum(sorted_r_3(1));
end
% [total_energy_2,max_indices]=maxk(all_list,2);
% sum_of_max_two_values
total_energy_1=sum_of_max_two_values_1;
total_energy_2=sum_of_max_two_values_2;
total_energy_3=sum_of_max_two_values_3;
% total_energy=max(all_list);
% total_energy=sum(total_energy_2);
% total_energy=sum(all_list);
% 计算整个信号总的能量
total_energy_X = sum(abs(f_hat_X).^2);

% 计算FCC参数
fcc_rank_value(1) = total_energy_1 / total_energy_X;
fcc_rank_value(2) = total_energy_2 / total_energy_X;
fcc_rank_value(3) = total_energy_3 / total_energy_X;

end


%% 函数区
% 1. 画模态分量的时域和频域图
function []= time_fre_fig(K,numpoint,Fs,u,u_hat,t_temp1)
for i = 1:K
    for j=1:numpoint
        t_temp1(j)=j/Fs;
    end
    subplot(K,2,2*i-1)
    plot(t_temp1,u(i,:));
    xlabel('Time/s','FontSize', 12, 'FontName', 'Arial') % 画出输入信号
    ylabel('Amplitude','FontSize', 12, 'FontName', 'Arial'); % 画出输入信号
    xLowerLimit = 0;  % 设置下限
    xUpperLimit = 0.834; % 设置上限
    xlim([xLowerLimit,xUpperLimit]);
    % 设置图像的宽度和高度（以厘米为单位）
    width_inches = 26;
    height_inches = 4*K;
    set(gcf, 'Units', 'centimeters', 'Position', [0, 0, width_inches, height_inches]);
end

% 频谱图
for i = 1:K
    fx = (0:length(u_hat)-1) * (Fs / length(u_hat));
    subplot(K,2,2*K+2-2*i)
    plot(fx,abs(u_hat(i,:)));
    xlabel('Frequency (Hz)','FontSize', 12, 'FontName', 'Arial') % 画出输入信号
    ylabel('Amplitude','FontSize', 12, 'FontName', 'Arial'); % 画出输入信号
    xLowerLimit = 0;  % 设置下限
    xUpperLimit =Fs/2; % 设置上限
    xlim([xLowerLimit,xUpperLimit]);
end
end

%% 绘制包络谱图
% 对每个IMF进行希尔伯特变换
% 对每个包络信号进行傅里叶变换
function [env_hat_return,f_return] = es_fig(K,u)
n = 10000; % n是信号的长度
fs = 12000; % 假设采样频率为1000 Hz
f = (0:n-1)*(fs/n); % f是一个频率向量
env_hat_return=zeros(K,1000);
for i = 1:K
    h(i,:) = hilbert(u(i,:)); % h是一个复数矩阵
    env(i,:) = abs(h(i,:)); % env是一个实数矩阵，每一列是一个IMF的包络信号
    env_hat(i,:) = fft(env(i,:)); % env_hat是一个复数矩阵，每一列是一个包络信号的傅里叶变换
    subplot(K,1,i); % 在第k个子图中绘制第k个IMF的包络谱
    plot(f,abs(env_hat(i,:))); % 绘制包络谱的幅度和频率
    xlabel('Frequency (Hz)'); % 设置x轴标签为频率
    ylabel(['Mode ' num2str(i)]); % 设置y轴标签为模式编号
    xLowerLimit = 1;  % 设置下限
    xUpperLimit =1000 ; % 设置上限
    xlim([xLowerLimit,xUpperLimit]);
    env_hat_return(i,1:1000)=abs(env_hat(i,1:1000)); % 包络谱数据
end
f_return=(0:1000-1)*1.2;
env_hat_return=env_hat_return;


% % 对每个包络信号进行傅里叶变换
% n = size(X_Gear,1); % n是信号的长度
% fs = 12000; % 假设采样频率为1000 Hz
% f = (0:n-1)*(fs/n); % f是一个频率向量
%
%
% % 绘制包络谱图
%
% for k = 1:K_SSA % K是分解出来的IMF的个数
%     subplot(K_SSA,1,k); % 在第k个子图中绘制第k个IMF的包络谱
%     plot(f,abs(env_hat(k,:))); % 绘制包络谱的幅度和频率
%     xlabel('Frequency (Hz)'); % 设置x轴标签为频率
%     ylabel(['Mode ' num2str(k)]); % 设置y轴标签为模式编号
% end

end

%
%
%
%     function h = hilbert(x)
% % HILBERT Compute the Hilbert transform of a real signal.
% %   H = HILBERT(X) computes the Hilbert transform of the real vector X
% %   using the discrete Fourier transform method. X must be a real vector.
%
%     % Check if the input is a real vector
%     if ~isvector(x) || ~isreal(x)
%         error('Input must be a real vector.');
%     end
%
%     % Compute the Fourier transform of the input signal
%     Xf = fft(x);
%
%     % Set the negative frequency components to zero (zero out the imaginary part)
%     n = length(x);
%     Xf((n/2 + 2):end) = 0;
%
%     % Reconstruct the complex signal from the Fourier transform
%     h = ifft(Xf);
%
%     % Ensure that the imaginary part is exactly zero due to numerical errors
%     h = real(h);
%
% end
%
% %     % 假设 modes 是 VMD 分解得到的子模态分量，time 和 freq 是对应的时频域信息
% % Fs=12000;
% % N=10000;
% % num_modes = size(u_SSA, 1);  % 子模态的个数
% % envelopes = zeros(size(u_SSA));  % 用于存储每个子模态的包络
% % time = (0:(N-1)) / Fs;
% %
% % % 计算每个子模态的包络
% % for i = 1:num_modes
% %     envelopes(i, :) = abs(hilbert(u_SSA(i, :)));  % 使用 Hilbert 变换计算包络
% % end
% % % 频率轴
% % freq = u_hat_SSA;
% % % 绘制子模态的包络谱图
% % figure(5);
% % for i = 1:num_modes
% %     subplot(num_modes, 1, i);
% %     plot(freq(i, :), envelopes(i, :));
% %     xlabel('时间');
% %     ylabel('包络');
% %     title(['子模态 ', num2str(i), ' 的包络谱图']);
% % end
%
%
%
%
% %     for i=1:num_point2
% %     %     x1(i)=signal(i+500);
% %           x22(i)=u_hat_SSA(i);
% %         t22(i)=i/fs2; % 采样前500个点
% %     end
% %     subplot(1,2,2)
% %     plot(t22,abs(x22));
% % %     xlabel('Time/s','FontSize', 12, 'FontName', 'Arial','LineWidth', 0.6) % 画出输入信号
% %        % title('信号的频域图');
% %     xlabel('Frequency(Hz)','FontSize', 12, 'FontName', 'Arial','LineWidth', 0.6);
% %     ylabel('Amplitude','FontSize', 12, 'FontName', 'Arial','LineWidth', 0.6);
%
%
%
%
%
% % %%
% % % 齿轮振动信号
% % load 276.mat  %
% % X_Gear = X276_FE_time(20001:30000); % 10000
% %   load ppg.mat % 加载PPG信号
% %     X_ppg_signal=signal(:,1);
% % %     X_ppg=X_ppg_signal(4001:6000);
% %         X_ppg=X_ppg_signal(4001:14000);
% %
% %  Fs = 1000; % 采样频率为1 kHz
% % T = 1; % 信号时长为1秒
% % t = 0:1/Fs:(T-1/Fs); % 时间向量
% %  f1 = 90; %
% % f2 = 150; %
% % f3 = 500; %
% % f4 = 270;
% % x_f1 = sin(2*pi*f1*t);
% % x_f2 = cos(2*pi*f2*t);
% % x_f3 = sin(2*pi*f3*t+cos(2*pi*f4*t));
% % % 生成高斯白噪声
% % SNR_db = 20; % 信噪比为20 dB
% % signal_power = max(x_f1)^2 + max(x_f2)^2 + max(x_f3)^2; % 信号功率
% % noise_power = signal_power / (10^(SNR_db/10)); % 噪声功率
% % noise = sqrt(noise_power) * randn(1, length(t)); % 生成高斯白噪声
% % signal_analog = x_f1 + x_f2 + x_f3+noise;
% %
% % tau=0;
% % DC=0;
% % init=1;
% % tol=1e-7;
% % K_SSA=2;
% % X=X_Gear;
% % % X=X_ppg;
% % [u_SSA, u_hat_SSA, omega_SSA] = VMD(X, 3836, tau, K_SSA, DC, init, tol);
% % % [u_ISSA, u_hat_ISSA, omega_ISSA] = VMD(input_signal, alpha_ISSA, tau, K_ISSA, DC, init, tol);
% % % [u_WOA, u_hat_WOA, omega_WOA] = VMD(input_signal, alpha_WOA, tau, K_WOA, DC, init, tol);
% % % [u_PSO, u_hat_PSO, omega_PSO] = VMD(input_signal, alpha_PSO, tau, K_PSO, DC, init, tol);
% %
% %
% % %% 画子模态和频谱曲线
% % Fs=12000;
% % T=1/Fs;
% % L=length(X); % 10000
% % t=(0:L-1)*T;
% % f=Fs*(0:L/2-1)/L;
% % % 子模态曲线——SSA
% % fig_SSA=figure(3);
% % set(fig_SSA,'name','VMD_SSA')
% % N_SSA=ls(u_SSA, u_hat_SSA, omega_SSA,K_SSA);
% % y_out=zeros();
% % t_out=zeros();
% % for i = 1:K_SSA
% %     for j=1:length(X)
% %         y_out(j)=u_SSA(N_SSA(i),j);% 到这里，将一个模态的2000个数据放进了y_out中，一行2000列
% %         t_out(j)=j/length(X)-1/length(X);
% %     end
% %     subplot(K_SSA,2,2*i-1);
% %     plot(t_out,y_out);
% %     ylabel(['IMF' num2str(i)]) % 添加纵坐标并编号
% %     if i==K_SSA
% %         xlabel('Time/s')
% %     end
% % end
% % %频谱曲线——SSA
% % for i =1:K_SSA
% %     P2=abs(u_hat_SSA(:,N_SSA(i))/L);
% % %     P1=P2(1:L/2+1);
% %     P1=3*P2(L/2+1:end);
% %     P1=3*P1;
% %     subplot(K_SSA,2,2*i)
% %     plot(f,P1);
% %     xlim([0, 500]); % 设置横轴范围为0到6000赫兹
% %     ylabel(['IMF' num2str(i)]) % 添加纵坐标并编号
% %     if i==K_SSA
% %         xlabel('f/Hz')
% %     end
% % end
% % % 画子模态和频谱曲线
% % Fs = 12000;
% % T = 1 / Fs;
% % L = length(X_Gear); % 使用X_Gear，因为之前的X是VMD分解后的结果，不再是原始信号
% % t = (0:L-1) * T;
% % f = Fs * (0:L/2-1) / L;
%
% % % 子模态曲线——SSA
% % fig_SSA = figure(3);
% % set(fig_SSA, 'name', 'VMD_SSA')
% % N_SSA = ls(u_SSA, u_hat_SSA, omega_SSA, K_SSA);
% % y_out = zeros(length(X_Gear), K_SSA); % 修正y_out的初始化
% % t_out = (0:length(X_Gear)-1) * T; % 直接计算t_out的值，无需使用循环
% % for i = 1:K_SSA
% %     y_out(:, i) = u_SSA(N_SSA(i), :);
% %     subplot(K_SSA, 2, 2*i-1);
% %     plot(t_out, y_out(:, i));
% %     ylabel(['IMF' num2str(i)]) % 添加纵坐标并编号
% %     if i == K_SSA
% %         xlabel('Time/s')
% %     end
% % end
%
%
% % 频谱曲线——SSA
% % 假设信号保存在名为u_hat_SSA的矩阵中，其中有10000行和2列
% % u_hat_SSA = [...]; % 在这里填入你的信号数据
%
% % % 信号长度
% % N = size(u_hat_SSA, 1);
% %
% % % 采样频率（假设为1000Hz，如果不清楚采样频率，可以根据实际情况调整这个值）
% % fs = 2000;
% %
% % % 遍历每一列信号，计算频谱并绘制频谱图
% % for col = 1:size(u_hat_SSA, 2)
% %     signal = u_hat_SSA(:, col);
% %
% %     % 计算信号的频谱
% %     frequency_spectrum = fft(signal);
% %     frequency_spectrum = abs(frequency_spectrum(1:N/2+1)); % 取单边频谱
% %
% %     % 计算频率坐标
% %     frequencies = (0:(N/2)) * (fs / N);
% %
% %     % 绘制频谱图
% %     figure;
% %     plot(frequencies, frequency_spectrum);
% %     xlabel('频率 (Hz)');
% %     ylabel('幅度');
% %     title(['信号', num2str(col), '的频谱图']);
% %     grid on;
% % %     xlim([0, 6000]); % 设置横轴范围为0到6000赫兹
% % end
%
%
%
%
%
%
%
%
%
% %% 函数区
%
% % 重心频率
% function F1=gc(f,result_fft)
%     result_fft=abs(result_fft);
%     F1=sum(f.*result_fft)/sum(result_fft);
% end
%
% %  中心频率
% function F2=fc(f,result_fft)
%     result_fft=abs(result_fft);
%     f=f-sum(f.*result_fft)/sum(result_fft);
%     F2=sum(f.*result_fft)/sum(result_fft);
% end
%
% % 傅里叶变换——源代码版
% function [f,result_FFt]=transToFFT(data,fs) %用于求解数据的快速傅里叶变换结果，以便于绘制频谱图
% N=length(data);
% %去均值
% data=data-mean(data);
% %求频率分辨率
% df=fs/(N-1);
% f=(0:N-1)*df;
% Y=fft(data)/N*2;
% result_FFt=abs(Y);
% %取一半
% result_FFt=result_FFt(1:ceil(N/2));
% f=f(1:ceil(N/2));
% end
%
% % 傅里叶变换——chatGPT版
% function [freqs, spectrum] = my_fft(signal, sampling_rate)
% % 计算信号的傅里叶变换并返回频谱
%
% % 采样时间间隔
% dt = 1/sampling_rate; % 1/12000
%
% % 信号长度
% N = length(signal); % 10000
%
% % 频率向量
% freqs = (0:N-1)*(1/dt)/N;
%
% % 对信号进行傅里叶变换并取变换结果的绝对值
% spectrum = abs(fft(signal))/N;
%
% % 对于双边频谱，只需要保留一半的频率和幅度信息
% spectrum = spectrum(1:N/2+1);
% freqs = freqs(1:N/2+1);
% end
%
% % 重心频率排序
% function [N]=ls(u, u_hat, omega, K)
%     U=zeros();
%     N=zeros();
%     for i=1:K
%     % [f1,result_FFt]=transToFFT(u(i,:),20000);
%     [f1,result_FFt]=my_fft(u(i,:),12000);
%     U(i)=gc(f1,result_FFt);
%     %这行MATLAB代码将变量 u 中的第 i 行作为输入信号，调用函数 transToFFT 进行傅里叶变换（FFT），
%     % 并将变换结果存储在变量 result_FFt 中。该函数的第二个参数 20000 表示信号的采样率（单位为 Hz），
%     % 用于将频率轴从正常化频率（单位为周期/样本）转换为实际频率（单位为 Hz）。
%     end
%     U
%     U1=U;
%     for i=1:K
%     [~,n]=min(U1);% m最小值，n最小值的序号
%     N(i)=n;
%     U1(n)=inf;
%     end
%     N
% end