matlab版本：2020b
simulink求解算法：Auto(ode3tb)，变步长运行
首先是Simulink中如何使用powergui进行FFT分析。
powergui在Simulink Library Browser中的路径为Simscape/Electrical/Specialized Power System/Fundammental Blocks。将powergui拖到Simulink中即可完成布置。
FFT在图中已用红框标出。
为了使用FFT Analysis，还需要一些其他设置。
a.获取信号
FFT中的信号来自于Scope，还需要对Scope进行设置。红框中为关键设置。变量名称随意。数据会同时存储到工作区中，也方便在工作区中进行二次分析。
2.模型参数
图中绿框给出了设置的路径。红框中不要选中。
经过设置后，重新运行仿真，即可获得数据，进入FFT分析。步骤如下。
使用powergui进行FFT分析，简单易上手，显示结果直观。而且可以看到THD。
但当需要获取一些特征频率的含量时，并不方便。因而想到使用m文件实现FFT分析。
注意到，在经过上述设置后，simulink文件运行后会在工作区中生成结构体变量。
Vabc_n内数据如下：
Vabc_n.time中存储了每个采样点的时间。Vabc_n.signals中存储了信号数据。
Vabc_n.signals.values中即为采样数据。
那么可以对该结构体变量进行FFT分析。
由于笔者对一直没有使用过fft函数。使用doc fft查看了例程。
其中关于余弦波的代码如下：
英文注释为例程中自带的注释，中文注释为笔者添加。
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
x1 = cos(2*pi*50*t); % First row wave
x2 = cos(2*pi*150*t); % Second row wave
x3 = cos(2*pi*300*t); % Third row wave
X = [x1; x2; x3];
n = 2^nextpow2(L); %为了优化fft性能，需要保证信号长度为2的幂次。
dim = 2; %对信号的每一行进行fft。默认情况下为1，即对每一列进行fft分析。
Y = fft(X,n,dim);
P2 = abs(Y/L); %双边频谱
P1 = P2(:,1:n/2+1); %单边频谱
P1(:,2:end-1) = 2*P1(:,2:end-1);
for i=1:3 %展示频谱分析结果
subplot(3,1,i)
plot(0:(Fs/n):(Fs/2-Fs/n),P1(i,1:n/2))
title(['Row ',num2str(i),' in the Frequency Domain'])
end
运行结果如下:
据此，笔者编写了如下的代码
function output = FFTofSignal(signals)
if class(signals)=="struct"
L = length(signals.time);
% 为了算法性能考虑，信号长度最好为2的幂次值
n = 2^nextpow2(L);
% 显然，n/2<L<n，截取信号长度可以为n/2或n/4时，3L可能超出索引，故取n/8。
L = n/8;
t = signals.time(L:3*L); %取信号的中间一段进行分析，这里是截取时间
Fs = (length(t)-1)/(t(end)-t(1)); %采样频率=数据点个数/时间长度
X = signals.signals.values(L:3*L,:); %截取信号
Y = fft(X,L,1); %FFT分析
P2 = abs(Y/L); %计算双侧频谱
P1 = P2(1:L/2+1,:); %单侧频谱
P1(2:end-1,:) = 2*P1(2:end-1,:);
f = Fs*(0:(L/2))/L;
output = figure;
subplot(2,1,1)
plot(t,X); %输出原信号
subplot(2,1,2)
plot(f,P1); %输出FFT分析结果
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
else
output = [];
disp("信号类型为结构体，请确认后再输入")
end
运行temp = FFTofSignal(Vabc_n);得到下图。原数据图和FFT分析结果均已放大局部区域。
可以看到分析结果并不好。原数据中谐波含量少，可以近似认为仅含有50Hz的基波分量。而FFT分析结果中，呈现了大量的其他频率。
但笔者并不清楚原因。在请教朋友后，发现了问题所在。
笔者在Simulink仿真中使用了变步长求解器，这意味着仿真得到的波形数据的采样频率不是固定的，而从fft的例程中可以看到，函数fft默认采样频率固定。实际采样频率不固定，导致了分析结果不理想。同时，DFT的理论建立在采样频率固定的基础上。如果采样频率不固定，则无法使用DFT。
为了对信号进行DFT，需要将信号采样频率固定，笔者想到了两种方法：
1.使用定步长仿真，包括连续仿真模式下的定步长模式和离散仿真。
2.对变步长仿真结果进行插值，获得定步长信号。
由于信号数据较多，为了减小运算量，采用线性插值。插值结果如下图所示。
图中从上到下依次为：
横坐标为时间，纵坐标为测量数据的原始信号
横坐标为自然数序列，纵坐标为测量数据的原始信号
横坐标为自然数序列，纵坐标为插值数据
完整代码如下
function output = FFTofSignal(t0,cycles,f0,fmax,signals)
% t0 起始时间
% f0 基频频率
% cycles 周期数
% fmax 最大频率
% signals 输入结构体信号
if class(signals)=="struct"
time = signals.time;
[~,I1] = min(abs(time-t0)); %获取与t0最近的时间点,即信号起点索引为I1；
t_end = t0 + cycles/f0;
[~,I2] = min(abs(time-t_end)); %信号终点索引为I2；
X = signals.signals.values; %截取信号
L = I2 - I1; %需要分析的信号长度。
n = 2^(nextpow2(L)-1); %为了算法性能考虑，信号长度最好为2的幂次值
t_step = (time(I2)-time(I1))/n;
t = t0:t_step:t_end;
for i = 1:length(t)
[~,I] = min(abs(time-t(i)));
X1(i,:) = X(I,:)+(X(I+1,:)-X(I,:))/((time(I+1)-time(I)))*(t(i)-time(I)); %使用线性插值获取数据点，分析得到的频谱会与实际频谱不相同,但影响应该不大。
end
Fs = (length(t)-1)/(t(end)-t(1)); %采样频率=数据点个数/时间长度
Y = fft(X1,n,1); %FFT分析
P2 = abs(Y/n); %计算双侧频谱
P1 = P2(1:n/2+1,:); %单侧频谱
P1(2:end-1,:) = 2*P1(2:end-1,:);
f = Fs*(0:(n/2))/n;
h_f = f0:f0:fmax;
for i = 1:length(h_f)
[~,I1] = min(abs(f-h_f(i)));
I(i) = I1;
end
P = P1(I,:);
THD = sqrt(sum(P(2:end,:).^2))/P(1);
THD1 = max(mean(THD),median(THD));
output = THD1;
figure;
subplot(2,1,1)
plot(time,X,'b')
hold on
plot(t,X1,'r')
hold off
cycles1 = floor(time(end)*f0);
title(['Selected signal:',num2str(cycles1),'cycles,FFT windows(in red):',num2str(cycles),'cycles'])
xlabel('Time(s)')
ylabel('Signal mag.')
subplot(2,1,2)
bar(f,P1);
title(['Fundamental(',num2str(f0),'Hz)=',num2str(P1(I(1))),',THD=',num2str(THD1*100),'%'])
xlabel('Frequency')
ylabel('Mag(% of Fundamental)')
xlim([0,fmax])
else
output = [];
disp("信号类型为结构体，请确认后再输入")
end
运行结果如下图所示。
从图中可以看到两者的分析结果相近，可以信任m文件的运行结果。
如果需要获取每一个频率的信息，可以将输出从THD更换为P1或P，P中明确给出了谐波信息。