Old Chinese version
在前一節中，我們可以使用「離散時間傅立葉轉換」（簡稱
DTFT）來將一段數位訊號轉換成各個頻譜的分量，但這是一個連續的函數，並不適合在電腦中處理，因此本節將介紹「離散傅立葉轉換」（Discrete Fourier
Transform），簡稱 DFT，其功能是將一段數位訊號轉換成其各個離散頻率的弦波分量，以便後續使用電腦進行各種處理。
如果我們的訊號可以表示成 x[n], n = 0~N-1，那麼 DFT 的公式如下：
X[k]=(1/N)*S n=0 N-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]=S k=0 N-1
X[k]*exp(j*2p *n*k/N), n=0, ..., N-1
Hint
DTFT 和 DFT
很類似，兩者都用來處理離散時間的訊號，只是前者產生一個角頻率的連續函數，而後者產生離散頻率的訊號，更適合使用電腦來進行處理。
這邊有幾點要說明：
如果原始訊號 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]p *n*1/N) +
X[N-1]*exp(j*2p *n*(N-1)/N)p *n*2/N) + X[N-2]*exp(j*2p *n*(N-2)/N)p *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)p *n*k/N) + X[N-k]*exp(j*2p *n)*exp(-j*2p *n*k/N)p *n*k/N) + X[N-k]*exp(-j*2p *n*k/N)p *n*k/N))
如果我們以 mk 表示 X[k] 的大小，以
pk 表示 X[k] 的相位，那麼上式可以化簡如下：
2*Re(mk *exp(j*pk )*exp(j*2p *n*k/N))k *exp(j*(2p *n*k/N + pk ))k *cos(2p *n*k/N + pk )
一般而言，N 是 2 的倍數，因此 x[n]
可以表示成
x[n] = X[0] + 2*S k=1 N/2-1 mk *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。
首先，我們使用 MATLAB 的 fft 指令來驗證實數訊號之 DFT 係數的共軛性，如下：
Example 1: fftSym01.m <![CDATA[% This example demonstrates the pair-wise conjugate of DFT (此範例展示實數訊號之 DFT 係數的共軛性)
N=64; % Length of vector
x=randn(N, 1);
z=fft(x);
plot(z, 'o'); grid on
%compass(z);
% Connect conjugate pairs (將上下對稱的共軛點連起來)
for i=2:N/2+1
twoPoint=z([i, N-i+2]);
line(real(twoPoint), imag(twoPoint), 'color', 'r');
end]]>
Hint
對於實數訊號 x[n] 而言：
當 N 為偶數時，只有 X[0] 和 X[N/2] 是單一實數，其餘均為成對出現之共軛複數。
當 N 為奇數時，只有 X[0] 是單一實數，其餘均為成對出現之共軛複數。
如果 x[n] 恰巧是這些基本弦波中的其中一個，那麼計算出來的雙邊頻譜，應該只有兩個係數不為零，我們可用 MATLAB 驗證如下：
Example 2: fft01.m <![CDATA[% This example demonstrates the two-side DFT of a sinusoidal function
% (此範例展示一個簡單正弦波的傅立葉轉換，以雙邊頻譜來顯示)
% Since the sinusoidal function has a frequency to be a multiple of fs/N, the two-side DFT have only two nonzero terms.
% (此正弦波的頻率恰巧是 freqStep 的整數倍，所以雙邊頻譜應該只有兩個非零點)
N = 256; % length of vector (點數)
fs = 8000; % sample rate (取樣頻率)
freqStep = fs/N; % freq resolution in spectrum (頻域的頻率的解析度)
f = 10*freqStep; % freq of the sinusoid (正弦波的頻率，恰是 freqStep 的整數倍)
time = (0:N-1)/fs; % time resolution in time-domain (時域的時間刻度)
y = cos(2*pi*f*time); % signal to analyze
Y = fft(y); % spectrum
Y = fftshift(Y); % put zero freq at the center (將頻率軸的零點置中)
% 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); % freq resolution in spectrum (頻域的頻率的解析度)
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 還是會將此訊號拆解成基本弦波的組合，因此能量頻譜就會分散在各個基本弦波，如下所示：
Example 3: fft02.m <![CDATA[% This example demonstrates the one-side DFT of a sinusoidal function (此範例展示一個簡單正弦波的傅立葉轉換，以雙邊頻譜來顯示)
% Since the sinusoidal function has a frequency not a multiple of fs/N, the two-side DFT smears. (此正弦波的頻率不是 freqStep 的整數倍，所以雙邊頻譜會「散開」(Smearing))
N = 256; % length of vector (點數)
fs = 8000; % sample rate (取樣頻率)
freqStep = fs/N; % freq resolution in spectrum (頻域的頻率的解析度)
f = 10.5*freqStep; % freq of the sinusoid (正弦波的頻率，不是 freqStep 的整數倍)
time = (0:N-1)/fs; % time resolution in time-domain (時域的時間刻度)
signal = cos(2*pi*f*time); % signal to analyze
[mag, phase, freq]=fftTwoSide(signal, fs, 1); % compute and plot the two-side DFT]]>
Example 4: fft03.m <![CDATA[% Same as fft02.m but use one-side DFT instead (同 fft02.m，但以單邊頻譜來顯示)
N = 256; % length of vector (點數)
fs = 8000; % sample rate (取樣頻率)
freqStep = fs/N; % freq resolution in spectrum (頻域的頻率的解析度)
f = 10.5*freqStep; % freq of the sinusoid (正弦波的頻率，不是 freqStep 的整數倍)
time = (0:N-1)/fs; % time resolution in time-domain (時域的時間刻度)
signal = cos(2*pi*f*time); % signal to analyze
[mag, phase, freq]=fftOneSide(signal, fs, 1); % Compute and plot one-side DFT]]>
Hint
一般實際世界中的訊號，都是實數訊號，因此我們只要使用 fftOneSide.m
來計算單邊頻譜即可。
下面這個範例，顯示一個語音音框的單邊頻譜：
Example 5: fft04.m <![CDATA[% This example demonstrates the DFT of a real-world audio signal (顯示一個語音音框的單邊頻譜)
[y, fs]=wavread('welcome.wav');
signal=y(2047:2047+237-1);
[mag, phase, freq]=fftOneSide(signal, fs, 1);]]>
Example 6: fftApproximate01.m <![CDATA[% 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('DFT of x[n]')
hold on
temp=ifftshift(1:frameSize);
ind=temp(remainIndex);
plot(ind, shiftedMagX(ind), 'or'); grid on
hold off
end]]>
Example 7: fftApproximate02.m <![CDATA[% 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('DFT 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
]]>
Example 8: fftRepeat01.m <![CDATA[% This example demos the effect of FFT for purely periodic signals
[y, fs]=wavread('welcome.wav');
x=y(2047:2126); % A full fundamental period
runNum=5;
for i=1:runNum
repeatedX = x*ones(1,i);
signal = repeatedX(:)+0.1;
% signal=zeros(runNum*length(x), 1); % Zero-padding version
% signal(1:length(repeatedX))=repeatedX(:); % Zero-padding version
[mag, phase, freq, powerDb]=fftOneSide(signal, fs);
mag=mag/length(signal); % Divided by vector length to normalize magnitude (due to the formula used by MATLAB)
subplot(runNum,2,2*i-1);
plot(signal, '.-'); grid on
title('x[n]'); set(gca, 'xlim', [-inf inf]);
subplot(runNum,2,2*i);
plot(freq, mag, '.-'); grid on;
% set(gca, 'yscale', 'log');
title('DFT of x[n]'); axis tight;
end]]>
Example 9: fftZeroPadding01.m <![CDATA[% This example demos the effect of zero-padding of DFT
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]]>
Hint
一般在進行音訊處理時，如果音框的點數不是 2n 的格式，我們通常就是使用補零的方式，使最後的點數變成
2n ，這樣在進行 DFT 時，才比較方便採取比較快速的 FFT 計算方式。
若在時域進行 Downsample，則在頻域所產生的效果如下：
Example 10: fftReSample01.m <![CDATA[% 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('DFT of x[n]')
axis tight;
end
]]>
Example 11: fftViaLse01.m <![CDATA[% FFT via least-squares estimate
N=8;
fs=1;
time=(0:N-1)'/fs;
x=rand(N,1)*2-1;
A=ones(N,1);
for i=1:N/2
A=[A, cos(2*pi*(i*fs/N)*time), sin(2*pi*(i*fs/N)*time)];
end
th=A\x;
plotNum=fix(N/2)+2;
subplot(plotNum, 1, 1);
N2=(N-1)*5+1; % Insert 4 points between every 2 points for better observation (兩點間插入四點，以便觀察波形)
timeFine=linspace(min(time), max(time), N2);
x2=th(1)*ones(N,1);
plot(timeFine, th(1)*ones(N2,1), '.-', time, x2, 'or'); % Plot the first term (畫出第一項)
ylabel('Term 0 (DC)'); axis([0 N/fs -1 1]); grid on
for i=1:N/2 % Plot terms 2 to 1+N/2 (畫出第二至第 1+N/2 項)
freq=i*fs/N;
y=th(2*i)*cos(2*pi*freq*time)+th(2*i+1)*sin(2*pi*freq*time); % a term (每一項)
x2=x2+y;
fprintf('i=%d, sse=%f\n', i, norm(x-x2)/sqrt(N));
subplot(plotNum, 1, i+1);
yFine=th(2*i)*cos(2*pi*(i*fs/N)*timeFine)+th(2*i+1)*sin(2*pi*(i*fs/N)*timeFine); % Finer verison for plotting
plot(timeFine, yFine, '.-', time, y, 'or'); ylabel(sprintf('Term %d', i));
axis([0 N/fs -1 1]); grid on
end
% Plot the original signal (畫出原來訊號)
subplot(plotNum, 1, plotNum)
plot(time, x, 'o-'); axis([0 N/fs -1 1]); grid on
xlabel('Time'); ylabel('Orig signals');
% Transform LSE result back to fft format for comparison (將 th 轉回 fft 並比較結果)
F=fft(x);
F2=[];
F2(1)=th(1)*N;
for i=1:N/2
F2(i+1)=(th(2*i)-sqrt(-1)*th(2*i+1))*N/2;
if (i==N/2)
F2(i+1)=2*F2(i+1);
end
end
% symmetric of DFT (DFT 的對稱性)
for i=N/2+2:N
F2(i)=F2(N-i+2)';
end
error1=sum(abs(F2-F.')); % F.' is simple transpose (F.' 是不進行共軛轉換的轉置)
error2=sum(abs(F2-F')); % F' is conjugate transpose (F' 是進行共軛轉換的轉置)
fprintf('Errors after transforming LSE to DFT coefficients (將 LSE 轉換成 DFT 係數的誤差)：error1=%f, error2=%f\n', error1, error2);
fprintf('Due to the symmetry of DFT, one of the above error terms should be zero. (由於 DFT 的對稱性，上述誤差應該有一項為零。)\n');]]> <![CDATA[i=1, sse=0.407551
i=2, sse=0.396878
i=3, sse=0.147817
i=4, sse=0.000000
Errors after transforming LSE to DFT coefficients (將 LSE 轉換成 DFT 係數的誤差)：error1=0.000000, error2=9.412318
Due to the symmetry of DFT, one of the above error terms should be zero. (由於 DFT 的對稱性，上述誤差應該有一項為零。)
]]>
Audio Signal Processing and Recognition (音訊處理與辨識)
No comments:
Post a Comment