Simulink中powergui的FFT分析怎么用m文件实现？

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;
    if cycles/f0>max(time)
        error("基频太小或周期数太多")
    end
    [~,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(P1(1: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*100/P(1));
    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中明确给出了谐波信息。