LMS.rar_lms均方误差_自适应梯度

%本程序LMS.m利用LMS算法模拟对工频干扰的单频自适应陷波过程 clear;close all f0=50; %模拟50Hz工频干扰 N=2000; %时间采样点个数 M=2048; n=[1:N]; t=[0:N-1]/10000; %采样点间隔1/10000＝0.0001s %产生有用信号及其频谱 fsig=300; %fsig是有用信号频率，一个信号周期1/300=0.0033s s=5*cos(2*pi*fsig*t); %假设有用信号是频率为fsig的余弦波。 figure(1); subplot(2,1,1); plot(t,s); xlabel('时间t(s)');ylabel('有用信号s'); subplot(2,1,2); S=fft(s,M); f=[0:M-1]*10000/M; stem(f,abs(S),'.'); xlabel('频率f(Hz)');ylabel('信号频谱S'); axis([0,900,0,4000]); %产生加入工频干扰的信号及其频谱 f1=[zeros(1,1000),ones(1,N-1000)]; %在第1000个采样点后加入工频波动f1＝1Hz x=s+[5*ones(1,1000),3*ones(1,N-1000)].*sin(2*pi*(f0+f1).*t); %叠加工频干扰的信号x figure(2); subplot(2,1,1); plot(t,x); xlabel('时间t(s)');ylabel('叠加干扰后信号x'); subplot(2,1,2); X=fft(x,M); stem(f,abs(X),'.'); xlabel('频率f(Hz)');ylabel('叠加干扰后信号频谱X'); axis([0,900,0,4000]); %产生参考输入信号及自适应过程初始化 x1=5*cos(2*pi*f0*t+pi/4); %参考输入 x2=5*sin(2*pi*f0*t+pi/4); %参考输入的90度相移 u=0.0003; %抽头权值系数调整步长 w=zeros(2,N+1); %自适应过程中，w记录所有抽头的变更，此处初始化为全0 y(1,1)=0; %自适应滤波的输出信号，提取工频干扰，此处初始化为0 yk=zeros(2,1); %两路LMS自适应输出的记录,初始化为0 e1=zeros(1,N-1); %自适应滤波的误差输出，提取有用信号，初始化为0 %自适应陷波过程 for(i=2:N) d=x(1,i-1); %原始输入 e=d-y(1,i-1);%产生误差信号 x3=e*x1(1,i); x4=e*x2(1,i); w(:,i+1)=w(:,i)+2*u*conj([x3;x4]);%抽头系数更新 yk=w(:,i).*[x1(1,i);x2(1,i)];%产生各路输出信号 y(1,i)=yk(1,1)+yk(2,1);%产生总的输出信号 e1(1,i-1)=e;%误差信号记录 end figure(3); subplot(2,1,1); plot(t(1:N-1),e1); xlabel('时间t(s)');ylabel('误差信号e1'); subplot(2,1,2); E1=fft(e1,M); stem(f,abs(E1),'.'); xlabel('频率f(Hz)');ylabel('误差信号频谱E1'); axis([0,900,0,4000]); figure(4); subplot(2,1,1); plot(t,y); xlabel('时间t(s)');ylabel('提取干扰信号y'); subplot(2,1,2); Y=fft(y,M); stem(f,abs(Y),'.'); xlabel('频率f(Hz)');ylabel('提取干扰信号频谱Y'); axis([0,900,0,4000]); figure(5); plot(t,s,'r'); hold on plot(t(1:N-1),e1); hold off xlabel('时间t(s)');ylabel('有用信号s与误差信号e1对比'); figure(6) plot(t(1:N-1),abs(s(1:N-1)-e1)); xlabel('时间t(s)');ylabel('有用信号与误差信号差|s-e1|'); grid; figure(7); plot(w(1,:),w(2,:)); xlabel('权值w1');ylabel('权值w2'); a=mean(abs(s(500:999)-e1(500:999))); %计算平稳时的误差均值 disp(sprintf('误差均值为： %2.4f',a)); b=std(w(1,500:999));%计算平稳时抽头权值w1标准差 disp(sprintf('w1的标准差为： %2.4f',b)); c=std(w(2,500:999));%计算平稳时抽头权值w2标准差 disp(sprintf('w2的标准差为： %2.4f',c));