「離散傅立葉轉換」（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)，這些弦波稱為「基本弦波」。
範例 1輸入（fft/fft01.m）：
% 此範例展示一個簡單正弦波的傅立葉轉換，以雙邊頻譜來顯示
% 此正弦波的頻率恰巧是 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)’);
由上圖可以看出：
- 能量頻譜只有在兩點不為零，其他各點理論上應該是零，但由於計算精度之故，實際值很接近於零，但不一定是零。
- 相位頻譜沒有一定規則，像是亂數。這是由於大部分的能量頻譜很小，因此實際的相位頻譜並沒有任何意義。
如果 x[n] 的頻率不是這些弦波中的其中一個，那麼 DFT 還是會將此訊號拆解成基本弦波的組合，因此能量頻譜就會分散在各個基本弦波，如下所示：
範例 2輸入（fft/fft02.m）：
% 此範例展示一個簡單正弦波的傅立葉轉換，以雙邊頻譜來顯示
% 此正弦波的頻率不是 freqStep 的整數倍，所以雙邊頻譜會「散開」(Smearing)
N = 256; % 點數
fs = 8000; % 取樣頻率
freqStep = fs/N; % 頻域的頻率的解析度
f = 10.5*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
subplot(3,1,2);
freq = freqStep*(-N/2:N/2-1); % 頻域的頻率刻度
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)’);
由於對於實數的 x[n] 而言，我們可以只看單邊頻譜，可用 MATLAB 來進行單邊頻譜的顯示：
範例 3輸入（fft/fft03.m）：
% 同 fft02.m，但以單邊頻譜來顯示
N = 256; % 點數
fs = 8000; % 取樣頻率
freqStep = fs/N; % 頻域的頻率的解析度
f = 10.5*freqStep; % 正弦波的頻率，不是 freqStep 的整數倍
time = (0:N-1)/fs; % 時域的時間刻度
y = cos(2*pi*f*time); % Signal to analyze
Y = fft(y); % Spectrum
freq = freqStep*(0:N/2); % 頻域的頻率刻度
Y = Y(1:length(freq));
% Plot time data
subplot(3,1,1);
plot(time, y, ‘.-‘);
title(‘Sinusoidal signals’);
xlabel(‘Time (seconds)’); ylabel(‘Amplitude’);
axis tight
% Plot spectral magnitude
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)’);
當 k 越大時，代表對應的弦波是高頻訊號，因此我們可以捨棄高頻訊號，只用幾個低頻弦波，來合成原始訊號，達到下列兩項目的：
- 資料壓縮
- 濾除高頻訊號
以下這個範例，顯示如何以弦波來合成一個方波：
範例 4輸入（fft/fftApproximate01.m）：
% This example demos the effect of square wave approximation by DFT
figure
L = 15; N = 25;
x = [ones(1,L), zeros(1,N-L)];
frameSize=length(x);
runNum=3;
for i=1:runNum,
pointNum=ceil(frameSize/(2*runNum)*i); % Actually 2*pointNum-1 coefs are taken
X = fft(x);
magX = abs(X);
remainIndex=[1:pointNum, frameSize-pointNum+2:frameSize];
X2=0*X;
X2(remainIndex)=X(remainIndex);
x2=ifft(X2);
x2=real(x2);
subplot(3,2,2*i-1);
stem(x);
hold on
plot(x2, ‘r’);
hold off
title(sprintf(‘x[n] and %d-points approximation’, 2*pointNum-1));
axis([-inf,inf,-0.5,1.5])
subplot(3,2,2*i);
shiftedMagX=fftshift(magX);
plot(shiftedMagX, ‘.-‘); axis tight
title(‘DFS of x[n]’)
hold on
temp=ifftshift(1:frameSize);
ind=temp(remainIndex);
plot(ind, shiftedMagX(ind), ‘or’); grid on
hold off
end
由此可見當我們所用的 DFT 係數越多，所合成的波形就會越接近原來的波形。
範例 5輸入（fft/fftApproximate02.m）：
% This example demos the effect of FFT approximation
[y, fs]=wavread(‘welcome.wav’);
x=y(2047:2047+237-1);
figure
frameSize=length(x);
runNum=3;
for i=1:runNum,
pointNum=ceil(frameSize/(8*runNum)*i); % Actually 2*pointNum-1 coefs are taken
X = fft(x);
magX = abs(X);
remainIndex=[1:pointNum, frameSize-pointNum+2:frameSize];
X2=0*X;
X2(remainIndex)=X(remainIndex);
x2=ifft(X2);
x2=real(x2);
subplot(3,2,2*i-1);
plot(x, ‘.-‘);
hold on
plot(x2, ‘r’);
hold off
title(sprintf(‘x[n] and %d-points approximation’, 2*pointNum-1));
set(gca, ‘xlim’, [-inf inf]);
subplot(3,2,2*i);
shiftedMagX=fftshift(magX);
plot(shiftedMagX, ‘.-‘);
title(‘DFS of x[n]’)
hold on
temp=ifftshift(1:frameSize);
ind=temp(remainIndex);
plot(ind, shiftedMagX(ind), ‘or’); grid on
hold off
set(gca, ‘xlim’, [-inf inf]);
end
由此可見在實際的聲音波形中，我們只要少數幾個低頻的係數，就可以合成類似原來的訊號，由此也可以知道很多高頻訊號是沒有作用的。
範例 6輸入（fft/fftRepeat01.m）：
% This example demos the effect of FFT approximation
[y, fs]=wavread(‘welcome.wav’);
x=y(2047:2126);
runNum=5;
for i=1:runNum,
repeatedX = x*ones(1,i);
repeatedX = repeatedX(:);
X = fft(repeatedX);
magX = abs(X);
frameSize=length(repeatedX);
subplot(runNum,2,2*i-1);
plot(repeatedX, ‘.-‘); grid on
title(‘x[n]’);
set(gca, ‘xlim’, [-inf inf]);
subplot(runNum,2,2*i);
freq = (0:frameSize/2)*fs/frameSize;
magX = magX(1:length(freq));
plot(freq, magX, ‘.-‘); grid on
title(‘DFS of x[n]’)
axis tight;
end
為什麼 DFT 會插入零點呢？若不使用詳細的數學分析，各位同學是否可以以直覺的方式來推導出這個結果呢？
範例 7輸入（fft/fftZeroPadding01.m）：
% This example demos the effect of zero-padding of DFT
figure
for i=1:3,
L = 5; N = 20*i;
x = [ones(1,L), zeros(1,N-L)];
subplot(3,3,i*3-2);
stem(x);
title(sprintf(‘x[n] with N=%d’,N));
set(gca, ‘xlim’, [-inf inf]);
omega=((1:N)-ceil((N+1)/2))*(2*pi/N);
X = fft(x);
magX = fftshift(abs(X));
subplot(3,3,i*3-1);
plot(omega, magX, ‘.-‘);
title(‘Magnitude of DFT of x[n]’)
set(gca, ‘xlim’, [-inf inf]);
phase=fftshift(angle(X));
subplot(3,3,i*3);
plot(omega, phase, ‘.-‘);
title(‘Phase of DFT of x[n]’)
set(gca, ‘xlim’, [-inf inf]);
set(gca, ‘ylim’, [-pi pi]);
end
換句話說，加上零點，在實質上並沒有增加訊號的資訊，但是 DFT 的點數必須隨之增加，因此在頻譜的效應則是進行內差以增加點數。
提示
範例 8輸入（fft/fftReSample01.m）：
% This example demos the effect of FFT approximation
[y, fs]=wavread(‘welcome.wav’);
x=y(2047:2126);
x=y(2047:2326);
n=length(x);
F = (0:n/2)*fs/n;
runNum=5;
for i=1:runNum,
newX=x(1:2^(i-1):length(x));
newFs=fs/(2^(i-1));
X = fft(newX);
magX = abs(X);
frameSize=length(newX);
subplot(runNum,2,2*i-1);
plot(newX, ‘.-‘);
title(‘x[n]’);
set(gca, ‘xlim’, [-inf inf]);
subplot(runNum,2,2*i);
freq = (0:frameSize/2)*newFs/frameSize;
magX = magX(1:length(freq));
M=nan*F;
M(1:length(magX))=magX;
plot(F, M, ‘.-‘);
title(‘DFS of x[n]’)
axis tight;
end
當我們 Downsample 進行越多次，曲線會越平滑，代表高頻部分已經慢慢不見了！