%% tests for spectro_tfe
% 
% This script creates a pairs of sinusoid functions to validate an
% absolute value and an argument of the complex gain estimation from
% function spectro_tfe.
%
%% 1) Length of sinusoids even
%
% Specify parameters of sinusoids and construct a time bases
%
fs = 100;            % a sampling frequency
t = 0:(1/fs):10;     % the first time base
f0 = 3;              % a sinusoid frequency 
fprintf('length of t: %5i\n', length(t)); 
%
%% g = 0.8, phi = - 3 * pi / 4
%
% Define a sinusoids with a different phases and different amplitudes 
%
phi = - 3 * pi / 4; 
x1 = sin(2 * pi * f0 * t);
x2 = 0.8 * sin(2 * pi * f0 * t + phi);
[g, phi] = spectro_tfe(x1, x2, f0, fs, t); 
fprintf('g: %5.2f, phi: %5.2f\n', g, phi); 
%
%% g = 0.2, phi = + 3 * pi / 4
%
% Define a sinusoids with a different phases and different amplitudes 
%
phi = + 3 * pi / 4; 
x1 = sin(2 * pi * f0 * t);
x2 = 0.2 * sin(2 * pi * f0 * t + phi);
[g, phi] = spectro_tfe(x1, x2, f0, fs, t); 
fprintf('g: %5.2f, phi: %5.2f\n', g, phi); 
%% phi varying from -pi to pi
%
% 
% 
n = 60; 
phi = zeros(n); 
x1 = sin(2 * pi * f0 * t);
G_est = zeros(n); 
Phi_est = zeros(n); 
for i = 1 : n, 
    phi(i) = - pi + 2 * pi * i / n; 
    x2 = sin(2 * pi * f0 * t + phi(i));
    [g_est, phi_est] = spectro_tfe(x1, x2, f0, fs, t); 
    G_est(i) = g_est; 
    Phi_est(i) = phi_est; 
end

figure
plot(phi, Phi_est, '.-');
xlabel('\phi (rad)'); 
ylabel('\phi_{est} (rad)'); 


%% 2) Length of sinusoids odd
%
% 
%
fs = 100;            % a sampling frequency
t = 0:(1/fs):10-1/fs;     % the first time base
f0 = 3;              % a sinusoid frequency 
fprintf('length of t: %5i\n', length(t)); 
%% g = 0.8, phi = - 3 * pi / 4
%
% Define a sinusoids with a different phases and different amplitudes 
%
phi = - 3 * pi / 4; 
x1 = sin(2 * pi * f0 * t);
x2 = 0.8 * sin(2 * pi * f0 * t + phi);
[g, phi] = spectro_tfe(x1, x2, f0, fs, t); 
fprintf('g: %5.2f, phi: %5.2f\n', g, phi); 
%
%% g = 0.2, phi = + 3 * pi / 4
%
% Define a sinusoids with a different phases and different amplitudes 
%
phi = + 3 * pi / 4; 
x1 = sin(2 * pi * f0 * t);
x2 = 0.2 * sin(2 * pi * f0 * t + phi);
[g, phi] = spectro_tfe(x1, x2, f0, fs, t); 
fprintf('g: %5.2f, phi: %5.2f\n', g, phi); 
%% phi varying from -pi to pi
%
% 
% 
n = 60; 
phi = zeros(n); 
x1 = sin(2 * pi * f0 * t);
G_est = zeros(n); 
Phi_est = zeros(n); 
for i = 1 : n, 
    phi(i) = - pi + 2 * pi * i / n; 
    x2 = sin(2 * pi * f0 * t + phi(i));
    [g_est, phi_est] = spectro_tfe(x1, x2, f0, fs, t); 
    G_est(i) = g_est; 
    Phi_est(i) = phi_est; 
end

figure
plot(phi, Phi_est, '.-');
xlabel('\phi (rad)'); 
ylabel('\phi_{est} (rad)'); 


