How to plot the output of multiple for loops?

1 view (last 30 days)
Andi
Andi on 17 Oct 2022
Answered: Vinayak Choyyan on 20 Oct 2022
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.
  3 Comments
Rik
Rik on 17 Oct 2022
Edited: Rik on 17 Oct 2022
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.

Sign in to comment.

Answers (1)

Vinayak Choyyan
Vinayak Choyyan on 20 Oct 2022
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

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by