lms.m

自适应信号处理

function [yn,w,e]=lms(xn,dn,belta,M,N,w0) 
     %输入，输出均为列向量
    n=length(xn);
    w=w0;%N为叠代次数
    e=zeros(n,1);
    yn=ones(n,1);
  if N>length(xn)
      error('迭代次数太大');
  elseif N<M
      error('迭代次数太小');
  end
  
  
    %R=max(eig(x*x'));
    %belta=rand()*(1/R);
         
          for k=M:N
             x0=xn(k:-1:k-M+1);
             y0=w(:,k-1)'*x0;
             e(k)=dn(k)-y0;
             w(:,k)=w(:,k-1)+2*e(k)*belta*x0;
          end   
          
          for k=M:length(xn)
              x1=xn(k:-1:k-M+1);
              yn(k)=w(:,end)'*x1;
          end
       clc;
clear;
clear all;
clf;
SN=40;
M=16;
N=1000;
t=linspace(0,100,N);
f0=5;
%d=2*sin(0.5*t);
phase=sin(0.25*pi)*cos(0.3*pi);

%d=exp(j*2*pi*f0*t+j*phase);
%d=real(d);


for k=1:N;%N
 sum=0;
for i=1:M%M天线元数目
sum=sum+cos(2*pi*f0*t(k))*exp(j*(i-1)*pi*phase);
end
dn(k)=real(sum);
end

%噪声幅度控制
sp=mean(dn)^2;%信号功率
np=sp*10^(-SN/20);%noise power
figure(1);
subplot(2,2,1);
plot(t,dn); 
title('输入信号');

%同相的高斯噪声

randn('state',randint(1,1)*5000000);
n1=randn(size(t));
n1=np.*exp(j*pi*phase*n1);
subplot(2,2,2);
plot(t,n1);title('同相噪声信号');


%其他方向的噪声

rand('state',randint(1,1)*5000000);
n2=rand(1,N);
alpha=rand(1,N)*pi;
beta=rand(1,N)*2*pi;
%[alpha,beta]=meshgrid(alpha,beta);
n2=n2.*cos(alpha).*sin(beta);
n2=np.*exp(j*pi*n2);

subplot(2,2,3);
plot(t,n2);title('各相噪声信号');

%合成输入信号
xn=dn+n1+n2;

xn=xn';
dn=dn';
R=max(eig(xn*xn'));
rand('state',randint(1,1)*10000)
belta=rand()*(1/R);
[yn,w,e]=lms(xn,dn,belta,M,N);
W=w(:,end)'
figure(2);
plot(t,xn,'b',t,yn,'r');title('输入输出信号');legend('xn','yn');
figure(3);
subplot(2,1,1);
plot(t,yn-dn,'r');title('误差曲线d-y');
subplot(2,1,2);plot(t,e);title('学习曲线');