%%
%% 例3.4  
%% 図3.14 MSNアダプティブアンテナの最急降下アルゴリズムによる
%% 初期状態から400イタレーションまでの指向性パターンの変化
%%

% 初期設定
clear all

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

ndata=400;

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));

theta_d=theta(1);
S=exp(-j*twpi*dx.'*sin(theta_d));  % ステアリングベクトル

W(:,1)=S;   % 初期ウエイト
mu=0.002;                    % ステップサイズ

% ウエイトの更新（最急降下法）
for m=1:ndata
    y(m)=W(:,m)'*X(:,m);
    W(:,m+1)=W(:,m)+mu*(S-X(:,m)*conj(y(m)));
end

% 指向性パターンの計算
for m=0:100:400;
    n=round((m-100)/100+2);
    D(n,:)=abs(W(:,m+1)'*exp(-j*twpi*dx.'*sin((-90:0.1:90)*derad))).^2;
    D(n,:)=D(n,:)/max(D(n,:));
end

% グラフプロット
plot(-90:0.1:90,10*log10(D))
axis([-90 90 -60 0]);
set(gca, 'XTick',[-90:30:90])
xlabel('角度(°)');
ylabel('アレー応答値(dB)');
title('MSNアダプティブアレーの最急降下アルゴリズムによる初期状態から400イタレーションまでの指向性パターンの変化')
h=legend('0回','100回','200回','300回','400回',4);
grid on
hold on
