g_fitting.m

使用正交多项式完成数据拟合。程序对读入的gps采样点完成曲线拟合。

function fit

%y=csvread('weidu1.dat');
%y=csvread('jingdu1.dat');
y=csvread('DongTaiWeiDu.dat');

%对纬度数据处理
%for i=1:numel(y)
%    if y(i)<34
%        break;
%    end
%end
%y(i)=34+y(i);
for i=1:numel(y)
    y(i)=y(i)-108;
    y(i)=y(i)*100;
end

x=[0.0:1:(numel(y)-1)*1];
size(x);
size(y);
y1=zeros(size(y));%定义逼近函数的输出 
m=numel(y);%数据个数
n=60;      %拟合的阶次
%书中算法1~3步
w=zeros(1,m);
P=zeros(1,m);
T=zeros(1,m);
S=zeros(1,m);
for i=1:m
    w(i)=1;
    P(i)=1;
    T(i)=1;
end
C=zeros(n+1,1);
E=zeros(n+1,1);%求误差用的变量
Y=zeros(n,m);%求逼近函数用的变量
A=zeros(n+1,1);
b=zeros(n,1);
%书中算法4步
e=sum(w);
g=0;
f0=0;
for i=1:m
    g=g+w(i)*y(i);
    f0=f0+w(i)*y(i)*y(i);
end
%书中算法5~6步
A(1)=g/e;
de=f0-A(1)*g;
E(1)=abs(de);
%书中算法7步
for j=1:n %j相当于书中的k+1
    %保存每步的拟合结果
    for i=1:m
        Y(j,i)=A(j)*P(i);
    end
    %(1)步
    b1=0;
    for i=1:m
        b1=b1+w(i)*x(i)*(P(i)^2);
    end
    b(j)=b1/e;
    %(2)步
    for i=1:m
        P(i)=(x(i)-b(j))*T(i)-C(j)*S(i);
    end
    %(3)步
    d=0;
    for i=1:m
        d=d+w(i)*(P(i)^2);
    end
    g=0;
    for i=1:m
        g=g+w(i)*y(i)*P(i);
    end
    %(4)~(6)步
    C(j+1)=d/e;
    A(j+1)=g/d;
    e=d;
    %(7)步
    for i=1:m
        S(i)=T(i);
        T(i)=P(i);
    end
    %(8)步
    de=de-A(j+1)*g;
    E(j+1)=abs(de);
end
%函数逼近求和（按列求和）
for i=1:m
    for j=1:n
        y1(i)=y1(i)+Y(j,i);
    end
end
%分离误差数据
fp1=fopen('e:\\error.dat','wt');
for i=1:m
    pointerm(i)=y1(i)-y(i);
    fprintf(fp1,'%f,',pointerm(i));
end


subplot(2,1,1);
plot(x,y,'*');
xlabel('时间（秒）');
ylabel('经度（度）');
hold on;
plot(x,y1,'-r', 'LineWidth',2);

%画误差
%x1=[0.0:1:n];
%subplot(2,1,2);
%plot(x1,E,'-b');
subplot(2,1,2);
plot(x,pointerm,'-r');
ylabel('经度误差（度）');
xlabel('时间（秒）(a)  (b)');


%新加代码
%求均值
s=0;
for i=1:m
   s=s+pointerm(i);
end
mean=s/m;  %均值
fp1=fopen('e:\\cha.dat','wt');
%零均值化
for i=1:m
   p(i)=pointerm(i)-mean;
   fprintf(fp1,'%f,',p(i));
end


pcount=0;
dcount=0;
for i=1:m
    if pointerm(i)< mean
        dcount =dcount + 1;
    end
    if pointerm(i)> mean
        pcount = pcount+1;
    end
end

%统计游程数
dlength=0; %负游程数
plength=0; %正游程数
if p(1)< 0
   dlength=1;
end
if p(1)> 0
   plength=1;
end
for i=2:m
    if p(i)<0 && p(i-1)>0
        dlength = dlength + 1;
    end 
    if p(i) > 0 && p(i-1)<0
            plength=plength+1;
    end        
end
length=plength+dlength
count=pcount+dcount;
u=(2*pcount*dcount/count)+1
x=(2*pcount*dcount*(2*pcount*dcount-count))/(count*count*(count-1))
x=sqrt(x);
z=(length-u)/x

mb=ar(p,2,'ls');
mc=ar(p,2,'burg');
mb
mc

[a,E]=arburg(p,2)