%dmod也可以產(chǎn)生模擬調(diào)制信號(hào)，見幫助
clear;clc;close all;
%----------------------------------測(cè)試信號(hào)------------------------------------
% %-------------------------BPSK信號(hào),數(shù)字調(diào)制-----------------------------------
% codes=6;                                %碼元個(gè)數(shù)，即輸入調(diào)制信號(hào)的長度
% sigma=1;                                %調(diào)制信號(hào)的幅度
% fs=600e3;                               %采樣率600KHz
% fb=1e3;                                 %波特率1KHz，fb表示對(duì)輸入調(diào)制信號(hào)的采樣率
% fc=100e3;                               %載頻100KHz
% Modulate=2;                             %為選擇調(diào)制方式
% N0=fs/fb;                               %一個(gè)碼元周期內(nèi)的采樣點(diǎn)數(shù)，一個(gè)輸入信號(hào)長度內(nèi)的采樣點(diǎn)數(shù)
% N=N0*codes;                             %總的采樣點(diǎn)數(shù)（已調(diào)信號(hào)的長度）
% p0=pi*rand(1,1);                        %產(chǎn)生初始相位
% symbols=randint(1,codes,[0,1]); 	      %產(chǎn)生基帶碼元
% x_B = ones(N0,1)*symbols;
% x_BB = x_B(:)';                         %根據(jù)波特率要求產(chǎn)生碼元
% signal_base = x_BB;                     %產(chǎn)生基帶信號(hào)
% signal=sigma*dmod(symbols,fc,fb,[fs p0],'psk',Modulate); 
%                                         %產(chǎn)生psk調(diào)制信號(hào)，p0是載頻的初始相位
% x=signal;

%----------------------ASK信號(hào),數(shù)字調(diào)制------------------------------------------
codes=6;                                %碼元個(gè)數(shù)，即輸入調(diào)制信號(hào)的長度
sigma=1;                                %調(diào)制信號(hào)的幅度
fs=600e3;                                 %采樣率1M
fb=1e3;                                 %波特率1KHz，fb表示對(duì)輸入調(diào)制信號(hào)的采樣率
fc=100e3;                               %載頻100KHz
N0=fs/fb;                               %一個(gè)碼元周期內(nèi)的采樣點(diǎn)數(shù)，一個(gè)輸入信號(hào)長度內(nèi)的采樣點(diǎn)數(shù)
N=N0*codes;                             %總的采樣點(diǎn)數(shù)（已調(diào)信號(hào)的長度）
M=4;                                    %調(diào)制方式選擇，MASK
symbols1=linspace(0,1,M);               %對(duì)應(yīng)的M個(gè)幅度
mm=randint(1,codes,[1,M]);
for i=1:codes
    symbols(i)=symbols1(mm(i));         %每個(gè)碼元隨機(jī)取指定值 
end
x=sigma*dmod(symbols,fc,fb,fs,'ask',M); 

%----------------------------FSK信號(hào),數(shù)字調(diào)制------------------------------------
%采用SVD去噪效果都調(diào)頻信號(hào)不太好
% codes=6;                                %碼元個(gè)數(shù)，即輸入調(diào)制信號(hào)的長度
% sigma=1;                                %調(diào)制信號(hào)的幅度
% fs=600e3;                               %采樣率600K
% fb=1e3;                                 %波特率1KHz，fb表示對(duì)輸入調(diào)制信號(hào)的采樣率
% fc=100e3;                               %載頻100KHz
% N0=fs/fb;                               %一個(gè)碼元周期內(nèi)的采樣點(diǎn)數(shù)，一個(gè)輸入信號(hào)長度內(nèi)的采樣點(diǎn)數(shù)
% N=N0*codes;                             %總的采樣點(diǎn)數(shù)（已調(diào)信號(hào)的長度）
% M=4;                                    %調(diào)制方式選擇，fsk的M個(gè)頻率
% symbols1=linspace(0,1,M);               %對(duì)應(yīng)的M個(gè)頻率值
% mm=randint(1,codes,[1,M]);
% for i=1:codes
%     symbols(i)=symbols1(mm(i));         %每個(gè)碼元隨機(jī)取指定值 
% end
% x=sigma*dmod(symbols,fc,fb,fs,'FSK',M);
% % ifr=instfreq(x');                     %測(cè)試頻率
% % ttt=3:N;
% % stem(ttt,ifr);
%-------------------------加入指定強(qiáng)度的噪聲--------------------------------
snr=5;
ps=sum(x.^2)/N;
sigma_noise=abs(sqrt(ps/(10^(snr/10))));
nn=randn(1,N);
enn=sum(nn)/N;                          %隨機(jī)數(shù)nn的均值
nn=nn-enn;                              %使nn均值為0 
noise=sigma_noise*nn;
s=x+noise;
%-------------------------原始信號(hào)信噪比及均方差分析-----------------
snrs=20*log10(norm(x)/norm(s-x));
mmses=mmse(s-x);
%-----------------------分解層數(shù)及構(gòu)成分解矩陣------------------------
%-----------------------矩陣構(gòu)造方法Hankel矩陣-----------------
LL=3;                                    %數(shù)據(jù)分段
NLL=N/LL;                                 %每段的長度
for i=1:LL
    s11(i,:)=s(((i-1)*NLL+1):(i*NLL));
    L=ceil(NLL/2)+1;                          %矩陣行數(shù)
    N1=NLL-L+1;                               %矩陣列數(shù)
    for j=1:L
        a(j,:)=s11(i,j:(N1+j-1));
    end
    MM=3;                                 %特征值個(gè)數(shù)為MM-1
    [u,s1,v]=svd(a);
    s2=zeros(M,N1);
    for k1=MM:L
        s2(k1,:)=s1(k1,:);
    end
    s3=s1-s2;
    %信號(hào)重構(gòu)
    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('純凈信號(hào)');
subplot(312);plot(s); title('含噪信號(hào)');
subplot(313);plot(ys);title('SVD去噪后信號(hào)');
% %-------------------------------SVD去噪后信號(hào)進(jìn)行小波去噪----------------------------------
% %效果不好，信噪比及均方差提高的幅度不大
% yout=wden(sy,'rigrsure','s','sln',6,'db7');
% %-------------------------------小波去噪后信號(hào)信噪比及均方差分析--------------------------
% snryout=20*log10(norm(x)/norm(yout-x))
% mmseyout=mmse(yout-x)
% figure,
% subplot(311);plot(x); title('純凈信號(hào)');
% subplot(312);plot(s); title('含噪信號(hào)');
% subplot(313);plot(yout);title('SVD及小波去噪后信號(hào)');
% %兩種方法誤差分析
% figure,
% subplot(211);plot(yout-x);title('SVD后采用小波');
% subplot(212);plot(sy-x);title('SVD');