%%
%% 例3.2  
%% 図3.4 LMSアルゴリズムによる2乗誤差の変化
%%

% 初期設定
clear all

% 入力パラメータ設定
derad=pi/180;   % deg->rad
twpi=2*pi; 

ndata=500;

kelm=2;                              % 素子数
dspc=0.5;                            % 素子間隔(*λ)
dx=linspace(0,dspc*(kelm-1),kelm);   % 素子位置

theta = [30,  -60]*derad;    % 所望波（第一成分），不要波の到来角
Pw    = [1.0, 1.0];          % 所望波（第一成分），不要波の入力電力
Pn=0.01;                     % 雑音電力
L = length(theta);           % 到来波数

V=exp(-j*twpi*dx.'*sin(theta));  % 所望波，不要波のアレー応答ベクトル
% V = [Vs, Vu1, Vu2,...]

rand('state',0)                 
randn('state',0)

for iL=1: L
    pha=rand(1,ndata)*twpi;
    sv(iL,:) = sqrt(Pw(iL))*exp(j*pha);
end

% 入力ベクトルのサンプル値：データ行列 (kelm x ndata)
X=V*sv+sqrt(Pn/2)*(randn(kelm,ndata)+j*randn(kelm,ndata));

% 相関行列の理論値計算
Rxx=V*diag(Pw)*V'+Pn*eye(kelm);

% 相関ベクトルの理論値計算
rxr=Pw(1)*exp(-j*twpi*dx.'*sin(theta(1)));

W(:,1)=[1 zeros(1,kelm-1)].';         % 初期ウエイト

for mu=0.01:0.02:0.07   % ステップサイズ

% ウエイトと２乗誤差の計算 最急降下法(LMSアルゴリズム)
        for m=1:ndata
            e(m)=sv(1,m)-W(:,m)'*X(:,m); 
            W(:,m+1)=W(:,m)+mu*(X(:,m)*conj(e(m)));
            n=round((mu-0.01)/0.02+1);
            E(n,m)=Pw(1)-W(:,m).'*conj(rxr)-W(:,m)'*rxr+W(:,m)'*Rxx*W(:,m);
        end
end

% グラフプロット
semilogy(1:ndata,abs(E))
axis([0 ndata 0.001 10])
xlabel('イタレーション回数');
ylabel('E[|e(m)|^2]');
title('LMSアルゴリズムによる2乗誤差の変化')
h=legend('\mu=0.01','\mu=0.03','\mu=0.05','\mu=0.07');
grid on
hold on
