
How to remove unwanted region in a graph
    2 views (last 30 days)
  
       Show older comments
    
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

0 Comments
Accepted Answer
  Mathieu NOE
      
 on 3 May 2024
        
      Edited: Mathieu NOE
      
 on 3 May 2024
  
      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 Comment
  Mathieu NOE
      
 on 3 May 2024
				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');
More Answers (0)
See Also
Categories
				Find more on Annotations in Help Center and File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

