加窗插值FFT源程序

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=