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中明确給出了諧波資訊。