%% DEMO on generated data

% -------------------------------------------------------------------------
% Generate Data
close all;
clear all;

fs=2000;
f1=900;
f2=500;
f3=400;
N=fs*2;

a=sin(2*pi*f1*(0:(N-1))/fs);
b=sin(2*pi*f2*(0:(N-1))/fs);
c=sin(2*pi*f3*(0:(N-1))/fs);

x=a+b+c;
y=awgn(x,0.1,'measured'); %add white noise with SNR of 0.1

number_lines=6; %3 sinusoids -> 6 spectral lines
%soundsc(y,fs);
disp(['-----------------------------']);
disp(['Results for generated data : ']);



% -------------------------------------------------------------------------
% MUSIC, compute the estimated spectral line using the MUSIC algorithm
est_music=music(y,number_lines,fs);
disp(['MUSIC res : ' num2str(est_music)]);



% -------------------------------------------------------------------------
% ESPRIT, compute the estimated spectral line using the ESPRIT algorithm
est_esprit=esprit(y,fs,number_lines);
est_esprit = reshape(est_esprit,[1,length(est_esprit)]);
disp(['ESPRIT res : ' num2str(est_esprit)]);



% -------------------------------------------------------------------------
% SAMV, compute the estimated spectral line using the SAMV algorithm
y=y(1:500); %cut the length to 500 since SAMV takes a lot of time
y=reshape(y, [length(y),1]);
[p_vectors, sigma, R]=samv(y, fs); %call SAMV

% Plot the frequencies with their magnitude (we then select the frequencies
% with the highest magnitudes)
f = linspace(0,1,fs/2)*fs/2;
figure;
plot(f,p_vectors(1:length(p_vectors)/2));
title("SAMV");

% Find the best frequencies and resize them between -fs/2 and fs /2
idxs=findmax_window(p_vectors,number_lines,3);
idxs=reshape(idxs,[1,length(idxs)]);
for i=1:length(idxs)
   if(idxs(i) >= fs/2)
       idxs(i)=idxs(i)-fs;
   end
end
disp(['SAMV res : ' num2str(idxs)]);







%% DEMO on real data

% Import bass line
[y,fs]=audioread('clean.wav');
[y_n,fs_n]=audioread('noisy.wav');

% We don't gain much precision by taking the whole line (we just want the
% four notes to be in our selection)
y=y(1:280000);
y_n=y_n(1:280000);

y=reshape(y,[1,length(y)]);
y_n=reshape(y_n,[1,length(y_n)]);

number_lines=24;
%soundsc(y,fs)

disp(['-----------------------------']);
disp(['Results for real data : ']);



% -------------------------------------------------------------------------
% MUSIC, compute the estimated spectral line using the MUSIC algorithm
est_music=music(y,number_lines,fs);
disp(['MUSIC res : ' num2str(sort(est_music))]);

est_music_n=music(y_n,number_lines,fs);
disp(['MUSIC res noisy : ' num2str(sort(est_music_n))]);



% -------------------------------------------------------------------------
% ESPRIT, compute the estimated spectral line using the ESPRIT algorithm
est_esprit=esprit(y,fs,number_lines);
est_esprit=reshape(est_esprit,[1,length(est_esprit)]);
disp(['ESPRIT res : ' num2str(sort(est_esprit))]);

est_esprit_n=esprit(y_n,fs,number_lines);
est_esprit_n=reshape(est_esprit_n,[1,length(est_esprit_n)]);
disp(['ESPRIT res noisy : ' num2str(sort(est_esprit_n))]);
