function [y,r,vr]=ssa1(x1,L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -----------------------------------------------------------------
% 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. Analisys 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step1 : Build trayectory matrix
N_ssa=length(x1);
if L>N_ssa/2;L=N_ssa-L;end
K=N_ssa-L+1;
X=zeros(L,K);
for i=1:K
X(1:L,i)=x1(i:L+i-1);
end
% Step 2: SVD
S=X*X';
[U,autoval]=eig(S);%求S的全部特征值,构成对角阵autoval,并求特征向量构成U的列向量。
[d,i]=sort(diag(autoval),'descend');%对autoval进行降序排序,i为d在降序前矩阵中的位置
d=-d;
lam=d.^2;
slam=sum(lam);
per=0;
n=1;
for i=1:length(d)-1
per(n)=(lam(i)+lam((i+1)))/slam;%per为贡献率
n=n+1;
end
U=U(:,i);
sev=sum(d);% 特征值的和,可根据占比选择有效信号
% figure(snr_idx)
% SNR_set=-20:2:20;
% snr=SNR_set(snr_idx);
% titlename1 = strcat('Singular Spectrum',num2str(snr),'db');
% plot((d./sev)*100),hold on,plot((d./sev)*100,'rx');
% title(titlename1);xlabel('Eigenvalue Number');ylabel('Eigenvalue (% Norm of trajectory matrix retained)');ylim([0,2]);
V=(X')*U;
rc=U*V';
% Step 3: Grouping
% I=input('Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL] ')
I=1;
Vt=V';
rca=U(:,I)*Vt(I,:);
% figure(snr_idx)
% plot(U(:,I));
% n1=512*(0:1024/4-1)/1024;
% figure(snr_idx)
% subplot 311
% mag11=abs(fft(rca));
% plot(n1,mag11(1:256));
% subplot 312
% mag111=abs(fft(U(:,I)));
% plot(n1,mag111(1:256));
% subplot 313
% mag1111=abs(fft(Vt(I,:)));
% plot(n1,mag1111(1:256));
% Step 4: Reconstruction
y=zeros(N_ssa,1);
Lp=min(L,K);
Kp=max(L,K);
for k=0:Lp-2
for m=1:k+1
y(k+1)=y(k+1)+(1/(k+1))*rca(m,k-m+2);
end
end
for k=Lp-1:Kp-1
for m=1:Lp
y(k+1)=y(k+1)+(1/(Lp))*rca(m,k-m+2);
end
end
for k=Kp:N_ssa-1
for m=k-Kp+2:N_ssa-Kp+1
y(k+1)=y(k+1)+(1/(N_ssa-k))*rca(m,k-m+2);
end
end
% for k=1:Lp-1
% for m=1:k
% y(k)=y(k)+(1/(k))*rca(m,k-m+1);
% end
% end
%
% for k=Lp:Kp-1
% for m=1:Lp
% y(k)=y(k)+(1/(Lp))*rca(m,k-m+1);
% end
% end
%
% for k=Kp+1:N
% for m=k-Kp+1:N-Kp+1
% y(k)=y(k)+(1/(N-k+1))*rca(m,k-m+1);
% end
% end
y1=y';
y=2*y1;
% figure(2)
% mag1=abs(fft(y));
% plot(n1,mag1(1:256));
% figure;subplot(2,1,1);hold on;xlabel('Data poit');ylabel('Original and reconstructed series')
% plot(x1,'b');grid on;plot(y,'r')
% plot(t,y)
r=x1-y1;
% subplot(2,1,2);plot(r,'g');xlabel('Data poit');ylabel('Residual series');grid on
vr=(sum(d(I))/sev)*100;
没有合适的资源?快使用搜索试试~ 我知道了~
MATLAB-多种降噪算法
共10个文件
m:10个
需积分: 2 9 下载量 115 浏览量
2023-03-21
20:33:53
上传
评论 1
收藏 6KB RAR 举报
温馨提示
内含10种降噪算法,包含小波变换、形态滤波、平滑滤波、奇异谱分析、卡尔曼滤波、中值滤波、EMD等等。
资源推荐
资源详情
资源评论
收起资源包目录
降噪算法.rar (10个子文件)
降噪算法
kalman_filter.m 311B
ssa1.m 4KB
emd.m 1KB
xingtailvbo.m 2KB
fushi2.m 592B
pengzhang2.m 587B
wavelet.m 2KB
middle_filte.m 1KB
FFTAnalysis.m 206B
smooth5_3.m 567B
共 10 条
- 1
资源评论
Sou^
- 粉丝: 0
- 资源: 5
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功