使用LMS滤去正弦干扰

LMS的公式

1. 最陡梯度下降法计算公式
$$H(n+1)=H(n)-\frac{1}{2}\delta V_G(n)\text{\hspace{1em}(1)}$$

其中$\delta$ 取0.4，H(0)=$[3-4{]}^T$
$$V_G(n)=2R_{xx}H(n)-2r_{yx}=2\left[\begin{array}{cc}{r}_{xx}(0)& r_{xx}(1)\\ r_{xx}(1)& r_{xx}(0)\end{array}\right]-2\left[22\begin{array}{c}{r}_{yx}(0)\\ r_{yx}(1)\end{array}\right]\text{\hspace{1em}(2)}$$
$$r_{xx}(k)=\frac{1}{16}\sum _{i=0}15\left(\sqrt{2}sin\frac{2\pi i}{16}\right)\left(\sqrt{2}sin\frac{2\pi (i-k)}{16}\right)=cos\frac{2\pi k}{16}\text{\hspace{1em}(3)}$$
$$r_{yx}(k)=\frac{1}{16}\sum _{i=0}15\left[\sqrt{2}sin\right(\frac{2\pi i}{16}+\frac{\pi }{10}\left)\right]\left(\sqrt{2}sin\frac{2\pi (i-k)}{16}\right)=\frac{1}{\sqrt{2}}cos(\frac{2\pi k}{16}+\frac{\pi }{10})\text{\hspace{1em}(4)}$$

得
$$V_G(n)=2\left[\begin{array}{cc}1& 0.9239\\ 0.9239& 1\end{array}\right]\left[\begin{array}{c}{h}_0(n)\\ h_1(n)\end{array}\right]-2\left[\begin{array}{c}0.6725\\ 0.5377\end{array}\right]$$
$$V_G(0)=\left[\begin{array}{c}-2.7362\\ -3.5320\end{array}\right]$$
2. LMS法计算公式
$$e(n+1)=y(n+1)-H^T(n)X(n+1)$$
$$H(n+1)=H(n)+\delta e(e+1)X(n+1),n=0,1,2\dots$$
其中$\delta$取0.4

误差性能曲面和等值曲线

matlab代码

[h0,h1]=meshgrid(-2:0.1:4,-4:0.1:2);
J=0.55+h0.*h0+h1.*h1+2*cos(pi/8)*h1.*h0-sqrt(2)*h0*cos(pi/10)-sqrt(2)*h1*cos(9*pi/40);%？？？
v=0:0.1:2;
figure;
surf(h0,h1,J);
xlabel('h0');
ylabel('h1');
title('误差性能曲面');
figure;
contour(h0,h1,J,v);
xlabel('h0');
ylabel('h1');
title('误差性能曲面等值曲线');

matlab的白噪声

使用函数mvnrnd(0,0.05,k)‘生成一列长度为K，均值为0，方差为0.05的白噪声
图像如下：

LMS算法程序

matlab代码如下

%LMS算法演示(matlab) 假如我年少有为
%设置参数，N为采样个数，u为步长
close all
clear,clc;
N=16;u=0.4;
%设置迭代次数k
k=1000;
t=1:k;
rk = mvnrnd(0,0.05,k)'; %方差为0.05的白噪声
%figure,plot(1:k,rk);%绘制白噪声
y=rk(t)+sin(2*pi*t/N+pi/10);
x=sqrt(2)*sin(2*pi*t/N);
figure;subplot(3,1,1);
plot(t,rk);
title('白噪声');
subplot(3,1,2);
plot(t,x);
title('相关噪声');
subplot(3,1,3);
plot(t,y);
title('叠加后噪声');
%设置起始权值
wk(1,:)=[3 -4];
%用LMS算法迭代求最佳权值
for i=2:k
    xk(i,:)=[x(i) x(i-1)];%输入信号
    yk(i-1)=xk(i-1,:)*wk(i-1,:)';%输出信号
    dk(i-1)=y(i-1);%期望信号
    err(i-1)=dk(i-1)-yk(i-1);%误差
    wk(i,:)=wk(i-1,:)+u*err(i-1)*xk(i-1,:);%权值迭代
end
Rxx=[[1 0.9239]',[0.9239 1]'];
Ryx=[0.6725 0.5377]';
wk2(:,1)=[3 -4]';
vg=[-2.736 -3.5320]';
for i=1:k-1
    vg=2*Rxx*wk2(:,i)-2*Ryx;
    wk2(:,i+1)=wk2(:,i)-1/2*u*vg;
end;
[x,y]=meshgrid([-2:0.1:4],[-4:0.1:2]);
%求性能表面
%生成曲面值
[h0,h1]=meshgrid(-2:0.1:4,-4:0.1:2);
J=0.55+h0.*h0+h1.*h1+2*cos(pi/8)*h1.*h0-sqrt(2)*h0*cos(pi/10)-sqrt(2)*h1*cos(9*pi/40);%？？？
v=0:0.1:2;
figure;
surf(h0,h1,J);
xlabel('h0');
ylabel('h1');
title('误差性能曲面');
figure;
contour(h0,h1,J,v);
xlabel('h0');
ylabel('h1');
title('误差性能曲面等值曲线');

%画迭代时权值的变化
hold on;plot(wk(:,1),wk(:,2),'r');
hold on;plot(wk2(1,:),wk2(2,:),'g');

%绘制误差与迭代次数的图
figure,plot(err);

多次LMS的统计平均

matlab代码如下

clear all;
close all;

g=100;
N=1000;
q=0.4;
k=2;

jn=zeros(1,N-k+1);
en0=zeros(1,N-k+1);
ejn = zeros(1,N-k+1,g);

for j=1:g
    
    rk = mvnrnd(0,0.05,N)';
    n=1:N;
    y=rk(n)+sin(2*pi*n/16+pi/10);
    x=sqrt(2)*sin(2*pi*n/16);
    
    h=[3 -4]';
    for i=k:N
        u=x(i:-1:i-1);
        en=y(i)-u*h;
        h=h+q*en*u';
        
        en0(i-k+1)=en;%保留最后一次的en
        jn(i-k+1)=en^2;%保留最后一次的jn
        ejn(:,i-k+1,j)=en^2;%保留所有的J(n)
    end
end

%计算统计平均值
for i=1:N-k+1
    jn(i)=sum(ejn(1,i,:))/g;
end

figure;
subplot(2,1,1);
plot(jn);
title('一次的J(n)');
subplot(2,1,2);
plot(en0);
title('一次的e(n)');

figure(2);
plot(jn);
title('100下的平均J(n)');

多次LMS下的平均误差性能曲线

matlab代码如下

clear all;
close all;
N=1000;q=0.4;
k=2;
t=1:N;
g=100;

wk =zeros(2,N,g);
for j=1:g
    rk = mvnrnd(0,0.05,N)'; %方差为0.05的白噪声
    y=rk(t)+sin(2*pi*t/16+pi/10);
    x=sqrt(2)*sin(2*pi*t/16);
    %设置起始权值
    h=[3 -4]';
    wk(:,1)=[3 -4]';
    hm=zeros(2,N);
    hm(:,1)=[-3,4]';
    %用LMS算法迭代求最佳权值
    for i=k:N
        u=x(i:-1:i-1);
        en=y(i)-u*h;
        h=h+q*en*u';
        hm(:,i)=h;
        wk(:,i,j)=h;
    end
    Rxx=[[1 0.9239]',[0.9239 1]'];
    Ryx=[0.6725 0.5377]';
    wk2(:,1)=[3 -4]';
    vg=[-2.736 -3.5320]';
    for i=1:N-1
        vg=2*Rxx*wk2(:,i)-2*Ryx;
        wk2(:,i+1)=wk2(:,i)-1/2*q*vg;
    end
end

for i=1:N-k+1
    wk_e(1,i+1)=sum(wk(1,i+1,:))/g;
    wk_e(2,i+1)=sum(wk(2,i+1,:))/g;
end
wk_e(:,1)=[3,-4]';

[h0,h1]=meshgrid(-2:0.1:4,-4:0.1:2);
J=0.55+h0.*h0+h1.*h1+2*cos(pi/8)*h1.*h0-sqrt(2)*h0*cos(pi/10)-sqrt(2)*h1*cos(9*pi/40);%？？？
v=0:0.1:2;
figure;
surf(h0,h1,J);
xlabel('h0');
ylabel('h1');
title('误差性能曲面');
figure;
contour(h0,h1,J,v);
xlabel('h0');
ylabel('h1');
title('误差性能曲面等值曲线');
hold on
plot(wk_e(1,:),wk_e(2,:),'r');
plot(wk2(1,:),wk2(2,:),'g');
title('100次LMS法平均轨迹');legend('等值曲线','最陡法','平均LMS法');