%spectrum_compressed.m

%different mat files can be used by just removing the comment below. Now, it is active for lowpass.mat



%Computation of Phi (The matrix of N_block and M_ruler) to compress the

%exisiting data.

%load('lowpasssignal.mat');
%load('multibandsignal.mat');
%load('highpasssignal.mat');
load('bandpasssignal.mat');

mn = dsp.Mean;
y_sig_new = y_sig;
y_sig_mean = mn(y_sig_new);
y_sig_updated = y_sig_new - y_sig_mean;
%y_sig_updated = y_sig_new;


Phi = zeros(M_ruler, N_block);
for i = 1:M_ruler
    Phi(i, ruler_index(i)) = 1;
end;

fs = 2* pi; %sampling of 2*pi becuase signal is complex valued

%computation of Ry - Autocorrelation of y in a M * M matrix and putting
%every element such that it becomes a M^2 * 1 matrix.

Ry_final = zeros(M_ruler^2, 1);


for i=1:M_ruler

    Ry_original_columns = zeros(M_ruler, 1);


    for l=1:M_ruler
        Ry_original_columns(l) = Ry_eval(l - 1, i - 1, y_sig_updated);
    end

    Ry_final((i - 1)*M_ruler + 1:i * M_ruler) = Ry_original_columns;
end




%Computation of T matrix to have Rx in a column matrix of size (N^2 * (2N -1))

T_final = zeros(N_block^2, 2*N_block - 1);

ID_mat = eye(2*N_block - 1);

for j=1:N_block^2
    T_final(j, :) = ID_mat(mod(j - 1 + (N_block - 2) * floor((j - 1) / N_block),(2 * N_block - 1))+1, :);
end


%Computation of RC matrix

Kron_Phi = kron(Phi, Phi);
R_c = Kron_Phi * T_final;

% Prediction of Rx, Rx_cap by Least Square:

Rx_cap = (R_c' * R_c)\R_c' * Ry_final;

%Retranslation of Rx from Rx_cap

Rx_cap_column = T_final * Rx_cap;

%R_x original form

Rx_cap_N_N = zeros(N_block);

for i = 1:N_block
   Rx_cap_N_N(:,i) = Rx_cap_column((i - 1) * N_block + 1: i * N_block);
end

%Direct Dft PSD

Px_cap = fft(Rx_cap);

freq_px = 0:fs/length(Px_cap):fs;

freq_px_updated = freq_px(1:length(freq_px)-1);

hold on
plot(freq_px_updated/pi, abs(Px_cap), 'LineWidth',2);
grid on
title('Power Spec Est of X -  Band Pass Signal')
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Power')


%PSD using Yule_Walker AR model.

%Solutions Yule_walker equations for AR using levinson algorithm:
%Trying number of coefficients as the lenghth from 1 to N_block - 1.

ang_freq = 0:fs/length(Rx_cap):fs;
ang_freq = ang_freq(1:length(ang_freq)-1);


% Yule_Walker equations to solve for A and B matrix:

Rx_Y = Rx_cap_N_N(2:end, 2:end);
rx_Y = Rx_cap_N_N(2:end, 1);

% Processes for AR 1 to AR max i.e N_block - 1

P = N_block - 1;

err = zeros(P, 1);

for l = 1:P
    Rx_l = Rx_Y(1:l, 1:l);

    rx_l = rx_Y(1:l);

    a_m = -Rx_l\rx_l; % A matrix with l co-efficients

    epsilon = Rx_l(1,1);
    denominator = 0;

    for m = 1:l
       epsilon = epsilon + a_m(m) * conj(rx_l(m));
       denominator = denominator + a_m(m) * exp(-1i * m* ang_freq);
    end
    b_0 = sqrt(epsilon);

    err(l) = length(y_sig) * log10(epsilon) + 2 * l; % To find out the best parametric solution
    Px_l = epsilon ./ (1+ denominator); % Power spectrum complex

    %plot(ang_freq/pi, abs(Px_l)) %This can be commented out to analyze all the plots from order 1 to 83

end

% For AR 12 process for example

l = 12;

Rx_l = Rx_Y(1:l, 1:l);

rx_l = rx_Y(1:l);

a_m = -Rx_l\rx_l; % a matrix with l co-efficients

epsilon = Rx_l(1,1);

denominator = 0;

for m = 1:l
       epsilon = epsilon + a_m(m) * conj(rx_l(m));
       denominator = denominator + a_m(m) * exp(-1i * m* ang_freq);
end
b_0 = sqrt(epsilon);

Px_l = epsilon ./ abs(1+ denominator); % Power spectrum complex


hold on
plot(ang_freq/pi, abs(Px_l), 'LineWidth',2)


grid on
%title('PSD | AR model')
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Power Spectrum | AR process')

% The minimum of variable err gives the best estimate with AR processe, the
% suitable number of parameters required in the AR process.

[v, index] = min(err); %value and index of err to find the best parameter value for the AR process.
l = index; % optimal index

Rx_l = Rx_Y(1:l, 1:l);

rx_l = rx_Y(1:l);

a_m = -Rx_l\rx_l; % a matrix with l co-efficients

epsilon = Rx_l(1,1);

denominator = 0;

for m = 1:l
       epsilon = epsilon + a_m(m) * conj(rx_l(m));
       denominator = denominator + a_m(m) * exp(-1i * m* ang_freq);
end
b_0 = sqrt(epsilon);

Px_l = (b_0)^2 ./ abs(1+ denominator); % Power spectrum complex
hold on
plot(ang_freq/pi, abs(Px_l), 'LineWidth',2)

grid on
%title('PSD | AR model')
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Power Spectrum | AR Process')

legend('Using DFT of Estimated Rx', 'AR Process with l = 12', 'AR with l = 52')