removewilddata.m

文件包含有5项内容： 一、扩展卡尔曼滤波EKF 二、去偏转换卡尔曼滤波CMKF 三、最小二乘拟和的方法 四、最小二乘、EKF、CMKF的比较 五、野值剔除算法 用MATLAB实现了这些具体

%以雷达仿真信号来验证野值判断算法   4.2

clear

 [Z,Z1,XTrue1,totaltime,T0]=source;
 
 times=totaltime/T0+1;

for i=1:times
    if i<7
%         a(i)=Z1(2,1,i);               %  azimuth
          a(i)=Z1(1,1,i);                           % range
    elseif i==7
       
%         a(7)=Z1(2,1,7);
        a(7)=Z1(1,1,7);
        b(1)=(32*a(1)+15*a(2)+3*a(3)-4*a(4)-6*a(5)-3*a(6)+5*a(7))/42;
        b(2)=(5*a(1)+4*a(2)+3*a(3)+2*a(4)+1*a(5)-a(7))/14;
        b(3)=(a(1)+3*a(2)+4*a(3)+4*a(4)+3*a(5)+a(6)-2*a(7))/14;
        b(4)=(-a(1)+3*a(2)+6*a(3)+7*a(4)+6*a(5)+3*a(6)-2*a(7))/21;
        b(5)=(-2*a(1)+a(2)+3*a(3)+4*a(4)+4*a(5)+3*a(6)+a(7))/14;
        b(6)=(-a(1)+a(2)+2*a(4)+3*a(5)+4*a(6)+5*a(7))/14;
        b(7)=(5*a(1)-3*a(2)-6*a(3)-4*a(4)+3*a(5)+15*a(6)+32*a(7))/42;
       
    else
    
%     a(i)=Z1(2,1,i);  %  azimuth
     a(i)=Z1(1,1,i);      % range
    kk1=i-6;kk2=i-5;kk3=i-4;kk4=i-3;kk5=i-2;kk6=i-1;
    b(i)=( 5*a(kk1)-3*a(kk2)-6*a(kk3)-4*a(kk4)+3*a(kk5)+15*a(kk6)+32*a(i))/42;
    
    sum=0;
%     summ=1;
%     summ1=1;
%     SUM=0;

    for k=i-6:i
        sum=sum+(a(k)-b(k))^2;
    end
     if abs(a(i)-b(i))>2.2*sqrt(sum/6)              %原公式此处为2.2可是对一些比较突出的跳点,不能进行处理
         i
%          for j=i-7:i
%              for h=i-7:i
%                  summ=summ*(i-h);
%                  
%                  summ1=summ1*(j-h);
%                 
%               end 
%               
%               if i==j
% %                   SUM=SUM;.3
%                  else
%                  SUM=summ*a(j)/((i-j)*summ1);
%               end
%           end
%           
%           a(i)=SUM;   
        h=i-7:i-1;
        bb=[a(i-7),a(i-6),a(i-5),a(i-4),a(i-3),a(i-2),a(i-1)];
%       bb=[b(i-7),b(i-6),b(i-5),b(i-4),b(i-3),b(i-2),b(i-1)];
        p= polyfit(h,bb,2);
%       hh=1:30
      a(i)=polyval(p,i);
%       
%        a(i)=p(1)*i+p(2);
%        a(i)=p(1)*i^2+p(2)*i+p(3);
%      MyData2(i)=a(i);

       end
   end

 end
   for i=1:times 
%        C(i)=Z1(2,1,i);
  C(i)=Z1(1,1,i);

   end
   
  x=1:times;
  subplot(2,1,1)
   
    plot(x,C) ;
    title('原始数据信号')
%     ylabel('弧度');
ylabel('米');
xlabel('k');

   subplot(2,1,2)
   
 plot(x,a);
 title('去野值后数据信号')
%  ylabel('弧度');
ylabel('米');

xlabel('k');

% %