「離散傅立葉轉換」（Discrete Fourier Transform）簡稱 DFT，其功能是將一段數位訊號轉換成其各個頻率的正弦波分量。如果我們的訊號可以表示成 x[n], n = 0~N-1，那麼 DFT 的公式如下：
X[k]=(1/N)*Sn=0N-1 x[n]*exp(-j*2p*n*k/N), k=0, ..., N-1 這些傅立葉係數 X[k] 所代表的資訊是 k 的函數，而 k 直接和頻率有正比關係，因此這些係數 X[k] 通稱為「頻譜」（Spectrum），而對於 X[k] 的分析，我們通稱為「頻譜分析」（Spectral Analysis）。我們也可以由這些傅立葉係數 X[k]，來反推原始訊號 x[n]，如下： x[n]=Sk=0N-1 X[k]*exp(j*2p*n*k/N), n=0, ..., N-1這邊有幾點要說明：
- 如果原始訊號 x[n] 有 N 點，那麼轉換出來的訊號 X[k] 也會有 N 點。
- 一般而言，X[k] 是一個複數，其大小是 |(X[k])| （abs(X[k]) in MATLAB），相位是 ∠X[k] (angle(X[k]) or atan(imag(X[k])/real(X[k]) in MATLAB)。
- 如果原始訊號 x[n] 都是實數，那麼 X[k] 和 X[N-k] 會是共軛複數，滿足 |(X[k])| = |(X[N-k])| 以及 ∠X[k] = -∠X[N-k]。
如果 x[n] 是實數，我們可以將 x[n] 表示如下： x[n] = X[0]
+ X[1]*exp(j*2p*n*1/N) + X[N-1]*exp(j*2p*n*(N-1)/N)
+ X[2]*exp(j*2p*n*2/N) + X[N-2]*exp(j*2p*n*(N-2)/N)
+ X[3]*exp(j*2p*n*3/N) + X[N-3]*exp(j*2p*n*(N-3)/N)
+ ... 對上述第 k 項而言，我們有 X[k]*exp(j*2p*n*k/N) + X[N-k]*exp(j*2p*n*(N-k)/N)
= X[k]*exp(j*2p*n*k/N) + X[N-k]*exp(j*2p*n)*exp(-j*2p*n*k/N)
= X[k]*exp(j*2p*n*k/N) + X[N-k]*exp(-j*2p*n*k/N)
= 2*Re(X[k]*exp(j*2p*n*k/N))
如果我們以 mk 表示 X[k] 的大小，以 pk 表示 X[k] 的相位，那麼上式可以化簡如下：
2*Re(mk*exp(j*pk)*exp(j*2p*n*k/N))
= 2*Re(mk*exp(j*(2p*n*k/N + pk))
= 2*mk*cos(2p*n*k/N + pk)
一般而言，N 是 2 的倍數，因此 x[n] 可以表示成 x[n] = X[0] + 2*Sk=0N/2mk*cos(2p*n*k/N + pk) + mN/2*cos(p*n + pN/2) 換句話說，我們可以將原始訊號拆解成一個直流訊號 X[0] 再加上 N/2 個弦波的組合，這些弦波的震幅就是 X[k] 的大小，而相位則是 X[k] 的相位。因此對於實數的 x[n] 而言，我們只需要看單邊的 X[k]，k = 0 ~ N/2，此種可稱為「單邊頻譜」。組成這些單邊頻譜的弦波共有 1+N/2 個，頻率由小到大分別是 fs/N*(0:N/2)，這些弦波稱為「基本弦波」。
若是套用上述公式來計算 DFT，所需要的複雜度是 O(n2)，但在 1965 年，有兩位學者提出來一套更精簡的演算法，所需的複雜度只有 O(n log n)，這一套演算法稱為「快速傅立葉轉換」（Fast Fourier Transform，簡稱 FFT），換句話說，FFT 是用來計算 DFT 的快速方法。若使用 MATLAB，相關的指令也是 fft。
如果 x[n] 恰巧是這些基本弦波中的其中一個，那麼計算出來的雙邊頻譜，應該只有兩個係數不為零，我們可用 MATLAB 驗證如下：
% 此範例展示一個簡單正弦波的傅立葉轉換，以雙邊頻譜來顯示 % 此正弦波的頻率恰巧是 freqStep 的整數倍，所以雙邊頻譜應該只有兩個非零點
N = 256; % 點數
fs = 8000; % 取樣頻率
freqStep = fs/N; % 頻域的頻率的解析度
f = 10*freqStep; % 正弦波的頻率，恰是 freqStep 的整數倍
time = (0:N-1)/fs; % 時域的時間刻度
y = cos(2*pi*f*time); % Signal to analyze
Y = fft(y); % Spectrum
Y = fftshift(Y); % 將頻率軸的零點置中
% Plot time data subplot(3,1,1); plot(time, y, '.-');
title('Sinusoidal signals');
xlabel('Time (seconds)');
ylabel('Amplitude');
axis tight % Plot spectral magnitude
freq = freqStep*(-N/2:N/2-1); % 頻域的頻率刻度
subplot(3,1,2);
plot(freq, abs(Y), '.-b');
grid on xlabel('Frequency)');
ylabel('Magnitude (Linear)'); % Plot phase
subplot(3,1,3);
plot(freq, angle(Y), '.-b');
grid on xlabel('Frequency)');
ylabel('Phase (Radian)');
請先 登入 以發表留言。