svdexp1.m

射频信号的去噪

%dmod也可以产生模拟调制信号，见帮助
clear;clc;close all;
%----------------------------------测试信号------------------------------------
% %-------------------------BPSK信号,数字调制-----------------------------------
% codes=6;                                %码元个数，即输入调制信号的长度
% sigma=1;                                %调制信号的幅度
% fs=600e3;                               %采样率600KHz
% fb=1e3;                                 %波特率1KHz，fb表示对输入调制信号的采样率
% fc=100e3;                               %载频100KHz
% Modulate=2;                             %为选择调制方式
% N0=fs/fb;                               %一个码元周期内的采样点数，一个输入信号长度内的采样点数
% N=N0*codes;                             %总的采样点数（已调信号的长度）
% p0=pi*rand(1,1);                        %产生初始相位
% symbols=randint(1,codes,[0,1]); 	      %产生基带码元
% x_B = ones(N0,1)*symbols;
% x_BB = x_B(:)';                         %根据波特率要求产生码元
% signal_base = x_BB;                     %产生基带信号
% signal=sigma*dmod(symbols,fc,fb,[fs p0],'psk',Modulate); 
%                                         %产生psk调制信号，p0是载频的初始相位
% x=signal;

%----------------------ASK信号,数字调制------------------------------------------
codes=6;                                %码元个数，即输入调制信号的长度
sigma=1;                                %调制信号的幅度
fs=600e3;                                 %采样率1M
fb=1e3;                                 %波特率1KHz，fb表示对输入调制信号的采样率
fc=100e3;                               %载频100KHz
N0=fs/fb;                               %一个码元周期内的采样点数，一个输入信号长度内的采样点数
N=N0*codes;                             %总的采样点数（已调信号的长度）
M=4;                                    %调制方式选择，MASK
symbols1=linspace(0,1,M);               %对应的M个幅度
mm=randint(1,codes,[1,M]);
for i=1:codes
    symbols(i)=symbols1(mm(i));         %每个码元随机取指定值 
end
x=sigma*dmod(symbols,fc,fb,fs,'ask',M); 

%----------------------------FSK信号,数字调制------------------------------------
%采用SVD去噪效果都调频信号不太好
% codes=6;                                %码元个数，即输入调制信号的长度
% sigma=1;                                %调制信号的幅度
% fs=600e3;                               %采样率600K
% fb=1e3;                                 %波特率1KHz，fb表示对输入调制信号的采样率
% fc=100e3;                               %载频100KHz
% N0=fs/fb;                               %一个码元周期内的采样点数，一个输入信号长度内的采样点数
% N=N0*codes;                             %总的采样点数（已调信号的长度）
% M=4;                                    %调制方式选择，fsk的M个频率
% symbols1=linspace(0,1,M);               %对应的M个频率值
% mm=randint(1,codes,[1,M]);
% for i=1:codes
%     symbols(i)=symbols1(mm(i));         %每个码元随机取指定值 
% end
% x=sigma*dmod(symbols,fc,fb,fs,'FSK',M);
% % ifr=instfreq(x');                     %测试频率
% % ttt=3:N;
% % stem(ttt,ifr);
%-------------------------加入指定强度的噪声--------------------------------
snr=5;
ps=sum(x.^2)/N;
sigma_noise=abs(sqrt(ps/(10^(snr/10))));
nn=randn(1,N);
enn=sum(nn)/N;                          %随机数nn的均值
nn=nn-enn;                              %使nn均值为0 
noise=sigma_noise*nn;
s=x+noise;
%-------------------------原始信号信噪比及均方差分析-----------------
snrs=20*log10(norm(x)/norm(s-x));
mmses=mmse(s-x);
%-----------------------分解层数及构成分解矩阵------------------------
%-----------------------矩阵构造方法Hankel矩阵-----------------
LL=3;                                    %数据分段
NLL=N/LL;                                 %每段的长度
for i=1:LL
    s11(i,:)=s(((i-1)*NLL+1):(i*NLL));
    L=ceil(NLL/2)+1;                          %矩阵行数
    N1=NLL-L+1;                               %矩阵列数
    for j=1:L
        a(j,:)=s11(i,j:(N1+j-1));
    end
    MM=3;                                 %特征值个数为MM-1
    [u,s1,v]=svd(a);
    s2=zeros(M,N1);
    for k1=MM:L
        s2(k1,:)=s1(k1,:);
    end
    s3=s1-s2;
    %信号重构
    anew=u*s3*v';
    ys1(i,:)=[anew(1,:),anew((2:L),N1)'];
end
for i=1:LL
    ys(((i-1)*NLL+1):(i*NLL))=ys1(i,:);
end

%----------------------SVD去噪后，信噪比及均方差分析--------------------------------
snry=20*log10(norm(x)/norm(ys-x))
mmses
mmsey=mmse(ys-x)
subplot(311);plot(x); title('纯净信号');
subplot(312);plot(s); title('含噪信号');
subplot(313);plot(ys);title('SVD去噪后信号');
% %-------------------------------SVD去噪后信号进行小波去噪----------------------------------
% %效果不好，信噪比及均方差提高的幅度不大
% yout=wden(sy,'rigrsure','s','sln',6,'db7');
% %-------------------------------小波去噪后信号信噪比及均方差分析--------------------------
% snryout=20*log10(norm(x)/norm(yout-x))
% mmseyout=mmse(yout-x)
% figure,
% subplot(311);plot(x); title('纯净信号');
% subplot(312);plot(s); title('含噪信号');
% subplot(313);plot(yout);title('SVD及小波去噪后信号');
% %两种方法误差分析
% figure,
% subplot(211);plot(yout-x);title('SVD后采用小波');
% subplot(212);plot(sy-x);title('SVD');