www.pudn.com > chebshev.zip > chebshev2.m, change:2013-04-12,size:1185b


function [ddr,hour1,hour2,hour3] = chebshev2 
% n介切比雪夫多项式,由拟合阶段的长度c决定, 
clc; 
clear; 
n=34;           %拟合阶数,n=10+(终止时间-初始时间) 
fprintf('\nStart processing:\n'); 
 
[GCmin,GPS_Time,S_PRN,data_SP3X,data_SP3Y,data_SP3Z,SP3clk]=refrinex; 
SP3time=GCmin;        SP3X=data_SP3X;       SP3Y=data_SP3Y;       SP3Z=data_SP3Z;  
 
t0=SP3time; 
x0=SP3X;        %x0为以分钟为单位的时间序列,y0为坐标 
y0=SP3Y; 
z0=SP3Z; 
 
t=min(t0):15:max(t0);       %按15分钟进行插值计算,时间间隔可调,单位分钟 
nn=length(t0);   %历元数 
m=length(t); 
c=t0(end)-t0(1);   %c:区间拟合长度,x0(1)起始历元时刻 
 
r0=2*(t0-t0(1))/c-1;         %变量变为(-1,1)区间 
T=zeros(n,n); 
for i=1:nn 
	T(i,1)=1; 
	T(i,2)=r0(i); 
	for j=3:n 
		T(i,j)=2*r0(i)*T(i,j-1)-T(i,j-2); 
	end 
end 
 
C1=inv(T'*T)*T'*x0;       %最小二乘求切比雪夫多项式系数 
C2=inv(T'*T)*T'*y0;  
C3=inv(T'*T)*T'*z0;  
%*————————求未知历元拟合坐标—————————— 
X=zeros(1,n);x=[];y=[];z=[]; 
 
for i=1:m 
    k=2*(t(i)-t0(1))/c-1; 
	X(1)=1; 
	X(2)=k; 
	for j=3:n 
		X(j)=2*k*X(j-1)-X(j-2); 
    end 
    x=[x;X*C1]; 
    y=[y;X*C2]; 
    z=[z;X*C3]; 
end 
save('doucment2','x','y','z'); 
format long g