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);
% s=awgn(x,3,'measured');
fs=180e3;
fc1=50;
t=[0:2/fs:0.04];
t=t(1:length(x))';
n1=cos(2*pi*fc1*t);
n2=cos(2*pi*fc1*t+pi/8);
s=x+n1;
% subplot(211);plot(n1);title('1');
% subplot(212);plot(n2);title('2');
%---------------------原始信號信噪比及均方差分析---------------------
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';
%     aaa=xxm           %測試程序
%---------------------RLS算法----------------------------
    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
y=conv(s,w);
y=y(1:length(s));
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(y) ; title('去噪后信號');%axis([0,4000,-1.5,1.5]);
% figure,
% plot();