请高人帮我解决基于MATLAB的FFT谐波分析问题 请写出详细步骤 谢谢！！（这里我提出了两个问题）

我来回答吧
1.首先建立信号模型并采样
%谐波分析仿真数据生成 CreatData.m
clc
fs=3000; f0=50;
N=1024;
n=1:N;
t=(n-1)/fs;
m=13;
Am=[380.4 0 2.035 0 3.156 0 1.042 0 4.291 0 16.13 0 2.016 ];
PH=[0 0 60 0 135 0 157.5 0 0 0 60 0 0];
x=zeros(1,N);
for k=1 : m
x=x+Am(k)*cos(2*pi*f0*k*t+PH(k)/180*pi);
end;
figure(1);
plot(x);
2.进行FFT谐波分析
%直接用FFT进行谐波计算 Harmonic_Analysis.m
N=length(x);
f=50;
m=13;
n=1:N;
Fn=zeros(9,3);% 谐波频率
An=zeros(9,3);% 幅值
Pn=zeros(9,3);% 初相角
for i=1 : m
Fn(i,1) = i; Fn(i,2) = f * i;
An(i,1) = i; An(i,2) = 0;
Pn(i,1) = i; Pn(i,2) = 0;
end
fsN=fs/N;
f0=50;
%1-FFT
X=fft(x);
X=X(1:N/2);
Xabs=abs(X);
for i= 1 : m
[Amax,index]=TriFind(Xabs,floor((i*f0-15)/fsN),ceil((i*f0+15)/fsN));
if(index==-1)
Fn(i,3) = 0;
An(i,3) = 0;
Pn(i,3) = 0;
else
if(Xabs(index-1) > Xabs(index+1))
a1 = Xabs(index-1) / Xabs(index);
r1 = 1/(1+a1);
k01 = index -1;
else
a1 = Xabs(index) / Xabs(index+1);
r1 = 1/(1+a1);
k01 = index;
end
Fn(i,3) = (k01+r1-1)*fs/N;
An(i,3) = 2*pi*r1*Xabs(k01)/(N*sin(r1*pi));
Pn(i,3) = phase(X(k01))-pi*r1;
Pn(i,3) = Pn(i,3)*180/pi;
Pn(i,3) = mod(Pn(i,3),180);
end
end
Fn
An
Pn

figure(2);
stem(Xabs(1:N/2)*2/N);

dlmwrite('An.txt',An,'delimiter', '\t','precision', '%.10f','newline', 'pc');
dlmwrite('Fn.txt',Fn,'delimiter', '\t','precision', '%.10f','newline', 'pc');
dlmwrite('Pn.txt',Pn,'delimiter', '\t','precision', '%.10f','newline', 'pc');

3.完全手打的，先运行CreatData.m，再运行Harmonic_Analysis.m，应该可行，有问题再补充
4.少了个函数文件
新建个.m,命名为 TriFind，将以下内容粘贴
function [ Amax,index ] = TriFind(A,first,last)
%TRIFIND find x[k] if x[k-1]x[k+1]
% return x[k] and the k
Amax=0;
index=-1;
for k=first+1:last-1
if A(k-1)A(k+1)
Amax=A(k);
index=k;
end
end

不太懂了，太专业了。