# APPENDIX

Here is our matlab code in order to study data from NSMB (FFT, Welch method...)

%---------------------------------------------------------------------------------------------------------------------------------------------------------

%----------------------------------------------------------------FFT Approach--------------------------------------------------------------------

%---------------------------------------------------------------------------------------------------------------------------------------------------------

% FFT approach for singal

% Author: Jingyan DUO-Christophe HEUDES, ENSEEIHT

% Creation: mars 2015

% Last modification: J. DUO, C. HEUDES

clear all

close all

%PARAMETERS

D=0.05715; %diameter

Uinf=44;   %velocity upstream

%fft on pressure data which is the 11th column

monitor=uigetfile

monitor=importdata(monitor);

n=length(monitor(:,11));

%---------------------------------------------------------------periodogram PSD estimate---------------------------------------------------

[psdper,Stt]=periodogram(monitor(:,11)-mean(monitor(:,11)),hann(n),2^(nextpow2(n)+6),1./(monitor(2,3)-monitor(1,3)));

frequence=Stt*Uinf/D;

[y x]=max(psdper)

St_mood1=Stt(x);

f_mood1=St_mood1*Uinf/D;

figure(1)

loglog(Stt,psdper,'k');

hold on;

line([St_mood1 St_mood1],[10e-20 1]);

hold off;

axis([1E-4 1E2 1e-20 1E3]);

title('periodogram PSD estimate as function of St on the monitor point 1 of the 1st cyliner','fontsize',14)

xlabel('Strouhal number','fontsize',12)

ylabel('PSD (Pa?/s)','fontsize',12)

legend('PSD as function of St with periodogram estimate','Strouhal number');

figure(2)

loglog(frequence,psdper,'k');

hold on;

line([f_mood1 f_mood1],[10e-20 1]);

line([1800 1800],[10e-20 1]);

hold off;

axis([1E-2 1E5 1e-20 1E3]);

legend('PSD as function of frequency','Kármán vortex frequency','shear layer vortex frequency')

title('periodogram PSD estimate  as function of frequency on the monitor point 1 of the 1st cyliner','fontsize',14)

xlabel('frequency (Hz)','fontsize',12)

ylabel('PSD (Pa?/s)','fontsize',12)

%----------------------------------------------------------------Welch PSD estimate----------------------------------------------------------

[psdwelch1,St] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(2000),1600,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));

[psdwelch2,St] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(1500),1200,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));

[psdwelch3,St] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(1000),800,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));

% %length of 2000 samples with 1700 overlapped samples etc. since here the fs=118.34

frequence=St*Uinf/D;

[y x]=max(psdwelch1)

St_mood1=St(x);

f_mood1=St_mood1*Uinf/D;

figure(3)

loglog(St,psdwelch1,'g');

hold on;

loglog(St,psdwelch2,'k');

loglog(St,psdwelch3,'r');

line([St_mood1 St_mood1],[10e-20 1]);

hold off;

xlabel('Strouhal')

ylabel('PSD on the monitor ponit 1')

legend('PSD as function of St with nfft=2000','PSD as function of St with nfft=1500','PSD as function of St with nfft=1000','Strouhal number St=0.23');

title('Welch PSD estimate as function of St on the monitor point 1 of the 1st cyliner','fontsize',14');

figure(4)

loglog(frequence,psdwelch1,'g');

hold on;

loglog(frequence,psdwelch2,'k');

loglog(frequence,psdwelch3,'r');

line([f_mood1 f_mood1],[10e-20 1]);

line([1800 1800],[10e-20 1]);

hold off;

legend('PSD as function of frequency with nfft=2000','PSD as function of frequency with nfft=1500','PSD as function of frequency with nfft=1000','Kármán vortex frequency f=174.13 Hz','shear layer vortex frequency f=1800Hz')

title('Welch PSD estimate as function of St on the monitor point 1 of the 1st cyliner','fontsize',14)

xlabel('frequency (Hz)','fontsize',12)

ylabel('PSD (Pa?/s)','fontsize',12)

%--------------------------------------------------------------------comparaision-----------------------------------------------------------------

%point number 1 of 1st and 2nd cylinder

display('Please select the file of the 1st ponit of 1st cylinder')

monitor1=uigetfile

monitor1=importdata(monitor1);

display('Please select the file of the 1st point of 2nd cylinder')

monitor2=uigetfile

monitor2=importdata(monitor2);

display('Please select the file of the 8th point of the 1st cylinder')

monitor3=uigetfile

monitor3=importdata(monitor3);

display('Please select the file of the 8th point of the 2nd cylinder')

monitor4=uigetfile

monitor4=importdata(monitor4);

figure(5)

plot(monitor1(:,5));

hold on;

plot(monitor2(:,5));

plot(monitor3(:,5));

plot(monitor4(:,5));

hold off;

legend('1st point of 1st cylinder','1st point of 2nd cylinder','8th point of 1st cylinder','8th point of 2nd cylinder');

title('horizontal speed of the flow u/Uinf');

xlabel('step');

ylabel('horizontal speed/Uinf');

%----------------------------------------------------------combinaision des deux methodes----------------------------------------------

clear all

close all

%PARAMETERS

D=0.05715; %diameter

Uinf=44;   %velocity upstream

%importer des donness

%fft on pressure data which is the 11th column

monitor=uigetfile

monitor=importdata(monitor);

n=length(monitor(:,11));

%pour bas frequence(f<7699), on utilise le methode de periodogram estimate

%pour haut frequence(f>7699), on utilise le methode de welch PSD estimate

N=10359;%pour amont, point 40

n1=88606;%logeur de bas frequence

n2=44304;%logeur de haut frequence

%diagramme de PSD en fonction de St

[psdper,St1]=periodogram(monitor(:,11)-mean(monitor(:,11)),hann(n),2^(nextpow2(n)+6),1./(monitor(2,3)-monitor(1,3)));

[y x]=max(psdper)

St_mood1=St1(x);

[psdwelch,St2] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(1000),800,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));

figure(6)

loglog(St1,psdper,'g.');

hold on;

loglog(St2,psdwelch,'r.');

psd1=psdper(1:n1,:);

psd2=psdwelch(n2:length(St2),:);

loglog(St1(1:n1,:),psd1,'k');

loglog(St2(n2:length(St2),:),psd2,'k');

line([10 10],[10e-20 1e4]);

line([St_mood1 St_mood1],[10e-20 1e4]);

line([2.34 2.34],[10e-20 1e4]);

hold off;

title('PSD estimate as function of St on the monitor point 1 of the 1st cyliner, with St_{cri}=10');

legend('Periodogram PSD estimate','Welch PSD estimate','Periodogram PSD estimate for low frquency','Welch PSD estimate for high frequency');

xlabel('Strouhal nomber St');

ylabel('PSD estimate (Pa^2/s)');

%diagramme de PSD en fonction de f

figure(7)

loglog(St1*Uinf/D,psdper,'g.');

hold on;

loglog(St2*Uinf/D,psdwelch,'r.');

psd1=psdper(1:n1,:);

psd2=psdwelch(n2:length(St2),:);

loglog(St1(1:n1,:)*Uinf/D,psd1,'k');

loglog(St2(n2:length(St2))*Uinf/D,psd2,'k');

line([7699 7699],[1e-20 1]);

line([1800 1800],[1e-20 1]);

line([St_mood1*Uinf/D St_mood1*Uinf/D],[1E-20 1]);

hold off;

title('PSD estimate as function of frequency on the monitor point 1 of the 1st cyliner, with f_{cri}=7699 Hz');

legend('Periodogram PSD estimate','Welch PSD estimate','Periodogram PSD estimate for low frquency','Welch PSD estimate for high frequency');

xlabel('frequency (Hz)');

ylabel('PSD estimate (Pa^2/s)');

%---------------------------------------------------------------------------------------------------------------------------------------------------------

%-------------------------------------------------------------------Wavelet transform-----------------------------------------------------------

%---------------------------------------------------------------------------------------------------------------------------------------------------------

% Compute CWT - Simple version

% Author: Damien SZUBERT - IMFT

% Creation: June 2013

% Last modification: D. SZUBERT  23II2015

%

% input :

% - tt: time vector (constant time step),

% - sig: signal (don't forget to remove mean!),

% - wname: wavelet name,

% - scales: scales of CWT (use scal2frq)

% - plotflag: plot type ('surf','imagesc','contourf')

%

% output :

% - cwt_data: coefficient of wavelet tranform,

% - pseudofreq: frequencies of CWT

%

% RMK: scales corresponds to periods in terms of samples (more or less)

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [cwt_data,pseudofreq] = monitorCWT_simple(tt,sig,wname,scales,plotflag)

pseudofreq = scal2frq(scales,wname,tt(2)-tt(1)); % Calculate frequency values from scales

cwt_data = cwt(sig,scales,wname); % Compute Continuous Wavelet Transform

%save('cwt_data.mat','cwt_data','tt','pseudofreq','-v7.3')

cwt_data_min = min(min(abs(cwt_data)));

cwt_data_max = max(max(abs(cwt_data)));

disp(['Min / max : ',num2str(cwt_data_min),' / ',num2str(cwt_data_max)])

%% Plot

if  strcmp(plotflag,'surf')

surf(tt,pseudofreq,abs(cwt_data),'edgecolor','none');view(2)

elseif  strcmp(plotflag,'imagesc')

imagesc(tt,pseudofreq,abs(cwt_data))

warning('Pay attention to the frequency scale! (log/linear axis/spacing)')

elseif strcmp(plotflag,'contourf')

contourf(tt,pseudofreq,abs(cwt_data),'edgecolor','none')

disp('Note: adapt number of contours if needed.')

end

%new_xtick_labels = cellfun(@(x) sprintf('%5.4f',x), num2cell(tt(get(gca,'XTick'))), 'uniformoutput', false);

%set(gca,'XTickLabel',new_xtick_labels)

%new_ytick_labels = cellfun(@(y) sprintf('%5.0f',y), num2cell(pseudofreq(get(gca,'YTick'))), 'uniformoutput', false);

%set(gca,'YTickLabel',new_ytick_labels,'YDir','reverse')

xlabel('Time (s)')

ylabel('Strouhal number')

xlim([tt(1) tt(end)])

ylim([pseudofreq(end) pseudofreq(1)])

%set(gca,'clim',[0 0.4*cwt_data_max])

end

%ondelette nous donne un diagramme de frequence en fonction de temps

%on a besoin le script 'monitorCWT_simple.m'

% % % % %n'a pas besoin de commander ''plot''

% % % % %faut calculer nbr de point, par example, dt=0.0084, St=0.23,

% % % % %period=1/St=4.34, nbr_point=period/dt=514, donc on prend 450:600

% % % % %'cmor1.5-1' c'est un caractere de ondelette, change pas ici!

[cwt_data pseudofreq]=monitorCWT_simple(monitor(:,3),burst,'cmor1.5-1',20:70,'contourf');

title('frequency modulation by wavelets');

colorbar;

%axis([0 80 1E-20 3]);

pseudofreq(end)

%---------------------------------------------------------------------------------------------------------------------------------------------------------

%----------------------------------------------AR-Autoregressive power spectral density estimate--------------------------------

%---------------------------------------------------------------------------------------------------------------------------------------------------------

clear all

clear y

close all

monitor=uigetfile

monitor=importdata(monitor);

[burst]=burst_process(monitor(:,11)-mean(monitor(:,11)),2.2,monitor(:,3));

sig=burst;

tt=monitor(:,3);

save_trace = 0;

font_size_label = 20;

font_name = 'Helvetica';

pAR = 200;                % AR model order

% pARvect = 5:10:100;

% col = hsv(length(pARvect));

%sigtemp = figure('name','Signal temporel');

nfft = 1000;            % FFT length calculation

fs = 1/(tt(2)-tt(1));   % Sampling frequency

dec_per = 4; % decalage periode

%per_ind = 12683;

per_ind =9172;

%deb = (12683*dec_per)+8000;

deb = 1;

sigtemp = figure(dec_per);

seg_length = 800;

nbr_seg = floor(per_ind/seg_length);

seg_length = floor(per_ind/nbr_seg); % longueur de segment ajustÃ©e pour prendre en compte le maximum de point

if seg_length < pAR

warning('seg_length < pAR')

end

col2 = hsv(nbr_seg);

%y(1,:) = sig_burst(deb:deb+per_ind,1)';

y(1,:) = sig(deb:deb+per_ind,1)';

ttloc = tt(deb:deb+per_ind)-tt(deb);

figure(1)

for k=1:nbr_seg

yseg(k,:) = y(1,(k-1)*seg_length+1:k*seg_length+1);

[Pxx(k,:),f(k,:)] = pyulear(yseg(k,:)-mean(yseg(k,:)),pAR,nfft,fs);

[PxxMax,indPxxMax] = max(Pxx(k,:));

fPxxMax(k) = f(k,indPxxMax);

loglog(f(k,:),Pxx(k,:))

pause

end

%% FFT

%  Y = fft(y(1,:),nfft);

%  Ymagn = abs(Y);

%  fmax = fs/2;

%  F1 = linspace(0,fmax,0.5*nfft);

%% Pwelch

%[Pxx2,F] = pwelch(sig_burst(:,1),hann(size(y,2)/2),65,nfft,fs);

%% Axes position

space = 0.01; % between graphs

margin = 0.1;

ax3 = [2*margin 2*margin 1-3*margin 0.75-2*margin-space];

ax1 = [2*margin 2*margin+ax3(4)+space ax3(3) 1-3*margin-space-ax3(4)];

%-ax2 = [margin margin 0.15 0.6];

%ax1 = [margin+ax2(3)+space margin+ax2(4)+space 1-2*margin-space-ax2(3) 1-margin-space-ax2(4)-0.08];

%ax3 = [margin+ax2(3)+space margin 1-2*margin-ax2(3)-space ax2(4)];

%% Figure finale

figure(sigtemp)

title('AR with filter Fc=2.2');

% Subplot signal

h1 = subplot(2,1,1);

plot(tt(deb:deb+per_ind)-tt(deb),y(1,:),'k','linewidth',1)

xlim([tt(deb) tt(deb+per_ind)]-tt(deb))

ymin = min(y(1,:));

ymax = max(y(1,:));

ampli = abs(ymin-ymax);

yliminf = ymin-0.1*abs(ampli);

ylimsup = ymax+0.1*abs(ampli);

ylim([yliminf ylimsup])

%title(['AR analysis - Windows length = ',num2str(seg_length), '  ;  Signal length = ',num2str(size(y,2)),'  ;  AR order = ',num2str(pAR)],...

%    'FontSize',10,'FontName','FreeSans','FontWeight','bold')

%xlabel('Temps (s)')

ylabel('P_{norm}','FontWeight','bold','fontsize',font_size_label,'fontname',font_name)

hold on

% Mise en Ã©vidence des segments

for k=1:nbr_seg

plot([ttloc((k-1)*seg_length+1) ttloc((k-1)*seg_length+1)],[yliminf ylimsup],'k--') % Extremite inf

%plot([ttloc(k*seg_length+1) ttloc(k*seg_length+1)],[yliminf ylimsup],'--','color',col2(k,:)) % Extremite sup

end

hold off

set(h1,'position',ax1,'XTick',[],'fontsize',16,'fontname',font_name)

% h2 = subplot(2,2,3);

fPxxMaxamp = abs(max(fPxxMax)-min(fPxxMax));

fliminf = 1E-20; %min(fPxxMax)-0.4*abs(fPxxMaxamp);

if fliminf <= 0

fliminf = 1;

end

flimsup = 5; %max(fPxxMax)+0.4*abs(fPxxMaxamp);

% indtemp = find(F<fliminf);

% indliminf = indtemp(end);

% indtemp = find(F<flimsup);

% indlimsup = indtemp(end);

% %semilogx(Ymagn(indliminf:indlimsup),F(indliminf:indlimsup))

% semilogx(Pxx2(indliminf:indlimsup),F(indliminf:indlimsup))

%

% ylim([fliminf flimsup])

% set(gca,'XDir','reverse','Ygrid','on','YTickLabel',{});

% %xlabel('Spectre d''amplitude','FontWeight','bold')

% xlabel('PSD (Hanning)','FontWeight','bold')

% set(h2,'position',ax2)

h3 = subplot(2,1,2);

for k=1:nbr_seg

hold on

plot([ttloc((k-1)*seg_length+1) ttloc((k-1)*seg_length+1)],[fliminf flimsup],'k--') % Extremite inf

plot([ttloc((k-1)*seg_length+1) ttloc(k*seg_length+1)],[fPxxMax(k) fPxxMax(k)],'k','LineWidth',5.4)

hold off

end

ylim([fliminf flimsup])

xlim([tt(deb) tt(deb+per_ind)]-tt(deb))

xlabel('Time (s)','FontWeight','bold','fontsize',font_size_label,'fontname',font_name)

ylabel('Strouhal Number','FontWeight','bold','fontsize',font_size_label,'fontname',font_name)

set(gca,'YGrid','on');

set(h3,'position',ax3,'fontsize',16,'GridLineStyle','--','box','on','fontname',font_name)

%[~,title_handle] = suplabel(['AR analysis - Windows length = ',num2str(seg_length), '  ;  Signal length = ',num2str(size(y,2)),'  ;  AR order = ',num2str(pAR)],'t');

%set(title_handle,'FontSize',10,'FontWeight','bold','fontname',font_name)

% A4 format

% xSize = 20.8; ySize = 29.5; % These are my size variables, width and height in cm

% xLeft = (21-xSize)/2; yTop = (30-ySize)/2; % Additional coordinates to center the figure on A4-paper

% set(gcf,'PaperUnits','centimeters') %This sets the units of the current figure (gcf = get current figure) on paper to centimeters

% set(gcf,'PaperPosition',[xLeft yTop xSize ySize])

if save_trace

print('-dpng','-r400',['/home1/emt2/dszubert/Signaux/Code_segmentation/AR/seg_length/',num2str(seg_length,'%06d'),'-',num2str(pAR,'%03d'),'.png'])

end

%clear yseg Pxx f fPxxMax

% figure(5)

% %

% % plot(Pxx2,F)

% % hold on

% % plot([min(Ymagn) max(Ymagn)],[F(indliminf) F(indliminf)],'--r')

% % plot([min(Ymagn) max(Ymagn)],[F(indlimsup) F(indlimsup)],'--r')

% % set(gca,'XScale','log','YScale','log','XDir','reverse','Ygrid','on');

% % hold off

% [cwt_data pseudofreq]=monitorCWT_simple(monitor(:,3),monitor(:,11)-mean(monitor(:,11)),'cmor1.5-1',20:200,'contourf');

% title('frequency modulation by wavelets');

% axis([0 80 pseudofreq(end) 3]);

% pseudofreq(end)

%---------------------------------------------------------------------------------------------------------------------------------------------------------

%-----------------------------------------------------------------------high-pass filter------------------------------------------------------------

%---------------------------------------------------------------------------------------------------------------------------------------------------------

function [burst] = burst_process(sig,Fc,tt)

close all

FlagPlot = 1;

FlagSavePlot = 0;

Fsize = 16; % taille titre graphs

% nbr_data = length(point);

% % frequence echantillonnage

% te = data_temp(2,1)-data_temp(1,1);

% fe = 1/te;

% % recherche frequence buffeting

% fenetre_freq = [70 90]; % freq buff = 78.8 Hz environ

% pmaxfft = monitormaxfft(k,point,data_fft,fenetre_freq);

% f_buff = sum(pmaxfft(:,3))/nbr_data;

% f_buff_std = std(pmaxfft(:,3));

% disp(['Frequence moyenne : ',num2str(f_buff), ' Hz (dispersion : ',num2str(f_buff_std),' Hz sur ',num2str(nbr_data),' valeurs)'])

% per_buff = 1/f_buff;

% clear pmaxfft

%% FFT

sig=sig-mean(sig);

yy=sig/sqrt(var(sig));

YY=fft(yy);

fe = 1/(tt(2)-tt(1));

dF=fe/length(yy);% length(yy)=length(a)

%F=0:dF:dF*length(yy)-dF;

%% EXTRACTION DES populations HF et BF BURSTS

disp('Extraction des populations HF et BF des bursts')

%     load('monitor_bursts_ifft.mat', 'tt', 'tab_sig_burst', 'F', 'tab_fft_burst', 'F_c_tremble', 'per_tremble')

%     load monitor_bursts_ifft tt tab_sig_burst F tab_fft_burst F_c_tremble per_tremble

%load tremble_burst_ifft  tab_sig tab_sig_burst tab_tremble tab_burst tab_fft_sig tab_fft_tremble tab_fft_burst t F rg_fin_tremble F_c_tremble

% DECOUPE (pour ne garder que la population HF BURSTS)

%     Fc = 58/per_tremble;  % 4.5726e+03

deb_burst = fix(Fc/dF)-3; % pour avoir var_i =0

fin_burst = length(sig)-deb_burst+3; % pour avoir var_i =0

disp(['  separation bur./burst - frequence de coupure = ',num2str(Fc),' Hz'])

%rg_fin_HF = fin_burst; % 866;

%     Ncc = size(sig,2);

%     for t=1:Ncc;

YY_c = YY;

YY_c([1:deb_burst-1 fin_burst:length(YY_c)]) = 0;

x_c = ifft(YY_c);

burst = real(x_c);

fft_HF = YY_c;

%         tab_burst_HF = [tab_burst_HF burst_HF];

%         tab_fft_burst_HF = [tab_fft_burst_HF YY_c];

%

%         Controle du caractere reel du signal

%         var_i = var(imag(x_c))

%     % Separation par DEFLATION pour ne garder que la population BF des BURSTS

%

%         fft_burst = tab_fft_burst(:,t);

%

%         YY_c = fft_burst-fft_HF;

%         x_c = ifft(YY_c);

%         q_c = real(x_c);

%         burst_BF = q_c;

%

%         tab_burst_BF = [tab_burst_BF burst_BF];

%         tab_fft_burst_BF = [tab_fft_burst_BF YY_c];

% Controle du caractere reel du signal

% var_i = var(imag(x_c))

%     end

%     disp('Save monitor_bursts_ifft.mat...')

%     save monitor_bursts_ifft tt tab_sig_burst tab_burst_HF tab_burst_BF F tab_fft_burst tab_fft_burst_HF tab_fft_burst_BF F_c_tremble Fc per_tremble

%save tremble_burst_ifft  tab_sig tab_sig_burst tab_burst_HF tab_burst_BF tab_fft_sig tab_fft_tremble tab_fft_burst tab_fft_burst_HF tab_fft_burst_BF t F rg_fin_tremble F_c_tremble rg_fin_HF F_c_HF

%% TRACES SEPARES

disp('Traces separes')

if FlagPlot

Ncc = size(sig,2);

for t=1:Ncc;

% Le burst BF

%             figure

%             plot(tt,tab_burst_BF(:,t))

%             v=axis;

%             v(3)=1.2*v(3);

%             v(4)=1.2*v(4);

%             axis(v)

%             title(['La partie BF du burst a F_c=', sprintf('%.0f Hz',F_c_HF),' au point ', num2str(t),' '],'fontsize',Fsize,'FontWeight','bold')

%             xlabel('Time (s)','fontsize',Fsize*0.8,'FontWeight','bold')

%             grid on

%

%             if FlagSavePlot

%                name = ['Sig_partie_BF_burst_m_ech_point_',num2str(t),'.jpg'];

%                print('-djpeg',name)

%             end

% Le burst HF

figure

plot(tt,burst(:,t))

title({'HF part of signal';['F_c = ', sprintf('%.2f ',Fc)]},'fontsize',Fsize,'FontWeight','bold')

xlabel('Time (s)','fontsize',Fsize*0.8,'FontWeight','bold')

ylabel('pressure variation (pa)')

grid on

if FlagSavePlot

name = ['Sig_partie_HF_burst_m_ech_point_',num2str(t),'.jpg'];

print('-djpeg',name)

end

end % for t

end % TRACES SEPARES