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
display('Please select the file')
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
% zeropadding with 2^(nextpow2(n)+6)
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
display('Please select the file')
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
display('Please select the file')
monitor=uigetfile
monitor=importdata(monitor);
[burst]=burst_process(monitor(:,11)-mean(monitor(:,11)),2.2,monitor(:,3));
% %load('D:\bei\fft-AR\burst_process.m','tt','sig_burst','burst_HF','burst_BF')
% %load('/home1/emt2/dszubert/Simu/sa_2D_dense_monit/monitor_point/resultats/monitor_sig_defl_222.mat','tt','sig_defl')
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_sig.mat','tab_sig','tab_fft_sig')
% load('monitor_bursts_ifft.mat', 'tt', 'tab_sig_burst', 'F', 'tab_fft_burst', 'F_c_tremble', 'per_tremble')
% load monitor_sig_ifft tab_sig tab_fft_sig
% 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