relax.m

来自「谐波检测的RELAX算法」· M 代码 · 共 29 行

%relax谱估计算法，k为所含谐波数，可以估算，本程序假定该值已知，且由于源程序是按幅值谐波计算，故k应为实际谐波个数2倍，且幅值减半
%th-----迭代阈值
%Nmax---迭代最大次数
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);