%本程序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));
