% Romain Monchaux - 12/10/2015
% M2 DFE
% Signal acquisition and processing
% Exercise 2: Analyse de signaux réels
%
clc
clear all;
close all;
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 1: 
% Charger et tracer le signal original et tester sa stationnarité
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% charger le signal:
file='gong_force.wav';
dirfile='../signaux\';
[signal,Fs,nbits] = wavread([dirfile,file]);
% signal: le signal, Fs: fréquence d'échantillonage, nbits: le nombre de bits
% du codage 
%
Ts=1/Fs; % pas de temps
Nsamples = length(signal); % nombre d'échantillons
time = (0:Nsamples-1)*Ts; % vecteur temps
duration = time(Nsamples); % durée du signal (s)
%
% Séparation forçage/réponse:
Fsig=signal(:,1);
Psig=signal(:,2);
%
disp(['Signal de durée ',num2str(duration),'s capté à Fs=',...
    num2str(Fs),'Hz sur ',num2str(nbits),' bits '])
%
% Tracer les signaux temporels:
figure(100)% le forçage
subplot(2,1,1)
plot(time,Fsig,'k-','LineWidth',2)
set(gca,'FontSize',15)
title(['Signal de durée ',num2str(duration),'s capté à Fs=',...
    num2str(Fs/1000),'kHz sur ',num2str(nbits),' bits '])
xlim([0 duration])
xlabel('t (s)')
ylabel('Af (UA)')
%
subplot(2,1,2)% la réponse en pression
plot(time,Psig,'k-','LineWidth',2)
set(gca,'FontSize',15)
title(['Signal de durée ',num2str(duration),'s capté à Fs=',...
    num2str(Fs/1000),'kHz sur ',num2str(nbits),' bits '])
xlabel('t (s)')
ylabel('p (Pa)')
xlim([0 duration])
saveas(gcf,'../figures/gong.png')
%
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcul des grandeurs statistiques sur des intervalles de longeur variable
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
Tint=?; % longueur des intervalles (s)
Nintervals=round(duration/Tint); % nombre d'intervalles
t1s =(0:Nintervals-1)*Tint+Tint/2; % vecteur temps
%
% initialisation des moyennes et écarts-type de pression et de forçage:
mean_p1s = zeros(1,Nintervals);
std_p1s = zeros(1,Nintervals);
mean_f1s = zeros(1,Nintervals);
std_f1s = zeros(1,Nintervals);
Nsamples_int = Tint*Fs; % nombre d'échantillons par intervalle
%
% boucle sur les intervalles
for ii = 1:Nintervals 
    mean_p1s(ii) = mean(Psig((ii-1)*Nsamples_int+1:ii*Nsamples_int));
    std_p1s(ii) = std(Psig((ii-1)*Nsamples_int+1:ii*Nsamples_int));
    mean_f1s(ii) = mean(Fsig((ii-1)*Nsamples_int+1:ii*Nsamples_int));
    std_f1s(ii) = std(Fsig((ii-1)*Nsamples_int+1:ii*Nsamples_int));
end
%
figure(110)
close(110)
figure(110)
subplot(2,1,1)
plot(t1s,mean_p1s,'ko-','LineWidth',2)
hold on
plot(t1s,std_p1s,'r^-','LineWidth',2)
set(gca,'FontSize',15)
xlabel('t_{1s} (s)')
ylabel('pression (Pa)')
% legend('moyenne','ecart-type',2)
grid on
%
subplot(2,1,2)
plot(t1s,mean_f1s,'ko-','LineWidth',2)
hold on
plot(t1s,std_f1s,'r^-','LineWidth',2)
set(gca,'FontSize',15)
xlabel('t_{1s} (s)')
ylabel('forçage (UA)')
% legend('moyenne','ecart-type','Location','SouthWest','FontSize',9)
grid on
saveas(gcf,['../figures/gong_moy',num2str(1000*Tint),'.png'])
%
return % à commenter quand la question 1 est finie
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 2: 
% calculer la densité spectrale de puissance du signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% on considere le signal complet (pas de moyenne):
Nfft = ;%Nsamples
Nwindow = Nfft; % nombre de points de la fenêtre
Noverlap = Nfft/2; % longueur de recouvrement
% Calcul du spectre
[PSD,freq] = pwelch(signal(1:Nwindow),Nwindow,Noverlap,Nfft,'psd',Fs);
df = Fs/Nfft % résolution fréquentielle
%
% tracé de la densité spectrale de puissance:
figure(200)
hold on
plot(freq,PSD,'k-','LineWidth',2)
set(gca,'FontSize',15)
xlabel('f (Hz)')
ylabel('PSD (Pa)^2/Hz)')
grid on
% xlim([0 1000]) % zoom en fréquence 10kHz (ajustable)
set(gca,'Yscale','Log') % passage en échelle log
saveas(gcf,'gong_spec_log.png')
% set(gca,'Xscale','Log') % passage en échelle log
set(gca,'Yscale','Lin') % passage en échelle linéaire
saveas(gcf,'../figures/gong_spec_lin.png')
% 
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 3: 
% calculer le spectrogramme en utilisant une FFT à temps court
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spectrogramme de la réponse en pression:
clc
df = ?; % resolution fréquentielle
Nfft = ?; % nombre de points de la FFT
Nwindow = Nfft; % nombre de points de la fenêtre
Noverlap = Nfft/2;
K = ?; % nombre de segments de la FFT
% calculer la DSP dans chaque segment et les temps correspondants
Specgram = zeros(Nfft/2+1,K);
tspecgram = zeros(1,K);
for ii = 1:K % boucle sur les segments
   [Specgram(:,ii),freq_specgram] = pwelch(Psig((ii-1)*Nwindow+1:ii*Nwindow),Nwindow,Noverlap,Nfft,'psd',Fs);
%[Specgram(:,ii),freq_specgram] = periodogram(Psig((ii-1)*Lwindow+1:ii*Lwindow),...
%        window,Nfft,Fs);
    tspecgram(ii) = (ii-0.5)*Nwindow*Ts; % temps correspondant au milieu de chaque segment
end
dt_specgram = Nwindow*Ts; % resolution temporelle
%
figure(31)
pcolor(tspecgram,freq_specgram(1:end)/1000,log(Specgram(1:end,:)))
shading interp
set(gca,'FontSize',15)
xlabel('t (s)')
ylabel('f (kHz)')
colorbar
ylim([0 4]) % zoom en fréquence
saveas(gcf,'../figures/gong_spectrogramme.png')
%
%%
% spectrogramme du forçage:
dt_specgram = 0.05; % temporal resolution
r = 1-dt_specgram*Fs/Nfft; % overlap (r between 0 and 1)
noverlap = r*Lwindow; %  % number of overlap points
window=hanning(Nwindow);
[STFT3,freq3,t3,PSD3] = spectrogram(Fsig,window,noverlap,Nfft,Fs);
%
figure(340)
imagesc(t3,freq3/1000,log(PSD3))
axis xy
colorbar
set(gca,'FontSize',15)
xlabel('t (s)')
ylabel('f (kHz)')
ylim([0 4])
saveas(gcf,'../figure/forcage_spectrogramme.png')
%
%%
figure(32)
subplot(3,1,1)
plot(freq_specgram(1:end)/1000,Specgram(:,50),'k')
set(gca,'Yscale','Log') % passage en échelle log
set(gca,'FontSize',15)
xlim([0 15])
xlabel('f (kHz)')
text(12,10^(-4),['t=',num2str(tspecgram(50))])
%
subplot(3,1,2)
plot(freq_specgram(1:end)/1000,Specgram(:,150),'b')
set(gca,'Yscale','Log') % passage en échelle log
set(gca,'FontSize',15)
xlabel('f (kHz)')
ylabel('PSD (hPa)^2/Hz)')
xlim([0 15])
text(12,10^(-4),['t=',num2str(tspecgram(150))])
%
subplot(3,1,3)
plot(freq_specgram(1:end)/1000,Specgram(:,225),'r')
set(gca,'Yscale','Log') % passage en échelle log
set(gca,'FontSize',15)
xlabel('f (kHz)')
set(gcf,'Position',[ 680   100   560   878])
xlim([0 15])
text(12,10^(-4),['t=',num2str(tspecgram(225))])
saveas(gcf,'../figures/gong_spectrogramme_discret.png')
%
%%


