使用奇异值分解对周期信号进行去噪

奇异值很多用来在预测系统上，感觉上是线性代数上的AU分解，不过高明得多，而且奇异值在周期信号的去噪效果上效果显著，我现在写的就是奇异值分解在周期信号上的应用，主体代码是一个师兄给我的，我对代码自己搞懂了再进行了部分修改，代码如下：

%==== 输入信号 ======
N=300; %生成300个点的信号
t=0:30*pi/N:30*pi;t(end)=[];

s=cos(0.1*pi*t)+sin(0.3*pi*t)+cos(0.5*pi*t)+sin(0.7*pi*t); % 原始信号
y=s+randn(1,N);

subplot(4,1,1); plot(s,’LineWidth’,1.5)
title(‘原始模拟信号 ‘,’FontSize’,16)
set(gca,’box’,’off’)
%set(gca,’linewidth’,1.5);
subplot(4,1,2); plot(y,’LineWidth’,1.5)
title(‘原始模拟信号+高斯噪声’,’FontSize’,16)
set(gca,’box’,’off’)
%set(gca,’linewidth’,1.5);
%=============================
%===== 傅里叶变换结果 ====
NFFT = 2^nextpow2(N); % Next power of 2 from length of y
Y = fft(y,NFFT)/N;
f = N/2*linspace(0,1,NFFT/2+1);
subplot(4,1,3);plot(f,2*abs(Y(1:NFFT/2+1)),’LineWidth’,1.5) %查看信号的频谱，因为奇异值去噪一般选择2*信号主频的个数
title(‘频率 ‘,’FontSize’,16)
set(gca,’box’,’off’)
%set(gca,’linewidth’,1.5);

n=8;

%=============================
%==== 奇异值分解 ====
for i=1:N/2+1
t=i:i+N/2-1;
for j=1:N/2
A(j,i)=y(t(j)); %把一维信号重构为矩阵做奇异值分解
end
end
[U,S,V] = svd(A);

%=============================
%==== 重构信号 =====
X=zeros(size(A));
for i=1:n; %选取前n个大奇异值
X=X+U(:,i)*S(i,i)*V(:,i)’;
end
JG=zeros(1,N);
for i=1:N
a=0;m=0;
for j1=1:N/2
for j2=1:N/2+1
if i+1==j1+j2
a=a+X(j1,j2); %把矩阵重构回一维信号
m=m+1;
end
end
end
JG(i)=a/m;
end
% JG(N/2+1:end)=X(:,end);
subplot(4,1,4);plot(JG,’LineWidth’,1.5)
title(‘奇异值降噪信号’,’FontSize’,16)
xlabel(‘时间’,’FontSize’,16)
set(gca,’box’,’off’)
%set(gca,’linewidth’,1.5);

set(gcf,’color’,’w’)

结果如下：

具体为什么要把一维信号变换成矩阵A，可以看一下这个文献，讲得还是比较清楚的：https://wenku.baidu.com/view/f75ffefcfab069dc50220132.html?qq-pf-to=pcqq.group
奇异值分解最好运算不要超过3000以上，因为运算量太大了，我的机子算大数据量（2600以上）的奇异值分解要卡好一会。