% Romain Monchaux - 12/10/2015
% M2 DFE
% Signal acquisition and processing
% Exercise 3: 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='signaljet';
dirfile='../signaux/';
A=load([dirfile,file]);% chargement du fichier
Fs=10000;% fréquence d'échantillonage
time = A(:,1); % vecteur temps
signal=A(:,2); % signal de vitesse longitudinale
Nsamples=length(time);
duration=?;
%
disp(['Signal de durée ',num2str(duration),'s capté à Fs=',...
    num2str(Fs/1000),'kHz'])
%
% Tracer le signal temporel:
figure(2011)% le forçage
plot(time,signal,'k-','LineWidth',2)
set(gca,'FontSize',15)
title(['Signal de durée ',num2str(duration),'s capté à Fs=',...
    num2str(Fs/1000),'kHz'])
xlabel('t (s)')
ylabel('U_x (m/s)')
xlim([? ?])
% enregistrement de la figure:
saveas(gcf,['../figures/jet','.png'])
%
% réaliser quelques zooms pour trouver "à l'oeil" des fréquences dans le
% signal :
%%
figure(2012)% la figure zoomée
plot(time,signal,'k-','LineWidth',2)
set(gca,'FontSize',15)
title(['Signal de durée ',num2str(duration),'s capté à Fs=',...
    num2str(Fs/1000),'kHz'])
xlabel('t (s)')
ylabel('U_x (m/s)')
% paramètres du zoom:
x0=?;
dX=?;
xlim([x0 x0+dX])
% enregistrement de la figure:
saveas(gcf,['../figures/jet',num2str(1000*dX),'.png'])
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcul des grandeurs statistiques sur des intervalles d'une seconde
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
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_U = zeros(1,Nintervals);
std_U = zeros(1,Nintervals);
Nsamples_int = Tint*Fs; % nombre d'échantillons par intervalle
%
% boucle sur les intervalles
for ii = 1:Nintervals 
    mean_U(ii) = mean(signal((ii-1)*Nsamples_int+1:ii*Nsamples_int));
    std_U(ii) = std(signal((ii-1)*Nsamples_int+1:ii*Nsamples_int));
end
%
figure(110)
plot(t1s,mean_U,'ko-','LineWidth',2)
hold on
plot(t1s,std_U,'r^-','LineWidth',2)
set(gca,'FontSize',15)
xlabel(['t_{',num2str(Tint),'s} (s)'])
ylabel('statistique de vitesse (m/s)')
legend('moyenne','ecart-type','Location','East')
grid on
saveas(gcf,['../figures/jet_moy',num2str(1000*Tint),'.png'])
%
% return % à commenter quand la question 1 est finie
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 2: 
% calculer la densité spectrale de puissance du signal en utilisant pwelch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% on considere le signal complet (pas de moyenne):
Nfft = ?;%Nsamples
Nwindow = ?; % 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 = ? % résolution fréquentielle
%
% tracé de la densité spectrale de puissance:
figure(201)
hold on
box on
plot(freq,PSD,'k-','LineWidth',2)
set(gca,'FontSize',15)
xlabel('f (Hz)','FontSize',15)
ylabel('PSD (m^2s^{-2}/Hz)','FontSize',15)
grid on
set(gca,'Yscale','Log')
set(gca,'Xscale','Log')
saveas(gcf,['../figures/jet_freq_',num2str(Nfft),'.png'])
% 
%%
% Calculer le spectre moyen à une résolution fréquentielle de 0.05 Hz
df = ?; % resolution fréquentielle
Nfft = ?;
Noverlap=Nfft/2;
%
% On utilise la fonction pwelch de Matlab:
Nwindow = Nfft; % nombre de points de la fenêtre
%
[PSD,freq] = pwelch(signal,Nwindow,Noverlap,Nfft,'psd',Fs);

% tracé:
figure(301)
hold on
box on
set(gca,'FontSize',15)
plot(freq,PSD,'b-','LineWidth',2)
set(gca,'Yscale','Log')
set(gca,'Xscale','Log')
xlabel('f (Hz)')
ylabel('PSD (m^2s^{-2}/Hz)')
grid on
set(gca,'XMinorTick','on')
set(gca,'YMinorTick','on')
%
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 3: 
% Calculer le periodogramme moyen à une résolution fréquentielle de 0.5 Hz
df = ?; % resolution fréquentielle
Nfft = ?;
Noverlap=Nfft/2;
%
% On utilise la fonction pwelch de Matlab:
Nwindow = Nfft; % nombre de points de la fenêtre
%
[PSD,freq] = pwelch(signal,Nwindow,Noverlap,Nfft,'psd',Fs);
%
% tracé:
figure(301)
hold on
box on
plot(freq,PSD,'r-','LineWidth',2)

%%
% Calculer le periodogramme moyen à une résolution fréquentielle de 5Hz
df = ?; % frequency resolution
Nfft = ?;
Noverlap=Nfft/2;
%
% On utilise la fonction pwelch de Matlab:
Nwindow = Nfft; % nombre de points de la fenêtre
%
[PSD,freq] = pwelch(signal,Nwindow,Noverlap,Nfft,'psd',Fs);
%
% tracé:
figure(301)
hold on
box on
plot(freq,PSD,'-','LineWidth',2,'Color',[.1 .5 .1])
legend('\Delta f = 0.05Hz','\Delta f = 0.5Hz','\Delta f = 5Hz',1)
set(gca,'FontSize',15)
set(gca,'Yscale','Log')
%
xlim([.01 10000])
ylim(10.^[-12 2])
saveas(gcf,['../figures/jet_freq_VarDf.png'])
%
%% Autocorrélation:
Sig=?;
%
[Cy,Lags] = xcorr(Sig,'coeff');
figure
hold on
box on
set(gca,'FontSize',18)
xlabel('\Delta t')
ylabel('Autocorrélation normalisée')
plot(Lags/Fs,Cy,'-o')
plot(Lags/Fs,zeros(size(Lags)),'-.k')
xlim([-0.00001 .2])
title('jet turbulent')
saveas(gcf,['../figures/Jet_COR.png'])
%% Tracé des densités de probabilité:
% choix du signal
Sig=?;
NbPt=length(Sig);
NbBin=?;
%
figure(1002)
hold on
box on
set(gca,'FontSize',18)
%
BG=linspace(-5,5,2^20);
Gauss=1/(sqrt(2*pi))*exp(-(BG.^2)/2);
[Bin,PDF]=Calc_PDF(Sig(1:NbPt),NbBin);
subplot(2,1,1)
set(gca,'FontSize',18)
hold on
plot(BG,Gauss,'-.k')
plot(Bin,PDF,'o','LineWidth',1.5)
xlabel('U (m/s)')
ylabel('PDF')
xlim([-6 6])
ylim(10.^[-6 0])
% set(gca,'Yscale','lin')
set(gca,'Yscale','log')
subplot(2,1,2)
set(gca,'FontSize',18)
hold on
plot(BG,Gauss,'-.k')
plot(Bin,PDF,'o','LineWidth',1.5)
xlabel('U (m/s)')
ylabel('PDF')
% set(gca,'Yscale','lin')
xlim([-2 2])
ylim([0 .5])
%
titre=(['Jet turbulent']);
text(-1.2 ,1.38,titre,'FontSize',18)
saveas(gcf,['../figures/jet_PDF.png'])

