%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Energy Compaction 評価プログラム I (一括変換版)
%
% いくつかの直交変換に対する Energy Compaction 特性 (少ない変換係数でどれだ
% けもとの信号を復元できるかを表す尺度) を、相関係数を変化させながら比較。
%
%                                                      2003.11.26 Katto Lab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

% 相関係数
p = 0.95;

% 試行回数
LOOP = 20;

% 入力信号ベクトル
n = 512;
nn = 1:n;
x = zeros(1,n);

% 変換係数ベクトル (DFT, HT, DCT, KLT)
XDFT = zeros(1,n);
XHT  = zeros(1,n);
XDCT = zeros(1,n);
XKLT = zeros(1,n);

% Energy Compaction 特性の格納ベクトル
ADFT = zeros(1,n);
AHT  = zeros(1,n);
ADCT = zeros(1,n);
AKLT = zeros(1,n);

% Energy Compaction 特性の累積結果の格納ベクトル
RDFT = zeros(1,n);
RHT  = zeros(1,n);
RDCT = zeros(1,n);
RKLT = zeros(1,n);

% LOOP 回数のシミュレーション試行
for loop=1:LOOP
    message = sprintf('LOOP: %d', loop);
    disp(message);
    
    % AR(1)モデルに従った入力信号の生成
    % x(n) = p * x(n-1) + z(n)
    % x(n): 入力信号
    % p:    相関係数
    % z(n): 平均0、分散1の乱数
    for i=1:n
        if(i==1) x(i) = randn;
        else x(i) = p * x(i-1) + randn;
        end;
    end;

    % 各種直交変換
    % 入力信号を直交変換し、変換係数ベクトルに格納
    %% DFT
    %% 変換係数の複素共役の関係を考慮して、実数、虚数の順番に並び替え
    xdft = (1.0/sqrt(n))*fft(x);
    XDFT(1) = xdft(1);
    XDFT(n) = xdft(n/2+1);
    for i=2:n/2
        XDFT(2*i-2) = sqrt(2) * real(xdft(i));
        XDFT(2*i-1) = sqrt(2) * imag(xdft(i));
    end;
    %% HT
    XHT = (1.0/sqrt(n)) * (x*hadamard(n));
    %% DCT
    XDCT = dct(x);
    %% KLT
    [R,RX] = corrmtx(x,n-1);    % 入力信号の自己相関行列の作成
    RX = (1.0/RX(1,1)) * RX;    % 対角要素を1に正規化 (不要)
    [L,lambda] = eig(RX);       % 固有値展開 (L: 固有ベクトル行列、lambda: 固有値ベクトル)
    L = fliplr(L);              % 固有値の大きな順番に固有ベクトルを並べ替え
    XKLT = x*L;
    
    % 変換係数の絶対値の大きな順番に並び替え
    %XDFT = fliplr(sort(abs(XDFT)));
    %XHT  = fliplr(sort(abs(XHT)));
    %XDCT = fliplr(sort(abs(XDCT)));
    %XKLT = fliplr(sort(abs(XKLT)));

    % Energy Compaction 特性の評価
    %% 変換係数の二乗値の累積和
    ADFT(1) = XDFT(1)*XDFT(1);
    AHT(1)  = XHT(1)*XHT(1);
    ADCT(1) = XDCT(1)*XDCT(1);
    AKLT(1) = XKLT(1)*XKLT(1);
    for i=2:n
        ADFT(i) = ADFT(i-1) + XDFT(i)*XDFT(i);
        AHT(i)  = AHT(i-1)  + XHT(i)*XHT(i);
        ADCT(i) = ADCT(i-1) + XDCT(i)*XDCT(i);
        AKLT(i) = AKLT(i-1) + XKLT(i)*XKLT(i);
    end;
    %% 0～1になるように正規化
    for i=1:n
        ADFT(i) = ADFT(i) / ADFT(n);
        AHT(i)  = AHT(i) / AHT(n);
        ADCT(i) = ADCT(i) / ADCT(n);
        AKLT(i) = AKLT(i) / AKLT(n);
    end;
    
    % Energy Compaction 特性の累積評価結果の格納
    for i=1:n
        RDFT(i) = RDFT(i) + ADFT(i);
        RHT(i)  = RHT(i) + AHT(i);
        RDCT(i) = RDCT(i) + ADCT(i);
        RKLT(i) = RKLT(i) + AKLT(i);
    end;
end;

% Energy Compaction 特性の累積評価結果が0～1になるように正規化
for i=1:n
    RDFT(i) = RDFT(i) / LOOP;
    RHT(i)  = RHT(i) / LOOP;
    RDCT(i) = RDCT(i) / LOOP;
    RKLT(i) = RKLT(i) / LOOP;
end;

% Energy Compaction 特性のグラフ表示
figure;
plot(nn,RDFT,'b'); hold on
plot(nn,RHT,'c');  hold on
plot(nn,RDCT,'g'); hold on
plot(nn,RKLT,'r');
axis([0, 64, 0.8, 1]);
legend('DFT', 'HT', 'DCT', 'KLT');
TITLE = sprintf('Energy Compaction: p=%f, n=%d', p, n);
title(TITLE);

%% KLT における変換係数の無相関性の3次元表示
%figure;
%[T,TX] = corrmtx(XKLT,n-1); % 変換係数の自己相関行列の作成
%TX = (1.0/TX(1,1)) * TX;    % 対角要素を1に正規化 (不要)
%[X,Y]=meshgrid(nn,nn);
%subplot(2,1,1);
%mesh(X,Y,RX);
%TITLE = sprintf('Correration of Input Signals: p=%f, n=%d', p, n);
%title(TITLE);
%subplot(2,1,2);   
%mesh(X,Y,TX);
%TITLE = sprintf('Correration of KLT Coefficients');
%title(TITLE);


