clc;clear all;close all;
% 純凈信號為數(shù)字調(diào)幅信號8ASK，載頻100KHz,波特率1KHz,采樣率600K
x = dmod([0,3,1,3,2,0],100e3,1e3,600e3,'ask',8)';
N = length(x);
%-------------------不隨機(jī)產(chǎn)生的相關(guān)噪聲-----------------------------------
s=awgn(x,3,'measured');
n1=wgn(1,N,1)';
kkk=0.8;
n2(1:(kkk*N))=n1(1:(kkk*N));
n3=wgn(1,N,1)';
n2((kkk*N+1):N)=n3((kkk*N+1):N);
n2=n2';
s=x+n1;
%-------------------隨機(jī)產(chǎn)生的相關(guān)噪聲-------------------------------------
% n1=wgn(1,N,1)';
% kk=300;
% loca=randint(kk,1,[1,N]);
% pp=sort(loca);
% n2=n1;
% n3=wgn(1,N,1)';
% n2(pp)=n3(pp);
% corr2(n1,n2)
% s=x+n1;
%-----用于測試n1與n2有幾個(gè)元素不同---------------------
% aaa=0;
% for iii=1:N
%     if (n1(iii)-n2(iii))
%         aaa=aaa+1;
%     end
% end
%---------------------原始信號信噪比及均方差分析---------------------
snrin=20*log10(norm(x)/norm(s-x));
msein=mmse(s-x);
%--------------------數(shù)據(jù)初始化---------------------------------------
itr=length(x);
m = 4;                  % 濾波器階數(shù)
lmda = 0.99;             % 遺忘因子
cm=35*eye(m);
ww(m)=0;
w=ww';
for i=1:N 
    xxm(m)=0;
    for j=1:1:m
        if i+1-j<=0    % 見自適應(yīng)信號處理，何振亞，P71
            xxm(j)=0;
        else
            xxm(j)=n2(i+1-j);
        end
    end
    xm=xxm';
%---------------------LMS算法----------------------------
    d=s(i);
    u=xm'*cm*xm;
    gm=cm*xm./(lmda+u);
    w=w+gm*(d-xm'*w);
    e(i)=d-xm'*w;
    cm=(cm-gm*xm'*cm)./lmda;
    y(i)=xm'*w;
end
snrout=20*log10(norm(x)/norm(y-x));
mseout=mmse(y-x);
snradd=snrout-snrin
mseadd=msein/mseout
figure,
subplot(3,1,1) ; plot(x) ; title('原始信號');
subplot(3,1,2) ; plot(s) ; title('含噪信號');axis([0,4000,-2,2]);
subplot(3,1,3) ; plot(e) ; title('去噪后信號');axis([0,4000,-1.5,1.5]);