【lssvm回归预测】基于变模态结合秃鹰算法优化最小二乘支持向量机VMD-BES-LSSVM实现数据预测模型（多输入单输出）含Matlab代码

% test-script for VMD
% authors: Dominique Zosso and Konstantin Dragomiretskiy
% zosso@math.ucla.edu --- http://www.math.ucla.edu/~zosso
% Initial release 2013-12-12 (c) 2013
%
% When using this code, please do cite our paper:
% -----------------------------------------------
% K. Dragomiretskiy, D. Zosso, Variational Mode Decomposition, IEEE Trans.
% on Signal Processing (in press)
% please check here for update reference: 
%          http://dx.doi.org/10.1109/TSP.2013.2288675

%--------------- Preparation
clear all;
close all;
clc;

% Time Domain 0 to T
fs=400;
T=60;
nn=fs*T;
t=1/fs:1/fs:T;
freqs = ((1/T:1/T:fs)-fs/2-1/T);
% nn = 1000;
% T=1;
% fs = nn/T;
% t = (1:nn)/fs;
% freqs = ((1/T:1/T:fs)-fs/2-1/T);
% freqs = 2*pi*(t-0.5-1/nn)*T*nn;

% center frequencies of components
f_1 = 20;
f_2 = 60;
f_3 = 100;
f_4 = 50;

% modes
v_1 = (cos(2*pi*f_1*t));
v_2 = 1/4*(cos(2*pi*f_2*t));
v_3 = 1/16*(cos(2*pi*f_3*t));
v_4=1/3*(cos(2*pi*f_4*t));
% 
% for visualization purposes
% fsub = {};
% wsub = {};
% fsub{1} = v_1;
% fsub{2} = v_2;
% fsub{3} = v_3;
% wsub{1} = 2*pi*f_1;
% wsub{2} = 2*pi*f_2;
% wsub{3} = 2*pi*f_3;

% composite signal, including noise
% f=csvread('17.txt',0,0,[0 0 nn-1 0]);
f = v_1 + v_2 + v_3 + v_4+0.1*randn(size(v_1));
f_hat = fftshift((fft(f)));

% some sample parameters for VMD
alpha = 2000;        % moderate bandwidth constraint
tau = 0.001;            % noise-tolerance (no strict fidelity enforcement)
K = 4;              % 3 modes
DC = 0;             % no DC part imposed
init = 1;           % initialize omegas uniformly
tol = 1e-7;






%--------------- Run actual VMD code

[u, u_hat, omega] = VMD(f, alpha, tau, K, DC, init, tol);

omega=omega*fs;










%--------------- Visualization

% For convenience here: Order omegas increasingly and reindex u/u_hat
[~, sortIndex] = sort(omega(end,:));
omega = omega(:,sortIndex);
u_hat = u_hat(:,sortIndex);
u = u(sortIndex,:);
linestyles = {'b', 'g', 'm', 'c', 'c', 'r', 'k'};

figure('Name', 'Composite input signal' );
plot(t,f, 'k');
set(gca, 'XLim', [0 1]);

% for sub = 1:length(fsub)
%     figure('Name', ['Input signal component ' num2str(sub)] );
%     plot(t,fsub{sub}, 'k');
%     set(gca, 'XLim', [0 1]);
% end

figure('Name', 'Input signal spectrum' );
plot(freqs(nn/2+1:end), abs(f_hat(nn/2+1:end)), 'k');
% loglog(freqs(nn/2+1:end), abs(f_hat(nn/2+1:end)), 'k');
% set(gca, 'XLim', [1 nn/2]*pi*2, 'XGrid', 'on', 'YGrid', 'on', 'XMinorGrid', 'off', 'YMinorGrid', 'off');
% ylims = get(gca, 'YLim');
% hold on;
% for sub = 1:length(wsub)
%     plot([wsub{sub} wsub{sub}], 'k--');
% %     loglog([wsub{sub} wsub{sub}], ylims, 'k--');
% end
% set(gca, 'YLim', ylims);

figure('Name', 'Evolution of center frequencies omega');
for k=1:K
    plot(omega(:,k), 1:size(omega,1), linestyles{k});
%      semilogx(2*pi/fs*omega(:,k), 1:size(omega,1), linestyles{k});
    hold on;
end
% set(gca, 'YLim', [1,size(omega,1)]);
% set(gca, 'XLim', [2*pi,0.5*2*pi/fs], 'XGrid', 'on', 'XMinorGrid', 'on');

figure('Name', 'Spectral decomposition');
plot(freqs(nn/2+1:end), abs(f_hat(nn/2+1:end)), 'k:');
% loglog(freqs(nn/2+1:end), abs(f_hat(nn/2+1:end)), 'k:');
% set(gca, 'XLim', [1 nn/2]*pi*2, 'XGrid', 'on', 'YGrid', 'on', 'XMinorGrid', 'off', 'YMinor服务器托管网Grid', 'off');
hold on;
for k = 1:K
    plot(freqs(nn/2+1:end), abs(u_hat(nn/2+1:end,k)), linestyles{k});
%     loglog(freqs(nn/2+1:end), abs(u_hat(nn/2+1:end,k)), linestyles{k});
end
% set(gca, 'YLim', ylims);


for k = 1:K
    figure('Name', ['Reconstructed mode ' num2str(k)]);
    plot(t,u(k,:), linestyles{k});   hold on;
%     if ~isempty(fsub)
%         plot(t, fsub{min(k,length(fsub))}, 'k:');
%     end
    set(gca, 'XLim', [0 1]);
end