%relax譜估計算法，k為所含諧波數(shù)，可以估算，本程序假定該值已知，且由于源程序是按幅值諧波計算，故k應(yīng)為實際諧波個數(shù)2倍，且幅值減半
%th-----迭代閾值
%Nmax---迭代最大次數(shù)
function [fre,ampitud,pha]=relax(x,fs,k);
data=zeros(k,length(x));f=zeros(1,k);amp=zeros(1,k);
t=1000;
N=length(x);th=0.01;w=0:2*pi/(N-1):2*pi;
for i=1:k
    t=1000;id=0;
    while t>th | id<10
        estdata=0;
       for j=i:-1:1
           m=1:i;m(m==j)=[];
           data(j,:)=x-sum(data(m,:),1);    %data(j,:)=yi
           fftdata=fft(data(j,:));
           maxdata=max(abs(fftdata));
           wi=w(abs(fftdata)==maxdata);f(j)=wi(1);
           amptest=fftdata(abs(fftdata)==maxdata);amp(j)=amptest(1);
           data(j,:)=exp(sqrt(-1)*wi(1)*[0:N-1].*amp(j));  %data(j,:)=amp(j)*exp(j*wj*n),與上面表示含義不同
           estdata=estdata+data(j,:);
       end 
           id=id+1;
           t=norm(x-estdata);
           
    end
    end
           fre=fs*f/(2*pi);ampitud=abs(amp)/2;pha=phase(amp);