clc
clear?all
R0=1;
M=8;????????????????%諧波次數
num=0;
fs=1000;
N=512;
n=0:N-1;
T=n/fs;
f0=50;
x=220*sin（2*pi*f0*T）+110*sin（4*pi*f0*T）+100*sin（6*pi*f0*T）+90*sin（8*pi*f0*T）+105*sin（10*pi*f0*T）+50*sin（12*pi*f0*T）+20*sin（14*pi*f0*T）+20*sin（16*pi*f0*T）;%不考慮諧波
z=windows（6xN）;????%加納托爾（Nuttall）窗
y=fft（zN）;
m=abs（y）*5/N;????????%時域幅值
mag0=abs（y）;?????????%求頻域幅值
xj=angle（y）;
for?i=0:N-1
????if（（i*fs/N）>=55*（R0-1）&&（i*fs/N）<55*R0&&R0<=M）
????????mag（R0i+1-（R0-1）*num）=mag0（i+1）;
????else?
????????if?R0==0
????????????num=i;
????????end
????????R0=R0+1;
????????mag（R0i+1-（R0-1）*num）=mag0（i+1）;
????end
end??
[k1U]=max（mag[]2）;
for?R0=1:M
mag（R0U（R0））=0;
end
[k2V]=max（mag[]2）;
for?R0=1:M
B（R0）=（k2（R0）-k1（R0））/（k1（R0）+k2（R0））;
A（R0）=2.5430095*B（R0）+0.