如何用MATLAB FFT畫出PSD？

發表於 2025-10-02 更新於 2025-10-03 分類於通訊技術 Disqus：

Why?

在MATLAB裡產生的訊號都是finite length、finite amplitude的，因此嚴格來說都是discrete-time energy signal。但是模擬power signal或是random process時，往往需要計算這些訊號的power spectrum或是power spectral density (PSD)。這篇文章說明用FFT畫出各種spectrum的作法。

How?

MATLAB裡的 fft 函數計算某個sequence的DFT，如下：

$$
X[k] = \sum_{n = 0}^{N_{fft} - 1}{x[n] e^{-j\frac{2\pi}{N_{fft}} kn}},\quad 0 \leq k \leq (N_{fft} - 1)
$$

其中 $N_{fft}$ 是DFT的點數，可以是sequence length，也可以自行指定，fft 函數會幫你做clipping或zero padding。由於算出來的頻譜，頻率是從 0 到 $2\pi$，也就是從低頻到高頻 ($\pi$) 再到低頻，fft 算出來的結果還要透過 fftshift 函數來把低頻移到中間、高頻移到兩邊。IDFT則是由 IFFT 函數實現：

$$
x[n] = \frac{1}{N_{fft}}\sum_{k = 0}^{N_{fft} - 1}{X[k] e^{j\frac{2\pi}{N_{fft}} kn}},\quad 0 \leq n \leq (N_{fft} - 1)
$$

實際上，sequence做DFT得到的結果，它的物理意義是什麼？我認為DFT大概呈現了sequence的頻譜成份，比如說若高頻的amplitude較高，則訊號的高頻成分較多。但是這個「多」是用什麼來衡量的？要更精確地描述DFT的結果，首先要確定這個sequence用來近似的訊號的類型。

Energy singal $\longrightarrow$ ESD
Periodic signal $\longrightarrow$ Power spectrum
Power signal/random process $\longrightarrow$ PSD

Energy Spectral Density (ESD)

Energy signal可以是finite length或是infinite length(如sinc函數)，這裡只考慮finite length。假設一finite length、finite amplitude的energy signal $x(t)$，在 $(0, T)$ 的區間非零，其餘時刻皆為零。這個訊號的總能量為：

$$
\begin{align}
E_{total} &= \int_{0}^{T}{\lvert x(t) \rvert^2 dt}\\
&= \int_{-\infty}^{\infty}{\lvert X_c(f) \rvert^2 df}\quad\text{(Parseval’s theorem)}
\end{align}
$$

可知 $\lvert X_c(f) \rvert^2$ 就是能量頻譜密度，ESD。我們將 $x(t)$ 以取樣頻率 $f_s = 1/T_s$ 取樣，得到sequence $x[n]$：

$$
x[n] = x(nT_s),\quad 0 \leq n \leq (T_{tot}/T_s)
$$

其中 $T_{tot}$ 為取樣時間總長度，$N = T_{tot}/T_s$ 為sequence length。

$x[n]$ 也是energy signal。根據impulse invariance的性質，$X(e^{j\omega}) = X_c(j\omega/T_s)/T_s$，以及 $\lvert X_c(j\Omega)\rvert^2 = \lvert X_c(f)\rvert^2,\quad\text{where } \Omega = 2\pi f$，得到

$$
\lvert X(e^{j\omega})\rvert^2 = \frac{1}{(T_s)^2} \lvert X_c(f)\rvert^2
$$

不過，我們透過 fft 得到的頻譜不是 $X(e^{j\omega})$，而是對 $X(e^{j\omega})$ 取樣的 $X(e^{j\omega})\rvert_{\omega = 2\pi k/N} = X[k]$，代入剛才的式子：

$$
\begin{align}
&\lvert X[k]\rvert^2 = \frac{1}{(T_s)^2} \lvert X_c(\frac{kf_s}{N})\rvert^2,\quad 0 \leq k \leq (\frac{N}{2} - 1)\\
&\lvert X[k]\rvert^2 = \frac{1}{(T_s)^2} \lvert X_c(\frac{(k - N)f_s}{N})\rvert^2,\quad \frac{N}{2} \leq k \leq (N - 1)
\end{align}
$$

所以如果我們想要呈現 $\lvert X_c(f)\rvert^2$ (或者 $X_c(f)$，一樣的意思)，就要對每個 $X[k]$ 做以下的scaling：

$$
X’[k] = T_s X[k] = \frac{1}{f_s} X[k]
$$

以 $x(t) = \sin{(2\pi t)},\quad 0 \leq t \leq 1$ 為例，可以看做是一個sine函數乘上一個square pulse，所以理論上 $X_c(f)$ 為

$$
\lvert X_c(f) \rvert = \frac{1}{2} \lvert sinc(f - 1) - sinc(f + 1)\rvert
$$

Sine乘上方波，增加FFT的長度以增加頻域解析度。

MATLAB程式碼如下，取樣頻率為10 kHz。下圖則是 $\lvert X_c(f)\rvert$ 和經過scaling的 fft 結果比較(取絕對值)。

fs = 10000; % sampling freq = 10 kHz
t1 = (0: fs - 1)/fs;
x1 = sin(2*pi*t1);
x = [zeros(1, round(4.5*length(t1))) x1 zeros(1, round(4.5*length(t1)))];
Nsample = length(x); % number of points of FFT
t = (-0.5*Nsample: 0.5*Nsample - 1)/fs;

% time domain
figure();
plot(t, x, "b-", "LineWidth", 2);
xlabel("time (sec)");
ylabel("x(t)");

% freq domain
figure();
f = (-0.5*Nsample: 0.5*Nsample - 1)*fs/Nsample;
X = fftshift(fft(x))/fs;
plot(f, abs(X), "b-", "LineWidth", 2);
hold on
plot(f, 0.5*abs(sinc(f - 1) - sinc(f + 1)), "r-", "LineWidth", 2);
hold off
legend("scaled fft", "theo");
xlabel("freq (Hz)");
ylabel("X_c(f)");
xlim([-6 6]);

兩條曲線重疊在一起，幾無差異。

除了amplitude之外，還要關注頻譜的總能量。理論上，$x(t)$ 的總能量為0.5(以剛才的例子來說)。我們可以用MATLAB嘗試計算取樣後訊號的總能量：

$$
\sum_{n = 0}^{N - 1}{\lvert x[n]\rvert^2} = 5000
$$

注意到上面的式子，並不是Riemann sum，而這麼做是錯的(所以算出來的結果是錯的)，因為計算總能量應該要用積分，對應到sequence就應該使用Riemann sum。因此正確的總能量應該這樣計算：

$$
\begin{align}
E_{total} &= \sum_{n = 0}^{N - 1}{\lvert x[n]\rvert^2 \frac{T_{tot}}{N}}\\
&= \sum_{n = 0}^{N - 1}{\lvert x[n]\rvert^2 T_s}
\end{align}
$$

又根據DFT版本的Parseval’s theorem：

$$
\sum_{n = 0}^{N - 1}{\lvert x[n]\rvert^2} = \frac{1}{N} \sum_{k = 0}^{N - 1}{\lvert X[k]\rvert^2}
$$

可知總能量也可以表示為：

$$
\begin{align}
E_{total} &= \frac{T_s}{N} \sum_{k = 0}^{N - 1}{\lvert X[k]\rvert^2}\\
&= \sum_{k = 0}^{N - 1}{\frac{1}{f_s N} \lvert X[k]\rvert^2}\\
&= \sum_{k = 0}^{N - 1}{\frac{f_s}{N} \lvert \frac{1}{f_s} X[k]\rvert^2}
\end{align}
$$

最後一行特別有意思，因為它正是頻域的Riemann sum！藉此也驗證了對FFT頻譜的scaling不僅保留amplitude，也保留總能量。

Power Spectrum

這裡的power spectrum指的是週期性訊號(因此也是power signal)的功率頻譜，看起來應該會是一根一根的，每一根的單位是Watt，而不是Watt/Hz(功率密度單位)。週期性訊號的功率這樣計算：

$$
P_{total} = \frac{1}{T} \int_{0}^{T}{\lvert x(t)\rvert^2 dt}
$$

$T$ 為訊號的週期。我們知道週期性訊號的Fourier series表示為

$$
x(t) = \sum_{k = -\infty}^{\infty}{C_k e^{j2\pi kf_0 t}}
$$

因此總功率的計算可以改寫為

$$
\begin{align}
P_{total} &= \frac{1}{T} \int_{0}^{T}{(\sum_{k = -\infty}^{\infty}{C_k e^{j2\pi kf_0 t}})(\sum_{l = -\infty}^{\infty}{C_l^{*} e^{-j2\pi lf_0 t}}) dt}\\
&= \sum_{k = -\infty}^{\infty}{\sum_{l = -\infty}^{\infty}{C_k C_l^{*} \delta[k - l]}}\\
&= \sum_{k = -\infty}^{\infty}{\lvert C_k\rvert^2}
\end{align}
$$

因此只要知道 $C_k$ 就可以畫出power spectrum。現在我們要找到FFT頻譜和 $C_k$ 的關係。假設取樣頻率為 $f_s$，取樣週期 $T_s$。用MATLAB來呈現週期性訊號有兩種方式：

取一個週期的訊號做FFT，相當於DFT頻譜即為DTFS的情況，當取樣點數夠多，$X[k]$ 就能夠充分表示 $C_k$ (比如說取樣頻率為基頻的整數倍，且滿足Nyquist sampling criterion，每個 $C_k$ 都會被表示到)。先從時域來計算總功率：

$$
\begin{align}
P_{total} &= \frac{1}{NT_s} \sum_{n = 0}^{N - 1}{\lvert x[n]\rvert^2 T_s}\\
&= \frac{1}{N} \sum_{n = 0}^{N - 1}{\lvert x[n]\rvert^2}
\end{align}
$$

因此

$$
\begin{align}
P_{total} &= \frac{1}{N^2} \sum_{k = 0}^{N - 1}{\lvert X[k]\rvert^2}\\
&= \sum_{k = 0}^{N - 1}{\lvert \frac{1}{N} X[k]\rvert^2} \approx \sum_{k = -\infty}^{\infty}{\lvert C_k\rvert^2}
\end{align}
$$

故

$$
X’[k] = \frac{1}{N} X[k]
$$

$N$ 為sequence length，也是FFT點數。
取數個週期的訊號做FFT，假設取了M個週期。一樣從時域來計算總功率：

$$
\begin{align}
P_{total} &= \frac{1}{MN_TT_s} \sum_{n = 0}^{MN_T - 1}{\lvert x[n]\rvert^2 Ts}\\
&= \frac{1}{MN_T} \sum_{n = 0}^{MN_T - 1}{\lvert x[n]\rvert^2}\\
&= \frac{1}{N} \sum_{n = 0}^{N - 1}{\lvert x[n]\rvert^2}\\
\end{align}
$$

其中 $N_T$ 為一個週期的取樣點數。這個結果和取一個週期計算的結果是一樣的，所以作法也一樣。

那麼，如果取樣的時間長度並非週期的整數倍呢？答案是DFT頻譜就無法準確反映週期訊號的power spectrum。如果無法取週期的整數倍，至少取較長的一段時間，如週期的100.3倍等等，可以減少頻譜的誤差。

以週期性訊號 $x(t) = \sin{(200\pi t) + 2\cos{(600\pi t)}}$ 為例，下面兩張圖，第一張圖拿一個週期的訊號做FFT，呈現完美的power spectrum；第二章圖拿10.5個週期的訊號做FFT，可以看到power spectrum出現一些誤差。

取一個週期做FFT。

取10.5個週期做FFT。

Power Spectral Density (PSD)

PSD大概是最難以呈現的，因為PSD描述的是非週期性的power signal或是random process，但是MATLAB只能寫出deterministic、finite length的訊號。假設一random process $X(t)$，其autocorrelation為 $R_X(\tau)$。則PSD可以對 $R_X(\tau)$ 做Fourier transform得到：

$$
S(f) = \int_{-\infty}^{\infty}{R_X(\tau) e^{-j2\pi f\tau} d\tau}
$$

很遺憾的，我們只能讓MATLAB生成sample function，所以 $R_X(\tau)$ 必須用time-average autocorrelation來估計，再拿來做Fourier transform。不過，我們也可以從PSD的定義找到線索。

根據Ziemer和Tranter的Principles of Communications，PSD可以定義為

$$
S(f) = \lim_{T \rightarrow \infty}{\frac{E[\lvert F\{X_{2T}(t)\}\rvert^2]}{2T}}
$$

其中 $F\{\cdot\}$ 為Fourier transform，$E[\cdot]$ 為ensemble average。由於我們只能生成deterministic signal，因此最外層的 $E[\cdot]$ 拿掉；再來，我們只能生成有限長度的sequence，故 $\lim$ 也拿掉。用 $\hat{S}(f)$ 來表示對 $S(f)$ 的估計：

$$
\begin{align}
\hat{S}(f) &= \frac{\lvert F\{x_{2T}(t)\}\rvert^2}{2T}\\
&= \frac{\lvert F\{x(t)\}\rvert^2}{T’}\\
&= \frac{\lvert X_c(f)\rvert^2}{T’}
\end{align}
$$

其中 $T’$ 為訊號 $x(t)$ 的長度。以 $\hat{S}(f)$ 估計總功率：

$$
\hat{P} = \int_{-\infty}^{\infty}{\hat{S}(f) df} = \int_{-\infty}^{\infty}{\frac{\lvert X_c(f)\rvert^2}{T’} df}
$$

將上式寫成Riemann sum的形式：

$$
\begin{align}
\int_{-\infty}^{\infty}{\hat{S}(f) df} &= \int_{-\infty}^{\infty}{\frac{\lvert X_c(f)\rvert^2}{T’} df}\\
&\approx \sum_{k = 0}^{N - 1}{\frac{\lvert \frac{X[k]}{f_s}\rvert^2}{NT_s}\cdot\frac{f_s}{N}}\\
&= \sum_{k = 0}^{N - 1}{\lvert \frac{1}{\sqrt{Nf_s}}X[k]\rvert^2\cdot\frac{f_s}{N}}
\end{align}
$$

當 $N$ 很大，上式便成立。上式將 $f_s/N$ 放在平方項外面，因為平方項的意義為功率「密度」，因此和ESD的情況相同，$f_s/N$ 作為頻率軸的微小變化量。最後得到對 $X[k]$ 的scaling constant：

$$
X’[k] = \frac{1}{\sqrt{Nf_s}} X[k]
$$

以Jake’s model的Doppler spectrum為例，若

$$
\begin{align}
E[\lvert X\rvert^2] &= E[{X_i}^2 + {X_q}^2] = 1\\
&= R_X(0) = \int_{-\infty}^{\infty}{S_D(f) df}
\end{align}
$$

當中 $X_i$ 為fading coefficient的實部，$X_q$ 則是fading coefficient的虛部，兩者為互相獨立的zero-mean Gaussian process。則理想的Doppler spectrum為

$$
S_D(f) = \frac{1}{\pi f_m\sqrt{1 - (\frac{f}{f_m})^2}}
$$

下圖為理想的Doppler spectrum和Jake’s model的Doppler spectrum比較。可以看到在頻譜U型區域附近，模擬值在理論值上下波動，足以驗證FFT scaling的正確性。

Jake's model Doppler spectrum。

Summary

最後總結一下將FFT轉換為能量/功率頻譜的作法，$X[k]$ 表示 fft 函數的輸出：

ESD

$$
X’[k] = \frac{1}{f_s} X[k]
$$
Power spectrum

$$
X’[k] = \frac{1}{N} X[k]
$$

其中 $N$ 為sequence length，也是FFT點數。
PSD

$$
X’[k] = \frac{1}{\sqrt{Nf_s}} X[k]
$$

$N$ 為sequence length，也是FFT點數。

另外，這些作法適用於以MATLAB模擬的連續時間的訊號。離散時間訊號的頻譜，從「DFT是對DTFT取樣」的概念發想，應該可以得到類似的作法。

參考資料

Ziemer和Tranter的Principles of Communications。下面這個網站也蠻有意思的，把傅立葉轉換的數學說明得很清楚的樣子，有時間再來慢慢看。

從傅立葉轉換到數位訊號處理