How to plot the output of multiple for loops?
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
Rik
2022-10-17
编辑:Rik
2022-10-17
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.
回答(1 个)
Vinayak Choyyan
2022-10-20
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Line Plots 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!