How to plot the output of multiple for loops?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Rik
le 17 Oct 2022
Modifié(e) : Rik
le 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.
Réponses (1)
Vinayak Choyyan
le 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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1162698/image.png)
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
0 commentaires
Voir également
Catégories
En savoir plus sur Line Plots dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!