How to remove unwanted region in a graph

9 次查看(过去 30 天)
Hi, I have this slight problem. consider the following attached code and the resultant graph also inserted as a picture. My problem is, the shoots on top of each wave is unwanted. I know they arise from points when sin(pi*t/T)=0. From my simulation, these are not important, since at these points are physically meaningless. Now how can I remove these values so the that the wave curves smoothly and peaks true maximum around the region circled(red). So that the curve smoothly peaks as illustrated with yellow marker. I thought I coluld use average around the vicinity wehere the peak appear and use it to replace the vauus that are shooting. How can I go about this?
clear
clc
r= 1e-7; %input('Enter the radius of each electron\n');
R=0.02; %input('Enter the the radius of the electron path\n');
T= 5; %input('Input the desire period\n');
t=0:0.01:30;
n=R/r;
h=6.626e-34;
me=9.109e-31;
acc1=900;
c=3e8;
f=3e15;
e=1.602e-19;
beta=1;
alpha=(h*f)/(me*c^2);
Keis=(me*c^2/(h*f))*(alpha-(alpha/(1+alpha*(1-cos(6*pi/6)))));
Z_i=beta.*sqrt(Keis);
A=Z_i*sqrt((h*f*e^2)/(R*acc1*me));
for i=1:1:numel(t)
m1=n*(abs(sin(pi*t(i)/T)));
m3=n*sqrt(abs(sin(pi*t(i)/T)));
m2=n*abs(cos(pi*t(i)/T));
if m1<2*cos(pi/6)
I(i,1)=-1000;
else
m=1:1:fix((m1)/(2*cos(pi/6)));
for j=1:1:numel(m)
A01=abs(sin(pi*t(i)/T));
A02=sqrt(A01);
I1(j,1)=A*(fix(m2/(2*A02))+fix((((m2*6)/(n*T)*180*(m(j))^2)/(A01^2*sqrt((n^2*A01^2)-(3*(m(j))^2))))+fix((m2*2*sin(acos(m(j)*sqrt(3)/(m1))))/(2*A02))));
if j==m(end);
mmax=m(end);
f10=m3*2*sin(acos((mmax*sqrt(3)/m1)));
f11=(fix(f10/2))/2;
f12=f11-fix(f11);
if f12==0
Na=1:1:((fix(f10/2)-2)/2);
for a=1:1:((fix(f10/2)-2)/2)
kon=2*180*R/(r*T);
kon1=4*180*r/(R*T);% define new constant
h=2*abs(csc(pi*t(i))/T)*Na(a)*r/R;
f13(a,1)=2*fix((2*(m2/n)*((kon*sin(acos(2*Na(a)*r/(m1*r))))+(kon1*Na(a)^2/(A01^3*sqrt(1-h^2))))));
%go to the next line
if a==numel(Na)
q=abs(csc(pi*t(i))/T)*r/R;
f14=2*fix(2*(m2/n)*((kon*sin(acos(r/(m1*r))))+(0.5*kon1*Na(a)^2/(A01^3*sqrt(1-q^2)))));
f15=A*fix((f14+sum(f13(1:a))));
I1(j,1)=f15;
end
end
else
Nb=1:1:((fix(f10/2)-1)/2);
for b=1:1:((fix(f10/2)-1)/2)
kon=2*180*R/(r*T);
kon1=4*180*r/(R*T);
h=2*abs(csc(pi*t(i))/T)*Nb(b)*r/R;
f16(b,1)=2*fix(2*(m2/n)*((kon*sin(acos(2*Nb(b)*r/(m1*r))))+(kon1*Nb(b)^2/(A01^3*sqrt(1-h^2)))));
if Nb==numel(Nb)
f17=fix(kon*(m2/n)*cos(pi*t(i)/T));
f18=A*fix(f17+sum(f16(1:b)));
I1(j,1)=f18;
end
end
end
f2=(sum(I1(1:j)));
I(i,1)=f2;
end
end
end
if i==numel(t)
for k=1:1:numel(t)
if I(k)==-1000
I(k)=max(I);%(I(2))*1;
if k==numel(t)
figure()
plot(t,I)
xlabel('time (s)')
ylabel('I(A)')
end
end
end
end
end

采纳的回答

Mathieu NOE
Mathieu NOE 2024-5-3
编辑:Mathieu NOE 2024-5-3
hello
see filloutliers , it was made for you
I changed a bit your code , I prefer to store the data (array tt and II) then plot it , instead of calling plot multiple times inside a for loop
makes also post processing easier
clear
clc
r= 1e-7; %input('Enter the radius of each electron\n');
R=0.02; %input('Enter the the radius of the electron path\n');
T= 5; %input('Input the desire period\n');
t=0:0.01:30;
n=R/r;
h=6.626e-34;
me=9.109e-31;
acc1=900;
c=3e8;
f=3e15;
e=1.602e-19;
beta=1;
alpha=(h*f)/(me*c^2);
Keis=(me*c^2/(h*f))*(alpha-(alpha/(1+alpha*(1-cos(6*pi/6)))));
Z_i=beta.*sqrt(Keis);
A=Z_i*sqrt((h*f*e^2)/(R*acc1*me));
tt = [];
II = [];
for i=1:1:numel(t)
m1=n*(abs(sin(pi*t(i)/T)));
m3=n*sqrt(abs(sin(pi*t(i)/T)));
m2=n*abs(cos(pi*t(i)/T));
if m1<2*cos(pi/6)
I(i,1)=-1000;
else
m=1:1:fix((m1)/(2*cos(pi/6)));
for j=1:1:numel(m)
A01=abs(sin(pi*t(i)/T));
A02=sqrt(A01);
I1(j,1)=A*(fix(m2/(2*A02))+fix((((m2*6)/(n*T)*180*(m(j))^2)/(A01^2*sqrt((n^2*A01^2)-(3*(m(j))^2))))+fix((m2*2*sin(acos(m(j)*sqrt(3)/(m1))))/(2*A02))));
if j==m(end);
mmax=m(end);
f10=m3*2*sin(acos((mmax*sqrt(3)/m1)));
f11=(fix(f10/2))/2;
f12=f11-fix(f11);
if f12==0
Na=1:1:((fix(f10/2)-2)/2);
for a=1:1:((fix(f10/2)-2)/2)
kon=2*180*R/(r*T);
kon1=4*180*r/(R*T);% define new constant
h=2*abs(csc(pi*t(i))/T)*Na(a)*r/R;
f13(a,1)=2*fix((2*(m2/n)*((kon*sin(acos(2*Na(a)*r/(m1*r))))+(kon1*Na(a)^2/(A01^3*sqrt(1-h^2))))));
%go to the next line
if a==numel(Na)
q=abs(csc(pi*t(i))/T)*r/R;
f14=2*fix(2*(m2/n)*((kon*sin(acos(r/(m1*r))))+(0.5*kon1*Na(a)^2/(A01^3*sqrt(1-q^2)))));
f15=A*fix((f14+sum(f13(1:a))));
I1(j,1)=f15;
end
end
else
Nb=1:1:((fix(f10/2)-1)/2);
for b=1:1:((fix(f10/2)-1)/2)
kon=2*180*R/(r*T);
kon1=4*180*r/(R*T);
h=2*abs(csc(pi*t(i))/T)*Nb(b)*r/R;
f16(b,1)=2*fix(2*(m2/n)*((kon*sin(acos(2*Nb(b)*r/(m1*r))))+(kon1*Nb(b)^2/(A01^3*sqrt(1-h^2)))));
if Nb==numel(Nb)
f17=fix(kon*(m2/n)*cos(pi*t(i)/T));
f18=A*fix(f17+sum(f16(1:b)));
I1(j,1)=f18;
end
end
end
f2=(sum(I1(1:j)));
I(i,1)=f2;
end
end
end
if i==numel(t)
for k=1:1:numel(t)
if I(k)==-1000
I(k)=max(I);%(I(2))*1;
if k==numel(t)
tt = [tt; t];
II = [II;I];
% figure()
% plot(t,I)
% xlabel('time (s)')
% ylabel('I(A)')
end
end
end
end
end
IIs = filloutliers(II,'linear','movmedian',25);
plot(tt,II,'b',tt,IIs,'r')
xlabel('time (s)')
ylabel('I(A)')
  1 个评论
Mathieu NOE
Mathieu NOE 2024-5-3
if you have an older matlab release , you can achieve the same result with interp1
load data.mat
% filloutliers METHOD
IIout = filloutliers(II,'linear','movmedian',25);
% OLD MATLAB METHOD (with interp1)
all_idx = 1:length(tt);
IIs = smoothdata(II,'movmedian',25);
outlier_idx = (abs(II./IIs) > 1.2); % Find outlier idx
IIold = II;
IIold(outlier_idx) = []; % remove outlier from dataset
IIold = interp1(all_idx(~outlier_idx), IIold, all_idx)'; % replace outliers by interpolated data
plot(tt,II,'b',tt,IIout,'*r',tt,IIold,'g')
xlabel('time (s)')
ylabel('I(A)')
legend('raw','filloutliers','interp1');

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Interpolation 的更多信息

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by