%%
%% 例3.7  
%% 図3.29 SD-CMA(p=q=2)による出力包絡線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, 0.5];          % 所望波（第一成分），不要波の入力電力
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));

W(:,1)=[1 zeros(1,kelm-1)].';         % 初期ウエイト
mu=0.005;                             % ステップサイズ
sigma=1;

% ウエイトの計算：最急降下法(2-2 CMA)
        for m=1:ndata
          y(m)=W(:,m)'*X(:,m); 
          W(:,m+1)=W(:,m)-4*mu*X(:,m)*conj(y(m))*(abs(y(m))^2-sigma^2);
            % 出力包絡線2乗誤差の計算
            ee(m)=0;
            k=m;
            while     (k>0) & (k>m-20)
                ee(m)=ee(m)+abs(abs(y(k))^2-sigma^2)^2;
                k=k-1;
            end    
            ee(m)=ee(m)/(m-k);
        end
        
% グラフプロット
semilogy(1:ndata,abs(ee))
axis([0 ndata 0.001 10])
xlabel('イタレーション回数');
ylabel('mean(||y(m)^2|-1|^2)');
title('SD-CMA(p=q=2)による出力包絡線2乗誤差の変化')
grid on
hold on