来源：微信公众号 EW Frontier
摘要
奇异谱分析（SSA）是一种信号分解技术，旨在将信号扩展为可解释的和物理上有意义的分量（例如正弦曲线，噪声等）。本文提出了新的理论和实践结果的分离SSA和介绍一种新的方法称为滑动SSA。首先，SSA与无监督分类算法相结合，提供了一个全自动的数据驱动的组件提取方法，我们调查的局限性，组件分离的理论研究。其次，详细的自动SSA方法被用来设计一种基于滑动分析窗口的方法，该方法在分析具有时变分量数的非平稳信号时，比经典的SSA方法提供更好的结果。最后，建议滑动SSA方法相比，经验模式分解（EMD）和同步压缩短时傅立叶变换（STFT），适用于合成和真实世界的信号。
I.介绍
尽管努力使谱方法对离群值更鲁棒[17]，[18]，但在强调制非平稳信号的情况下，经典SSA方法无法正确重建分量。基于奇异谱分析的分组方法不再有效，因为每个奇异向量只能捕获每个分量的一部分。为了克服这一限制，引入了一种称为滑动SSA的新算法，并提出了一种基于帧的SSA（如短时傅立叶变换（STFT）），在相邻帧之间跟踪所提取的确定性分量。这种新的方法相比，经典的SSA和平稳和非平稳信号提供了显着更好的结果。本文件的结构安排如下。下一节简要回顾了基本的SSA算法，并提出了一种基于层次聚类的无监督组件分组方法。第三节研究了分离能力的SSA考虑几个理论和实际情况下，通过模拟验证。第四节介绍了一种新的算法，称为滑动SSA，它超越了经典SSA算法在处理非平稳多分量信号时的局限性。第五节对四种模式分解技术之间的性能比较作出了贡献：SSA，EMD，同步压缩和新提出的滑动SSA，这是适用于合成和真实世界的信号。最后，在第六部分对本文进行了总结，并对未来的研究进行了展望.
算法流程
II.仿真结果
III.结论
从理论和实践的角度研究了SSA方法，并提出了几种增强建议，用于自动分量分组和非平稳信号处理。一个主要的贡献是引入了一种新的算法，称为滑动SSA，它增强了非平稳信号的分析由于滑动分析窗口的STFT。将该算法与经典的SSA方法以及其他最先进的模式提取方法：EMD方法和同步压缩方法进行了比较。我们的数值实验表明，SSA明显优于EMD和许多情况下，滑动SSA可以获得可比的结果同步压缩方法。在提高新算法的鲁棒性方面，获得了更好的跟踪分量。这可能会导致开发更先进的应用程序，为现实世界的信号。
% Automatic separation between a sinusoid and a chirp merged in a Gaussian white noise using SSA+CAH
%
% Author: D. Fourer and J. Harmouche
% Date: Dec. 2016
% Contact: dominique@fourer.fr
%
% Refs:
% [J. Harmouche, D. Fourer, P. Flandrin, F. Auger and P. Borgnat. One
% or Two Components ? The Singular Spectrum Analysis answers. Proc. SLRA'2015. Grenoble, France.June 2015]
% [J. Harmouche, D. Fourer, P. Flandrin, F. Auger and P. Borgnat. Une ou deux composantes:
% la réponse de l'analyse spectrale singulière. Proc. GRETSI'15. Lyon, France]
clear all
close all
N = 300; %% signal length
L = 40; %% embedded dimension SSA parameter
nc = 2; %% number of components
%% 1: Create signals
time=(1:N)';
s=zeros(N,2);
% s1
lambda0=0.03;
s(:,1) = cos(2*pi*lambda0*time);
s(:,1) = s(:,1) - mean(s(:,1));
% s2
lambda1 = 0.06;
deltalambda = (0.15-0.06);
s(:,2) = 0.8 * cos(2*pi * (lambda1 * (time) + (deltalambda/(2*N)) * (time).^2));
s(:,2) = s(:,2) - mean(s(:,2));
% mix: x= s1 + s2
x = s(:,1)+s(:,2);
% noise with arbitrary std value
noise = 0.05 * randn(N,1);
% mix + noise
xn = x+noise;
%% 2 Components extraction using SSA_auto
epsilon = 1e-3; %% singular spectrum thresholding parameter (increase for more robustness to noise)
Y = ssa_decomp(xn, L, nc, epsilon);
%% 3 Display results
figure
plot(xn)
title(sprintf('observed mixture - SNR=%.2f dB', RQF(x, xn)))
saveas(gcf,'../output/observation.png')
figure
for i = 1:nc
subplot(1,nc, i)
plot(s(:,i))
hold on
plot(Y(:,i), 'r-.')
legend('reference', 'reconstruction')
title(sprintf('component %d - RQF=%.2f dB', i, RQF(s(:,i), Y(:,i))))
end
saveas(gcf,'../output/components.png')function value = RQF(s, s_hat)
% function value = RQF(s, s_hat)
%
% Compute the Reconstruction Quality Function expressed in dB
% and defined by 10 log10(|s|^2 / |s-s_hat|^2)
%
% INPUT:
% s: Reference signal
% s_hat: Estimated signal
%
%
% Author: D. Fourer (dominique@fourer.fr)
% Date: Sept 2015
% Ref: [D.Fourer, F. Auger and P.Flandrin. Recursive versions of the Levenberg-Marquardt
% reassigned spectrogram and of the synchrosqueezed STFT. IEEE Proc. ICASSP 2016.]
value=20*log10(norm(s)/norm(s-s_hat));
endfunction modes = ssa_decomp(x, L, nc, epsilon, force_supk)
% function modes = ssa_decomp(x, L, nc, epsilon, force_supk)
%
% SSA with hierarchical clustering for components grouping
% In the hierarchical clustering, the last components whose energy contribution is relatively
% null are eliminated from the clustering procedure. epsilon serves as threshold.
%
% Input:
% x: input signal
% L: SSA window length
% nc: number of components to extract
% epsilon: threshold applied on the values in the singular spectrum (default: 3e-2)
% force_supk: constraint the number of singular values to use (default: supk depends on epsilon)
%
%
% Authors: J. Harmouche and D. Fourer
% Date: Apr. 2016
%
% Refs:
% [J. Harmouche, D. Fourer, P. Flandrin, F. Auger and P. Borgnat. One
% or Two Components ? The Singular Spectrum Analysis answers. Proc. SLRA'2015. Grenoble, France.June 2015]
% [J. Harmouche, D. Fourer, P. Flandrin, F. Auger and P. Borgnat. Une ou deux composantes:
% la réponse de l'analyse spectrale singulière. Proc. GRETSI'15. Lyon, France]
if ~exist('epsilon', 'var')
epsilon = 3e-3; %0.5; %
end
if ~exist('force_supk', 'var')
force_supk = inf;
end
Nx=length(x);
modes = zeros(nc, 1);
if nc < 1
return;
end
[~,d]=ssa(x,L,1);
% Remove negligible components to improve clustering
supk = length(find( (d/max(d))>epsilon));
if supk < nc
warning('epsilon is too high, trying with eps=%.2f', d(nc)/max(d))
supk = nc;
end
if force_supk < inf
supk = force_supk;
end
y = zeros(Nx,supk);
for I = 1:supk;
[y(:,I),~]=ssa(x,L,I);
end
wcorr = wCorrMat2(x,L,supk);
if ~isempty(wcorr)
Z = linkage(wcorr);
T = cluster(Z,'maxclust', nc);
modes = zeros(Nx, nc);
for i = 1:nc
idx = find(T == i);
if length(idx) > 1
modes(:, i) = sum(y(:, idx).');
else
modes(:, i) = y(:, idx).';
end
end
end
end
function wCorr = wCorrMat2(x,L,supk)
N=length(x);
w=zeros(N,1);
wCorr=zeros(supk,supk);
for i=1 : supk
[Y(:,i),d] = ssa(x,L,i);
end
for i=1:N
w(i)=min([i L N-i+1]);
end
for i =1:supk
for j= 1 : supk
wCorr(i,j)= sum(w.*Y(:,i).*Y(:,j))/(sqrt(sum(w.*Y(:,i).^2)*sum(w.*Y(:,j).^2)));
end
end
endfunction [y,d]=ssa(x,L,I)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -----------------------------------------------------------------
% Author: Francisco Javier Alonso Sanchez e-mail:fjas@unex.es
% Departament of Electronics and Electromecanical Engineering
% Industrial Engineering School
% University of Extremadura
% Badajoz
% Spain
% -----------------------------------------------------------------
%
% SSA generates a trayectory matrix X from the original series x1
% by sliding a window of length L. The trayectory matrix is aproximated
% using Singular Value Decomposition. The last step reconstructs
% the series from the aproximated trayectory matrix. The SSA applications
% include smoothing, filtering, and trend extraction.
% The algorithm used is described in detail in: Golyandina, N., Nekrutkin,
% V., Zhigljavsky, A., 2001. Analysis of Time Series Structure - SSA and
% Related Techniques. Chapman & Hall/CR.
% x1 Original time series (column vector form)
% L Window length
% y Reconstructed time series
% r Residual time series r=x1-y
% vr Relative value of the norm of the approximated trajectory matrix with respect
% to the original trajectory matrix
% The program output is the Singular Spectrum of x1 (must be a column vector),
% using a window length L. You must choose the components be used to reconstruct
% the series in the form [i1,i2:ik,...,iL], based on the Singular Spectrum appearance.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1: Build trajectory matrix
N=length(x);
if L>N/2; L=N-L; end
K=N-L+1;
X=zeros(L,K); % [L,K]=size(X)
for i=1:K
X(1:L,i)=x(i:L+i-1);
end;
% Step 2: SVD
S=X*X';
[EigenVectors,EigenValues]=eig(S);
[d,i]=sort(-diag(EigenValues)); % sort(X) sorts the elements of X in ascending order.
d=-d; EigenVectors=EigenVectors(:,i);
%
%sev=sum(d);
%figure(1); plot(1:L, (d/sev)*100, 1:L, (d/sev)*100, 'rx');
%title('Singular Spectrum');
%xlabel('Eigenvalue Number'); ylabel('Eigenvalue (% Norm of trajectory matrix retained)');
V=(X')*EigenVectors;
% rc=EigenVectors*V';
% Step 3: Grouping
% I=input('Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL]');
Vt=V';
rca=EigenVectors(:,I)*Vt(I,:);
% Step 4: Reconstruction
y=zeros(N,1);
Lp=min(L,K);
Kp=max(L,K);
for k=1:Lp-1
for m=1:k
y(k)=y(k)+rca(m,k-m+1)/k;
end
end
for k=Lp:Kp
for m=1:Lp
y(k)=y(k)+rca(m,k-m+1)/Lp;
end
end
for k=Kp+1:N
for m=k+1-Kp:N+1-Kp
y(k)=y(k)+rca(m,k+1-m)/(N-k+1);
end
end
end
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
