Question

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

以一个实际工业电力系统为例，列举了电力传动装置中的相电流奇次谐波测量平均值，基波为50HZ。

表为谐波信号参数
谐波次数 幅值（A） 相位（°）
1 380.4 0
3 2.035 60
5 3.156 135
7 1.042 157.5
9 4.291 0
11 16.13 60
13 2.016 0

根据表设信号表达式如图
（谐波次数记为n，n=2k+1）

1. 请问，如何用MATLAB进行仿真
仿真出 谐波电流信号图

2. 当取采样频率为3000HZ，数据长度为1024采样点时，利用FFT编写MATLAB语言程序得到信号频谱
请问，是怎么得出的谐波幅度谱

如果可以的话，请用word发给我，我的邮箱是847787358@qq.com
不胜感激！！

Answer 1

我来回答吧
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] and x[k]>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) && A(k)>A(k+1)
Amax=A(k);
index=k;
end
end

Answer 2

不太懂了，太专业了。