单谱线插值FFT的MATAB实现

f0=50.1; % 基波频率
fs=1500; % 采样频率
N=512; % 数据长度
n=0:N-1; % 数据索引
rad=pi/180; % 角度和弧度的转换因子
xb=[1,0.02,0.1,0.01,0.05,0,0.02,0,0.01]; % 谐波幅值
Q=[-23,115.6,59.3,52.4,123.8,0,-31.8,0,-63.7]rad; % 谐波初始相位
s=zeros(1,N); % 初始化
M=9; % 谐波个数
for i=1:M % 产生谐波信号
s=s+xb(i)cos(2pif0in./fs+Q(i));
end

w=0.5-0.5cos(2pi*n./N); % 海宁窗
x=s.*w; % 信号乘以窗函数
v=fft(x,N); % FFT
u=abs(v); % 取频谱的幅值
k1=zeros(1,M); % 初始化
k2=zeros(1,M);
A=zeros(1,M);
ff=zeros(1,M);
Ph=zeros(1,M);
df=fs/N; % 频率分瓣率

for i=1:M % 计算基波和各阶谐波的参数
if i==1 % 若计算基波,在40-60Hz区间中寻找最大峰值
n1=fix(35/df); n2=fix(65/df); % 求出40Hz和60Hz对应的索引号
else % 若计算谐波,从该谐波理论值-10和+10的区间中寻找最大值
n1=fix((iff(1)-15)/df); % 求出区间对应的索引号
n2=fix((iff(1)&#