clear;clc;

%采樣點(diǎn)數(shù)
 N=25000;

%系統(tǒng)采樣頻率fs=8kHz
Ts=125e-6;

%以下為網(wǎng)絡(luò)響應(yīng)h(n)
n=10;hn=2:n;hn=1./hn;

A=1;
Choice=menu('請(qǐng)選擇信號(hào)形式','　　1KHz正弦信號(hào)　　','　[0,1]均勻分布信號(hào)　','　　　　　EXIT　　　　　');
if Choice==1;
    %輸入為正弦信號(hào)u(n)
    f=1000;
    t=1:N;
    v=A*sin(2*pi*f*t*Ts);
  
elseif Choice==2;
    %輸入為隨機(jī)信號(hào)
    v=rand(1,N);
else
    return;
end;

%濾波器參數(shù)初始化
M=9;%濾波器階數(shù)
step=0.09;%步長(zhǎng)

SNR=0:5:20;
MSE=zeros(length(SNR),N-M+1);
dn_cap=zeros(length(SNR),N-M+1);
ERLE=zeros(length(SNR),N-M+1);
for i=1:length(SNR);
	%加入噪聲
	E=1/2*A^2;
	snr=10.^(SNR(i)/10);
	sigma=sqrt(E/snr);
	randn('state',10);
	Noise=sigma*randn(size(v));
	u=v+Noise;
	%網(wǎng)絡(luò)輸出d(n)
	dn=conv(u,hn);
	%系數(shù)初始值	
	w=zeros(M,1);
	s2=0;e2=0;
	%LMS算法
	for num=M:N;
        un=u(num:-1:num-M+1);
        un=un.';
        dn_cap(i,num)=w'*un;
        en=dn(num)-dn_cap(i,num);
        MSE(i,num-M+1)=en.^2;
        %迭代
        w=w+step.*un.*en;
        %求ERLE 
        s2=s2+dn(num).^2;
        e2=e2+en.^2;
        ERLE(i,num-M+1)=10*log10(s2/e2);
	end;
end;

figure
semilogy(1:20000,abs(MSE(:,1:20000)))
axis([0 20000 10^(-40) 1])
xlabel('n')
ylabel('MSE')
if Choice==1;
    title('輸入為f=1000Hz正弦信號(hào)時(shí)的LMS回波抵消器性能');
else
    title('輸入為[0,1]均勻分布信號(hào)時(shí)的LMS回波抵消器性能');
end;
legend('SNR=0dB','SNR=5dB','SNR=10dB','SNR=15dB','SNR=20dB')
grid on

figure
plot(1:20000,ERLE(:,1:20000))
xlabel('n')
ylabel('ERLE(dB)')
axis([0 20000 10 40])
legend('SNR=0dB','SNR=5dB','SNR=10dB','SNR=15dB','SNR=20dB',4)
grid on

% figure
% plot(u);
% axis([0 500 -1.5 1.5])
% 
% figure
% subplot(2,1,1)
% plot(dn)
% axis([0 500 -1.5 1.5])
% 
% subplot(2,1,2)
% plot(dn_cap)
% axis([0 500 -1.5 1.5])

% figure
% axis([M 500  10^(-3) 1]);