How to plot the output of multiple for loops?

Hi everyone,

My script works well and give output as print in the end. However, i want to plot the output crossponding to the successful iteration number. Here is my script:

clear all
clc

% BELOW SCRIPT IS FOR MULTIPLE DAYS CALCULATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd_ev=readmatrix('CE_F.csv'); % 
ev_time=datenum(cd_ev(:,1),cd_ev(:,2),cd_ev(:,3),cd_ev(:,4),cd_ev(:,5),cd_ev(:,6));
CE_M=ev_time';
for jj=1:15
     b=CE_M(:,jj);
     aa(jj)= addtodate(b, 10, 'day');
     bb(jj)= addtodate(b, -10, 'day');
end 
CE_U=aa;
CE_L=bb;
% % ------------ Data Set 2: -------------%
s_tid=dlmread('TT_test.csv', ',', 1, 1);
t1=datenum(s_tid(:,1),s_tid(:,2),s_tid(:,3),s_tid(:,4),s_tid(:,5),s_tid(:,6)); 
t_d=[t1 s_tid(:,7)];

%------------- Data Set 3: --------------% 
s_ven=dlmread('VV_test.csv', ',', 1, 1);

ts=datenum(s_ven(:,1),s_ven(:,2),s_ven(:,3),s_ven(:,4),s_ven(:,5),s_ven(:,6)); 
v_d=[ts s_ven(:,7)];

% ---------- Data selection ---------%

for kk=5:5
    s2=t_d(t_d(:,1)>= CE_L(kk) & t_d(:,1)<= CE_U(kk),:);
    s1=v_d(v_d(:,1)>= CE_L(kk) & v_d(:,1)<= CE_U(kk),:);
    ser1=s1(:, 2)-mean(s1(:,2)); % series 1 minus its mean
    ser2=s2(:, 2)-mean(s2(:,2)); % series 2 minus its mean
    ta=s1(:, 1);
    tb=s2(:, 1);
    for i=0:1:17;
             t1_f= addtodate(CE_L(kk), i, 'day');
             t2_f= addtodate(CE_L(kk), i+2, 'day');
             sf2=t_d(t_d(:,1)>= t1_f & t_d(:,1)<= t2_f,:);
             sf1=v_d(v_d(:,1)>= t1_f & v_d(:,1)<= t2_f,:);
             if length(sf1) ~= length(sf2)
                continue;
            end
            S=[sf1 sf2];
            tfa=sf1(:, 1);
            tfa=(tfa-tfa(1));  %time (days) (initial offset removed)
            sff1=S(:, 2)-mean(S(:,2)); % series 1 minus its mean
            sff2=S(:, 4)-mean(S(:,4)); % series 2 minus its mean

            N=length(tfa);
            dt=(tfa(end)-tfa(1))/(N-1); %sampling interval (days)
            fs=1/dt;  fNyq=fs/2;    %sampling frequency, Nyquist frequency (cycles/day)
            f1=1; f2=4;             %bandpass corner frequencies (cycles/day)
            [z,p,k]=butter(4,[f1,f2]/fNyq,'bandpass');  %4th order Butterworth bandpass filter coefficients
            [sos,g]=zp2sos(z,p,k); %sos representation of the filter
            y1=filtfilt(sos,g,sff1);  %apply zero phase filter to ser1
            y2=filtfilt(sos,g,sff2);  %apply zero phase filter to ser2


            %     % %--------- Zero crossing time --------%
% 
            zci = @(v) find(diff(sign(v))>0);  %define function: returns indices of +ZCs
            izc1=zci(y1);                      %find indices of + zero crossings of y1
            izc2=zci(y2);                      %find indices of + zero crossings of y2

            ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
            ZT1 = ZeroX(tfa(izc1),y1(izc1),tfa(izc1+1),y1(izc1+1));
            ZT2 = ZeroX(tfa(izc2),y2(izc2),tfa(izc2+1),y2(izc2+1));

                if length(ZT1)==length(ZT2)
            tzcd=ZT1-ZT2;          %zero crossing time differences
                elseif length(ZT1)==length(ZT2)-1
            tzcd=ZT1-ZT2(1:end-1); %zero crossing time differences
                elseif length(ZT1)-1==length(ZT2)
            tzcd=ZT1(2:end)-ZT2;   %zero crossing time differences
                end
            fdom=(length(ZT2)-1)/(ZT2(end)-ZT2(1));
            phz=360*fdom*tzcd;          %phase lag of y1 relative to y2 (degrees)
            ww=1:length(tzcd);
            ph=abs(phz);
            pha=mean(ph);
            med=median(ph);
            st=std(ph);
            fprintf('\n');
            fprintf(' %.1f\t',i, med,st,pha); 

    end

end 

fprintf command in the last line prints all the required parameters alongwith the iteration id (i). I want to 2 subplots plot (i, med) and (i, pha).

Data is also attached here.

[Geoff Hayes]

@Adnan Barkat - you could use plot and hold to plot on each successful iteration. Else store those values (from each successful iterations) into several arrays and plot those arrays outside of the for loop.

[Andi]

@Geoff Hayes, I have tried but still oonly one value

[Rik]

Use functions. Use one function to calculate what ever you want to plot for a given input. Then write a different function that will plot those data. Now you can write a third function that calls the first function several times to calculate each iteration separately.

Scripts are not for serious work, especially not with a clear all statement.

Answer 1

[Vinayak Choyyan]

Hi,

As per my understanding, you are trying to create 2 subplots of the data you have been able to successfully compute. Please try the following code:

clear all
clc
% BELOW SCRIPT IS FOR MULTIPLE DAYS CALCULATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd_ev=readmatrix('CE_F.csv'); % 
ev_time=datenum(cd_ev(:,1),cd_ev(:,2),cd_ev(:,3),cd_ev(:,4),cd_ev(:,5),cd_ev(:,6));
CE_M=ev_time';
for jj=1:15
     b=CE_M(:,jj);
     aa(jj)= addtodate(b, 10, 'day');
     bb(jj)= addtodate(b, -10, 'day');
end 
CE_U=aa;
CE_L=bb;
% % ------------ Data Set 2: -------------%
s_tid=dlmread('TT_test.csv', ',', 1, 1);
t1=datenum(s_tid(:,1),s_tid(:,2),s_tid(:,3),s_tid(:,4),s_tid(:,5),s_tid(:,6)); 
t_d=[t1 s_tid(:,7)];
%------------- Data Set 3: --------------% 
s_ven=dlmread('VV_test.csv', ',', 1, 1);
ts=datenum(s_ven(:,1),s_ven(:,2),s_ven(:,3),s_ven(:,4),s_ven(:,5),s_ven(:,6)); 
v_d=[ts s_ven(:,7)];
% ---------- Data selection ---------%
ploti=[];               %change here
plotpha=[];             %change here
plotmed=[];             %change here
for kk=5:5
    s2=t_d(t_d(:,1)>= CE_L(kk) & t_d(:,1)<= CE_U(kk),:);
    s1=v_d(v_d(:,1)>= CE_L(kk) & v_d(:,1)<= CE_U(kk),:);
    ser1=s1(:, 2)-mean(s1(:,2)); % series 1 minus its mean
    ser2=s2(:, 2)-mean(s2(:,2)); % series 2 minus its mean
    ta=s1(:, 1);
    tb=s2(:, 1);
    for i=0:1:17;
             t1_f= addtodate(CE_L(kk), i, 'day');
             t2_f= addtodate(CE_L(kk), i+2, 'day');
             sf2=t_d(t_d(:,1)>= t1_f & t_d(:,1)<= t2_f,:);
             sf1=v_d(v_d(:,1)>= t1_f & v_d(:,1)<= t2_f,:);
             if length(sf1) ~= length(sf2)
                continue;
            end
            S=[sf1 sf2];
            tfa=sf1(:, 1);
            tfa=(tfa-tfa(1));  %time (days) (initial offset removed)
            sff1=S(:, 2)-mean(S(:,2)); % series 1 minus its mean
            sff2=S(:, 4)-mean(S(:,4)); % series 2 minus its mean
            N=length(tfa);
            dt=(tfa(end)-tfa(1))/(N-1); %sampling interval (days)
            fs=1/dt;  fNyq=fs/2;    %sampling frequency, Nyquist frequency (cycles/day)
            f1=1; f2=4;             %bandpass corner frequencies (cycles/day)
            [z,p,k]=butter(4,[f1,f2]/fNyq,'bandpass');  %4th order Butterworth bandpass filter coefficients
            [sos,g]=zp2sos(z,p,k); %sos representation of the filter
            y1=filtfilt(sos,g,sff1);  %apply zero phase filter to ser1
            y2=filtfilt(sos,g,sff2);  %apply zero phase filter to ser2
            %     % %--------- Zero crossing time --------%
% 
            zci = @(v) find(diff(sign(v))>0);  %define function: returns indices of +ZCs
            izc1=zci(y1);                      %find indices of + zero crossings of y1
            izc2=zci(y2);                      %find indices of + zero crossings of y2
            ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
            ZT1 = ZeroX(tfa(izc1),y1(izc1),tfa(izc1+1),y1(izc1+1));
            ZT2 = ZeroX(tfa(izc2),y2(izc2),tfa(izc2+1),y2(izc2+1));
                if length(ZT1)==length(ZT2)
            tzcd=ZT1-ZT2;          %zero crossing time differences
                elseif length(ZT1)==length(ZT2)-1
            tzcd=ZT1-ZT2(1:end-1); %zero crossing time differences
                elseif length(ZT1)-1==length(ZT2)
            tzcd=ZT1(2:end)-ZT2;   %zero crossing time differences
                end
            fdom=(length(ZT2)-1)/(ZT2(end)-ZT2(1));
            phz=360*fdom*tzcd;          %phase lag of y1 relative to y2 (degrees)
            ww=1:length(tzcd);
            ph=abs(phz);
            pha=mean(ph);
            med=median(ph);
            st=std(ph);
            fprintf('\n');
            fprintf(' %.1f\t',i, med,st,pha);
            ploti=[ploti,i];        %change here
            plotpha=[plotpha,pha];  %change here
            plotmed=[plotmed,med];  %change here
    end
end
subplot(2,1,1);                     %change here
plot(ploti,plotmed);                %change here
subplot(2,1,2);                     %change here
plot(ploti,plotpha);                %change here

The image above shows the resulting plot. Data needed for plotting has been saved in an array and then plotted at the end of code. You can read more about ‘subplot’ here

Create axes in tiled positions - MATLAB subplot - MathWorks India