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

% 相関係数
p = 0.95;

% 試行回数
LOOP = 20;

% ブロック分割サイズ
size = 64;
ss = 1:size;

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

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

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

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


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

    % size x size KLT 行列の作成
    [R,RX] = corrmtx(x,size-1); % 入力信号の自己相関行列の作成
    RX = (1.0/RX(1,1)) * RX;    % 対角要素を1に正規化 (不要)
    [L,lambda] = eig(RX);       % 固有値展開 (L: 固有ベクトル行列、lambda: 固有値ベクトル)
    L = fliplr(L);              % 固有値の大きな順番に固有ベクトルを並べ替え

    % 入力信号をブロックに分割し、直交変換実行、および Energy Compaction 特性評価
    for j=1:n/size
        for k=1:size
            y(k) = x((j-1)*size + k);
        end;
        
        % 各種直交変換
        % 入力信号を直交変換し、変換係数ベクトルに格納
        %% DFT
        %% 変換係数の複素共役の関係を考慮して、実数、虚数の順番に並び替え
        xdft = (1.0/sqrt(size))*fft(y);
        XDFT(1) = xdft(1);
        XDFT(size) = xdft(size/2+1);
        for i=2:size/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(size)) * (y*hadamard(size));
        %% DCT
        XDCT = dct(y);
        %% KLT
        XKLT = y*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:size
            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:size
            ADFT(i) = ADFT(i) / ADFT(size);
            AHT(i)  = AHT(i) / AHT(size);
            ADCT(i) = ADCT(i) / ADCT(size);
            AKLT(i) = AKLT(i) / AKLT(size);
        end;
    
        % Energy Compaction 特性の累積評価結果の格納
        for i=1:size
            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;
end;

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

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

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


