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('學習曲線');