test_rjmcmc_sinu.m

%clear;

clc , close all;

echo off;


boites_dialogue = 1;

font_size       = 13;

font_name       = 'times';

nb_classe       = 80;

%%--------- Parameters -----------%%

N               =  64;             %Number of data samples

sig             =  1;


if(~boites_dialogue)

    %param_sinu     =  [-0 , -0 , -0 ; 0.1*2*pi , (0.2 + 0/N)*2*pi  , (0.3 + 0/N)*2*pi ];

    param_sinu          =  [-3  -3 ; 0.1*2*pi 0.3*2*pi ];

else

    question1            = {'N' , 'Amplitude (dB)'  , 'Frequencies [0 0.5]' , '\sigma' };

    if (~exist('def1'))

        def1                 = {'64' , '[-3 , -3]' , '[ 0.1 , 0.3]' , '1' };

    end

    Titreboite1          = 'Measurements Parameters';

    nbligne1             = 1;

    AddOpts1.Resize      = 'on';

    AddOpts1.WindowStyle = 'normal';

    AddOpts1.Interpreter = 'tex';

    reponse1             = inputdlg(question1 , Titreboite1 , nbligne1 , def1 , AddOpts1);


    N                    = str2num(char(reponse1(1)));

    param_sinu           = [str2num(char(reponse1(2))) ; str2num(char(reponse1(3)))*2*pi];

    sig                  = str2num(char(reponse1(4)));

    def1                 = reponse1;


end


%%---------- Generate the data ----------%%

y              = generate_y_sinu(N , param_sinu , sig);

z              = y + sig*randn(N , 1);

%%---------------------------------------%%

hyper.T        = 5000;

hyper.burnin   = hyper.T/2;

hyper.kmax     = 12;

hyper.c        = 0.5;

hyper.O        = 12;


hyper.upsilon  = 0.1;

hyper.gamma    = 0.1;

hyper.epsi1    = 10e-4;

hyper.epsi2    = 10e-4;

hyper.alpha    = 2;

hyper.beta     = 5;

hyper.lambda   = 0.1;

hyper.merge    = 0.1;

hyper.sRW      = (1/(5*N))^2;

hyper.walk     = 0.0;

if(boites_dialogue)

    question2            = {'T' , 'burnin'  , 'k_{max}' , 'c : c=0 <=> no RJ, c=0.25 <=> 4 RJ, c=0.5 <=> 2 RJ' , 'O' , '\upsilon' , '\gamma' , '\epsilon_1' , '\epsilon_2' , '\alpha' , '\beta' , '\lambda' , 'merge' , '\sigma_{RW}' , 'walk'  };

    if (~exist('def2'))

        def2                 = {'5000' , '5000/5' , '12' , '0.5' , '14' , '0.1' , '0.1' , '10e-4' , '10e-4' , '2' , '5', '0.3' , '0.1' , '(1/(5*N))^2' , '0' };
    end

    Titreboite2          = 'RJ-MCMC Parameters';

    nbligne2             = 1;

    AddOpts2.Resize      = 'on';

    AddOpts2.WindowStyle = 'normal';

    AddOpts2.Interpreter = 'tex';

    reponse2             = inputdlg(question2 , Titreboite2 , nbligne2 , def2 , AddOpts2);

    hyper.T              = str2num(char(reponse2(1)));

    hyper.burnin         = str2num(char(reponse2(2)));

    hyper.kmax           = str2num(char(reponse2(3)));

    hyper.c              = str2num(char(reponse2(4)));

    hyper.O              = str2num(char(reponse2(5)));

    hyper.upsilon        = str2num(char(reponse2(6)));

    hyper.gamma          = str2num(char(reponse2(7)));

    hyper.epsi1          = str2num(char(reponse2(8)));

    hyper.epsi2          = str2num(char(reponse2(9)));

    hyper.alpha          = str2num(char(reponse2(10)));

    hyper.beta           = str2num(char(reponse2(11)));

    hyper.lambda         = str2num(char(reponse2(12)));

    hyper.merge          = str2num(char(reponse2(13)));

    hyper.sRW            = eval(char(reponse2(14)));

    hyper.walk           = str2num(char(reponse2(15)));

    def2                 = reponse2;


end


%%--------- Run the algorithm ------------%%

t0            = clock;

[out , param] = rjmcmc_sinu(z , hyper);

t1            = etime(clock , t0);

disp( ['Run time  = ' num2str(t1) 's for T = ' int2str(hyper.T)] );


%%------------ Parsing output -----------------%%

ak_est_db     = 10*log10(out.ak_est.^2/(2*out.sig_est^2));

y_est         = generate_y_sinu(N , [ak_est_db' ; out.wk_est'] , out.sig_est);


if (~isempty(param_sinu))

    disp(['fk_true = ' num2str(sort(param_sinu(2 , :))/(2*pi)) ] );

end

disp(['fk_est = ' num2str(sort(out.wk_est')/(2*pi)) ]);


disp(['sigma_true = ' num2str(sig')]);


disp(['sigma_est = ' num2str(out.sig_est')]);

if (~isempty(param_sinu))

    disp(['ak_true = ' num2str(sort(param_sinu(1 , :))) ] );

end

disp(['ak_est = ' num2str(sort(ak_est_db')) ]);


%%===================== Figures ===================%%



figure(1)

plot(param.k)

title('\bf k' , 'fontname' , font_name , 'fontsize' , font_size);


figure(2)

edge              =   0 : hyper.kmax;

histo             =   histc(param.k , edge)/hyper.T;

bar(edge , histo);

title('p({\bfk}|{\bfy})' , 'fontname' , font_name , 'fontsize' , font_size)


figure(3)

[histo , bin] = hist(param.delta(2:end) , nb_classe);

bar(bin , histo/(hyper.T - 1));

title('p(\bf{\delta}|\bf{y})' , 'fontname' , font_name , 'fontsize' , font_size);


figure(4)

[histo , bin] = hist(param.sigma(2:end) , nb_classe);

bar(bin , histo/(hyper.T - 1));

title('p(\bf{\sigma}|\bf{y})', 'fontname' , font_name , 'fontsize' , font_size)



figure(5)


p1 = plot((1 : N)' , y , (1 : N)' , z  , 'r-+' , (1 : N)' , y_est , 'k-^');

h = legend(p1 , ' {\bf y}' , ' {\bf z}' , ' {\bf y_{est}}');


figure(6)

S                      = hyper.O*N;

w                      = (0 : 1/(S -1 ) : 0.5);

temp_fft_y             = (1/N)*abs(fft(y , S)).^2;

perio_y                = temp_fft_y(1 : S/2);

temp_fft_z             = (1/N)*abs(fft(z , S)).^2;

perio_z                = temp_fft_z(1 : S/2);


p1                     = plot(w , perio_y , w , perio_z , 'r' , 'linewidth' , 2);

axis tight

xlabel('F')


h                      = legend(p1 , 'periodogram {\bf y}' , 'periodogram {\bf z}');


figure(7)

if (out.k_est > 0)

    temp  = find(param.k == out.k_est);

    plot((param.wk(1:out.k_est , temp)/(2*pi))')

    axis([0 length(temp) 0 0.5])

    ylabel('F')

end