下面是我在之前论坛里找到的程序里进行改编的，但还是弄不来，希望懂的人能指导下

f0=50;

fs=50*256;

N=2048;%2048/256=8

rad=180/pi;

n=[0:1:N-1];

xb=[220,0.5,25,0.4,6,0.3,4,0.2,2];

Q=[pi/3,pi/3,pi/6,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3];

s=zeros(1,N);

%s=220*sin(2*pi*f0*n./fs+pi/6)+0.5*sin(4*pi*f0*n./fs+pi/5)+0.2*sin(8*pi*f0*n./fs+pi/7)+0.07*sin(16*pi*f0*n./fs+pi/5)+0.5*sin(32*pi*f0*n./fs+pi/5);

for i=1:9

s=s+xb(i)*sin(2*pi*f0*i*n./fs+Q(i));

end

w=0.338946-0.481973*cos(2*pi*n./N)+0.161054*cos(4*pi*n./N)-0.018027*cos(6*pi*n./N);

r=s.*w;

v=fft(r,N);

u=abs(v);

vf=abs(v)/1024;

%为了求出谐波，那么必须找出基波，在35-65中频谱最大值所对应的就是离基波频率最近的那个点

y=0;

for i=fix(35/(fs/N)):fix(65/(fs/N))

if u(i+1)>y

y=u(i+1);

num=i+1;

end

end

%break

%谱线最大值已经找到，下面是要确定k1和k2还有y1和y2

k1=zeros(1,9);

k2=zeros(1,9);

k0=zeros(1,9);

b=zeros(1,9);

aa=zeros(1,9);

y1=zeros(1,9);

y2=zeros(1,9);

y3=zeros(1,9);

k1(1)=num;

k0(1)=f0/(fs/N);

aa(1)=k0(1)-k1(1)-0.5;

y1(1)=u(k1(1));%

k2(1)=k1(1)+1;

y2(1)=u(k2(1));

y3(1)=u(k1(1)-1);

max=y2(1);

if y3(1)>y2(1)

max=y3(1);

end

if max==y3(1);

t=y1(1);

y1(1)=max;

y2(1)=t;

k1(1)=k1(1)-1;%根据公式，谱线谁在前谁是y1，对应的就是k1

k2(1)=k1(1)+1;%根据公式，谱线谁在后谁是y2，对应的就是k2

end

%%k1和k2还有y1和y2确定完毕

%%确定αβ，在这里用a和b表示，确定基波频率

A=zeros(1,63);

ff=zeros(1,63);

Ph=zeros(1,63);

b(1)=(y2(1)-y1(1))/(y2(1)+y1(1));

for i=2:9

k1(i)=fix(i*50/(fs/N));

k2(i)=k1(i)+1;

k0(i)=f0*i/(fs/N);

y1(i)=u(k1(i)+1);

y2(i)=u(k2(i)+1);

b(i)=(y2(i)-y1(i))/(y2(i)+y1(i));

aa(i)=k0(i)-k1(i)-0.5;

end

n=3;

p=polyfit(aa,b,n);

plot(aa,b,0:0.01:10,polyval(p,0:0.01:10))