close all;
clear all;
clc;
% 周期信号的产生
fs=5000;
num=10000;
N=100;
u1=0.005;
% u2=0.001;
% u3=0.007;
n=1:num;
t=n/fs;
f1=50;
f2=100;
f3=250;
f4=500;
f5=750;
f6=1250;
% s=1.0*sin(2*pi*2.7*t)+0.35*sin(2*pi*5.6*t);
s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t)+sin(2*pi*f4*t)+sin(2*pi*f5*t)+sin(2*pi*f6*t);
figure;
subplot(2,1,1);
plot(t,real(s),'k');grid;
ylabel('幅值');
title('输入周期性信号');
% 噪声信号的产生
xn=awgn(s,10);
subplot(2,1,2);
plot(t,xn,'k');grid;
ylabel('幅值');
xlabel('时间');
title('随机噪声信号');
% 信号滤波
y1=zeros(1,num);
y1(1:N)=xn(1:N);
w1=zeros(1,N);
% y2=zeros(1,num);
% y2(1:N)=xn(1:N);
% w2=zeros(1,N);
% y3=zeros(1,num);
% y3(1:N)=xn(1:N);
% w3=zeros(1,N);
% % e=zeros(1:num) ; % 预期结果序列
% rho_max = max(eig(xn*xn.')); % 输入信号相关矩阵的最大特征值
% mu = rand()*(1/rho_max) ; % 收敛因子 0 < mu < 1/rho
for i=(N+1):num
XN=xn((i-N+1):i);
y1(i)=w1*XN';
e1(i)=s(i)-y1(i);
w1=w1+u1*e1(i)*XN;
% y2(i)=w2*XN';
% e2(i)=s(i)-y2(i);
% w2=w2+u2*e2(i)*XN;
% y3(i)=w3*XN';
% e3(i)=s(i)-y3(i);
% w3=w3+u3*e3(i)*XN;
end
% 绘制滤波器输入信号
figure(2);
% subplot(2,1,2);
plot(t,real(y1),'k');
% grid;hold on;
% plot(t,real(y2),'-.k');hold on;
% plot(t,real(y3),'k*');legend('u1=0.001','u2=0.005','u3=0.007');
ylabel('幅值');
xlabel('时间/s');
title('自适应滤波器输出信号');
% for b=1:num-N;
% bi(b)=sum(pp(:,b))/g;
% end
% 绘制自适应滤波器输出信号,预期输出信号和两者的误差
figure(3)
% t=1:num-N;
% plot(10*log10(bi));
plot(abs(e1),'k');
% grid;hold on;
% plot(abs(e2),'k-.');hold on;
% plot(abs(e3),'k*');legend('u1=0.001','u2=0.005','u3=0.007');
ylabel('幅值');
xlabel('频率/Hz');
title('误差信号');
nfft=1024;
window=hamming(1024);
noverlap=256;
range='oneside';
[Pxx1,f]=pwelch(y1,window,noverlap,nfft,5000,range);
% [Pxx2,f]=pwelch(y2,window,noverlap,nfft,100,range);
% [Pxx3,f]=pwelch(y3,window,noverlap,nfft,100,range);
plot_Pxx1=10*log10(Pxx1);
% plot_Pxx2=10*log10(Pxx2);
% plot_Pxx3=10*log10(Pxx3);
% subplot(222);
figure(4);
plot(f,plot_Pxx1/10,'k');
% hold on;
% plot(f,plot_Pxx2/10,'k-.');hold on;plot(f,plot_Pxx3/10,'k*');legend('u1=0.001','u2=0.005','u3=0.007');
% set(gca,'Xtick',[0,1,2,3,4,5,6,7,8,9,10]);
% set(gca,'Ytick',[0.02,0.04,0.06,0.08]);
% axis([0,10,-inf,inf]);
xlabel('频率/Hz');
ylabel('功率谱密度');
title('PSD');
% set(gca,'FontSize',15)