%%FFT

f0=50;

X=fft(x);%x是訊號離散值

N=length(x);

X=X(1:N/2);

Xabs=abs(X);

fori=1:m%m是諧波次數

[Amax,index]=TriFind(Xabs,floor((i*f0-15)/fsN),ceil((i*f0+15)/fsN));%如果只求基波可以改成求最大值的[Amax,index]=max(Xabs);

if(index==-1)

Fn(i)=0;

An(i)=0;

Pn(i)=0;

else

if(Xabs(index-1)>Xabs(index+1))

a1=Xabs(index-1)/Xabs(index);

r1=1/(1+a1);

k01=index-1;

else

a1=Xabs(index)/Xabs(index+1);

r1=1/(1+a1);

k01=index;

end

Fn(i)=(k01+r1-1)*fs/N;%頻率校正fs為取樣頻率

An(i)=2*pi*r1*Xabs(k01)/(N*sin(r1*pi));%幅值校正

Pn(i)=phase(X(k01))-pi*r1;%相位校正

%Pn(i,3)=mod(Pn(i),pi);

end

end

這是我寫的一個FFT求取諧波引數的程式，Pn中存的就是各次諧波的相位，要是做相位差，分別對兩個訊號求取相位，再差就行了;如果是非同步取樣，可能誤差較大，需要加窗改進；

有問題可以追問。