%
clc;
clear all
N = 1500;
n = 1:N;
Fs = 5000;
fb = 950/Fs;
fb0 = 900/Fs;
h = dfilt.delay(200);
w = 10*rand(1,1500);
s = 2*sin(2*pi*fb*n) + 5*sin(2*pi*fb0*n);
% s = 2*sin(2*pi*fb*n);
xn = w + s;
xn_delay = filter(h,xn);
figure
subplot(311);plot(n,s);title('正弦信号');grid
subplot(312);plot(n,w);title('大功率噪声');grid
subplot(313);plot(n,xn);title('输入信号');grid
Nfft = 512;
nfft = 1:Nfft;
xk = abs(fft(xn,Nfft));
f = nfft/Nfft*Fs;
% subplot(211);plot(n,xn);grid
% title('xn-时域图');
% subplot(212),plot(f,xk);grid
% title('xn-幅谱图');
%***********自适应滤波*********
M = 30; %number of the weights
N2 = 500;
n1 = 1:M;
n2 = 1:N2;
u = 0.00005;
k = 1/10*M;
%***********Initialization of w1、e、y*********
w1 = zeros(1,M); %weight
en = zeros(1,N2); %error
y = zeros(1,N2); %output of Adaptive Filtering
%*******************LMS迭代********************
for i = 1:N2
xz(M:-1:1) = xn_delay(i:M+i-1);
y(i) = w1*xz';
en(i) = xn(i+M) - y(i);
w1 = w1+2*u*en(i)*xz;
end
% figure
% subplot(211);plot(n2,en);grid
% title('en-时域图');
ek = abs(fft(en,Nfft));
% subplot(212),plot(f,ek);grid
% title('en-幅谱图');
figure(2)
subplot(211);plot(f,xk);grid
title('xn-幅谱图');ylabel('FFT|x(n)|');
subplot(212);plot(f,ek);grid
title('en-幅谱图');ylabel('FFT|e(n)|');
figure(3)
ASE = k*en.^2;
plot(n2,ASE);grid
title('短时平均误差ASE');
