clear;
close all;

% returns sampling frequency in Hz and data
[y,Fs] = audioread('myson.wav');

% Fs = sampling frequency, 32000 for fluteircam.wav
lenW = 0.04*Fs; % window of lenW samples, i.e. 40 ms

df=0.9765625;      %%% la dsp est calcul\'ee tous les df Hz
ff=-Fs/2:df:Fs/2; % length 32769 for fluteircam.wav
ffsub = ff(1:11:end); % compression x11, length 2979

tt = (0:length(y)-1)/Fs;
ttsub = (lenW/2:lenW:length(y))/Fs;

dsps = zeros(length(ffsub), floor((length(y)-lenW+1)/lenW));

% dsp fft
T = 1/Fs; % Sampling period
L = lenW; % Length of signal
t = (0:L-1)*T; % Time vector
fftsp = zeros(L/2+1, floor((length(y)-lenW+1)/lenW));
fftspCpx = zeros(L/2+1, floor((length(y)-lenW+1)/lenW));
f = Fs*(0:(L/2))/L;

for i = 0:floor((length(y)-lenW+1)/lenW)
    
    % compute dsp AR
    [~, ~, ~, ~, mydsp] = mylevinsondurbin(y(lenW*i+1:lenW*(i+1))',200,Fs);
    dsps(:,i+1) = mydsp(1:11:end)'; % compression x11, length 2979
    
    % compute dsp fft
    myfft = fft(y(lenW*i+1:lenW*(i+1)));
    P2 = abs(myfft/L);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);
    fftsp(:,i+1) = P1;
end

% take only positive frequencies for dsp
ffsubp = ffsub(1,(length(ffsub)-1)/2+1:end);
dspsp = dsps((length(dsps)-1)/2+1:end,:);

% plot
figure()
plot(tt,y)
xlabel('temps (s)')
ylabel('amplitude (u.a.)')
title('signal fluteircam')

% figure()
% surf(ttsub,ffsub,dsps,'EdgeColor','None');
% xlabel('temps (s)')
% ylabel('fréquences (Hz)')
% zlabel('amplitudes (u.a.)')
% title('Full DSP AR signal fluteircam')

figure()
surf(ttsub,ffsubp,dspsp,'EdgeColor','None');
xlabel('temps (s)')
ylabel('fréquences (Hz)')
zlabel('amplitudes (u.a.)')
title('DSP AR signal fluteircam')

figure()
imagesc(ttsub,ffsubp,dspsp)
xlabel('temps (s)')
ylabel('fréquences (Hz)')
title('Amplitude DSP AR signal fluteircam')

figure()
surf(ttsub,f,fftsp,'EdgeColor','None');
xlabel('temps (s)')
ylabel('fréquences (Hz)')
zlabel('amplitudes (u.a.)')
title('DSP FFT signal fluteircam')

figure()
imagesc(ttsub,f,fftsp)
xlabel('temps (s)')
ylabel('fréquences (Hz)')
title('Amplitude DSP FFT signal fluteircam')


% take max amplitude frequency
[maxDspsp, maxIndDspsp] = max(dspsp);
maxFfsubp = ffsubp(maxIndDspsp);
figure
plot(ttsub,maxFfsubp)
xlabel('temps (s)')
ylabel('fréquences (Hz)')
title('Frequency max DSP AR signal fluteircam')

[maxFftsp, maxIndFftsp] = max(fftsp);
maxF = f(maxIndFftsp);
figure
plot(ttsub,maxF)
xlabel('temps (s)')
ylabel('fréquences (Hz)')
title('Frequency max DSP FFT signal fluteircam')