clear all
close all
clc


signals = loadAllSignals('Supelec_2012_SIR_Spectral_Analysis_EA_v001.mat');
Fe = 1000000.0;
signal = signals(:,1);
t = (0: 1: length(signal)-1)/Fe;

%filtering smooth
nb_coeff = 10;
coeff = ones(1, nb_coeff)/nb_coeff;
signalm = filter(coeff, 1, signal);

%envelope time
signalenv = filter(coeff, 1, abs(hilbert(sgolayfilt(signal, 1, 3))));
plot_param = {'Color', [0.6 0.1 0.2],'Linewidth',1}; 

% savitzky-golay
signalsg = sgolayfilt(signal, 1, 3);

%fft
t0pad=2^12;%0-padding
f=((1: t0pad)-1)*Fe/t0pad-Fe/2;%echelle des frequences centrée en 0
TFx0pad=filter(coeff, 1, fft(signal, t0pad));%fft avec 0-padding
ModulTFx0pad=abs(TFx0pad);%module
CenterModulTFx0pad=fftshift(ModulTFx0pad);%shift zero-frequency component to centrer of spectrum
signalF=20*log10(CenterModulTFx0pad);%log
%fft smooth
TFx0pad=filter(coeff, 1, fft(signalm, t0pad));%fft avec 0-padding
ModulTFx0pad=abs(TFx0pad);%module
CenterModulTFx0pad=fftshift(ModulTFx0pad);%shift zero-frequency component to centrer of spectrum
signalmF=20*log10(CenterModulTFx0pad);%log
%fft hilbert
TFx0pad=filter(coeff, 1, fft(signalenv, t0pad));%fft avec 0-padding
ModulTFx0pad=abs(TFx0pad);%module
CenterModulTFx0pad=fftshift(ModulTFx0pad);%shift zero-frequency component to centrer of spectrum
signalenvF=20*log10(CenterModulTFx0pad);%log
%savitzky-golay
TFx0pad=filter(coeff, 1, fft(signalsg, t0pad));%fft avec 0-padding
ModulTFx0pad=abs(TFx0pad);%module
CenterModulTFx0pad=fftshift(ModulTFx0pad);%shift zero-frequency component to centrer of spectrum
signalsgF=20*log10(CenterModulTFx0pad);%log

%plot
figure
subplot(221)%%%%%%%%%%%%%%%%%%%%%%%
plot(t, signal)
xlabel('Time (s)')
ylabel('Amplitude (a.u.)')
hold on
plot(t, signalm)
xlim([0 length(signal)/Fe])
title('smooth filtering')
legend('signal', 'avg signal')
subplot(223)%%%%%%%%%%%%%%%%%%%%%%
plot(t, signal)
xlabel('Time (s)')
ylabel('Amplitude (a.u.)')
hold on
plot(t,[-1,1].*signalenv,plot_param{:})
plot(t, signalsg)
xlim([0 length(signal)/Fe])
title('Hilbert Envelope')
legend('signal', 'Hilbert Envelope', 'savitzky-golay')
subplot(222)%%%%%%%%%%%%%%%%%%%%%%
plot(f, signalF)
xlabel('f(Hz)');
ylabel('Module de la TF en dB');
grid on
hold on
plot(f, signalmF)
legend('signal', 'avg signal')
title('FFT')
xlim([0 Fe/2])
subplot(224)%%%%%%%%%%%%%%%%%%%%%%
plot(f, signalF)
xlabel('f(Hz)');
ylabel('Module de la TF en dB');
grid on
hold on
plot(f, signalenvF)
plot(f, signalsgF)
legend('signal', 'signal Hilbert', 'savitzky-golay')
title('FFT')
xlim([0 Fe/2])

figure
ax1 = subplot(3,1,1);
plot(t,signal)
ax2 = subplot(3,1,2);
pspectrum(signal,Fe,'spectrogram','OverlapPercent',99, 'Leakage',1,'MinThreshold',-150)
colorbar(ax2,'off')
ax3 = subplot(3,1,3);
pspectrum(signal,Fe,'spectrogram','OverlapPercent',99, 'Leakage',1,'MinThreshold',-150,'TimeResolution', 10e-6)
colorbar(ax3,'off')
linkaxes([ax1,ax2,ax3],'x')

figure
[sp,fp,tp] = pspectrum(signal,Fe,'spectrogram','OverlapPercent',99, 'Leakage',1);
mesh(tp,fp,sp)
view(-15,60)
xlabel('Time (s)')
ylabel('Frequency (Hz)')