%用正交多項式作最小二乘曲線擬合
%用正交的多項式函數系{Qj(x)(j=0,1,……,m)}作為基函數作最小二乘擬合，其系數為qj(j=0,1,……,m)
%用三個向量b,t與s分別存放多項式Qj-1(x),Qj(x)與Qj+1(x)的系數
function P=spir(X,Y,n,m)    %(X,Y)是n 個數據點的坐標,a返回m-1次擬合多項式的m個系數,m是擬合多項式的項數(多項式的次數為m-1)
for i=1:m
    a(i)=0.0;
end
if m>n
    m=n;
end

%計算X(i)平均值
z=0.0;
for i=1:n
    z=z+X(i)/(1.0*n);
end

%（1）構造Q1(x)  設Q1(x)=b1,由遞推公式,有b1=1
b(1)=1.0;
d1=1.0*n;
p=0.0;c=0.0;
for i=1:n
    p=p+(X(i)-z);
    c=c+Y(i);
end
q=c/d1; p=p/d1;  %計算q1與p1
a(1)=q*b(1);      %即得a1

%（2）構造Q2(x) 設Q2(x)=t1+t2x,由遞推公式,有t1=-p,t2=1
if m>1
    t(2)=1.0;t(1)=-p;
    d2=0.0;c=0.0;g=0.0;
    for i=1:n
        Q2=(X(i)-z)-p;d2=d2+Q2*Q2;
        c=c+Y(i)*Q2;
        g=g+(X(i)-z)*Q2*Q2;
    end
    q=c/d2; p=g/d2; B=d2/d1;  %計算q2與p2與B2
    d1=d2;
    a(2)=q*t(2); a(1)=q*t(1)+a(1);
end

%（3）對于j=3,4,……,m-1,逐步遞推Qj(x) 
for j=3:m
    s(j)=t(j-1);
    s(j-1)=-p*t(j-1)+t(j-2);
    if j>=4
        for k=j-3:1
            s(k)=-p*t(k)+t(k-1)-B*b(k);
        end
    end
    s(1)=-p*t(1)-B*b(1);
    d2=0.0;c=0.0;g=0.0;
    for i=1:n
        q=s(j);
        for k=j-1:1
            q=q*(X(i)-z)+s(k);
        end
        d2=d2+q*q; c=c+Y(i)*q;
        g=g+(X(i)-z)*q*q;
    end
    c=c/d2; p=g/d2; q=d2/d1;
    d1=d2;
    a(j)=c*s(j);t(j)=s(j);
    for k=j-1:1
        a(k)=c*s(k)+a(k);
        b(k)=t(k); t(k)=s(k);
    end
end
%P=a(1);
%x=X(1:300);
%syms x;
%for i=1:m-1
 %   for k=1:300
  %  P=P+a(i+1)*x^i;
   % end
%end
P=polyval(a,X)