clc;
clear?all;
close?all;
fs=1250;
N=1024;
n=0:N-1;
t=n/fs;
s=0.006*310*sin（2*pi*25*t+pi/3）+310*sin（2*pi*50*t+pi/4）+0.002*310*sin（2*pi*165*t+pi/6）+0.003*310*sin（2*pi*270*t+pi/3）+0.04*310*sin（2*pi*350*t+pi/4）+0.05*310*sin（2*pi*450*t）;
w=0.5-0.5*cos（2*pi*n./N）;%Hanning窗
r=s.*w;
v=fft（rN）;
vz=abs（v）/512*2;%Hanning窗的修正系數為2;vz存放未經插值的諧波系數
stem（abs（v）/512*2）;%繪圖
u=abs（v）;%為下面插值修正做準備，注意插值修正公式用到的是未經過幅值修正的FFT結果
A=zeros（16）;%生成1行5列0矩陣以存儲各次諧波幅值
%插值
%25Hz
y1=u（21）;
y2=u（22）;
y3=u（20）;
max=y2;
if?y3>y2
???max=y3;
????end
????if?max==y3
????????t=y1;
????????y1=max;
????????y2=t;
????end
????b=（y2-y1）/（y2+y1）;%相當于參數β
????a=1.5*b;A（1）=（y1+y2）*（2.35619403+1.15543682*a^2+0.32607873*a^4+0.07891461*a^6）/1024;
%50Hz
????y1=u（41）;
????y2=u（42）;
????y3=u（40）;
????max=y2;
????if?y3>y2
????????max=y3;
????end
????if?max==y3
????????t=y1;
????????y1=max;
????????y2=t;
????end
????b=（y2-y1）/（y2+y1）;%相當于參數β
????a=1.5*b;A（2）=（y1+y2）*（2.35619403+1.15543682*a^2+0.32607873*a^4+0.07891461*a^6）/1024;
%165
????y1=u（136）;
????y2=u（137）;
????y3=u（135）;
????max=y2;
????if?y3>y2
????????max=