Simulate-Fourier Transform

上一篇有提到關於Discrete Fourier Transform和Continuous Fourier Transform的差別，在這裡要做一些MATLAB的模擬

首先先來回憶一下Discrete Fourier Transform的公式\[\chi(m\Delta f)=\sum_{n=0}^{N-1}\chi(n\Delta t)e^{-i\frac{2\pi mn}{N}\Delta t},m=0,1,\ldots,N-1\]

其中\[\Delta f=\frac{1}{T}\] \[\Delta t=\frac{T}{N}\]

觀察上面的公式，我在MATLAB做了一點小模擬

% sampling delta t
t = 0:.001:1;

% origin signal
x = sin(2*pi*90*t) + sin(2*pi*400*t);

% add noise to signal
y = x + 2*randn(size(t));

figure;
plot(y);
xlabel('Time');
% discrete fourier transform
Y = fft(y);
Pyy = Y.*conj(Y);

figure;
plot(Pyy);
title('Power spectral density');
xlabel('Frequency');

利用上面的MATLAB code可以得到下面兩張圖

上圖為在Time domain的Wave長相，很不幸的，我們無法從Time domain得到什麼資訊

上圖為Frequency Domain的頻譜圖，可以發現我們的波是由90Hz和400Hz所組成的，而比較高頻的部分為對稱的部分。

接下來就是有趣的地方，如果我們增加取樣時間為3秒

% sampling delta t
t = 0:.001:3;

% origin signal
x = sin(2*pi*90*t) + sin(2*pi*400*t);

% add noise to signal
y = x + 2*randn(size(t));

figure;
plot(y);
xlabel('Time');
% discrete fourier transform
Y = fft(y);
Pyy = Y.*conj(Y);

figure;
plot(Pyy);
title('Power spectral density');
xlabel('Frequency');

可以得到如下圖的結果

好了我們發現怎麼90跑到270的位置，400跑到1200的位置喏？不要忘了$\Delta f$代表的是每$\Delta f$為一個值，也就是說上圖的270代表的是index，也就是第270個點的值，那到底真正的訊號是多少Hz呢？很簡單，$\Delta f$代表index和index之間的間隔，所以代表我們每$\Delta f$計算一個值，所以真正的頻率為$index*\Delta f$啦！也就是說270為index乘上$\Delta f$為$\frac{1}{3}$就得到真正的頻率90Hz，所以頻率仍然為90Hz和400Hz啦！所以不要忘記一件是那就是把Time domain的signal丟入fft計算之後出來的為index，需要再乘以$\Delta f$才能夠得到真正的frequency！

如果改變time的區間呢？

% sampling delta t
t = 5:.001:6;

% origin signal
x = sin(2*pi*90*t) + sin(2*pi*400*t);

% add noise to signal
y = x + 2*randn(size(t));

figure;
plot(y);
xlabel('Time');
% discrete fourier transform
Y = fft(y);
Pyy = Y.*conj(Y);

figure;
plot(Pyy);
title('Power spectral density');
xlabel('Frequency');

可以發現到我們得到的frequency仍然為90Hz和400Hz

如果改變sampling rate呢？

% sampling delta t
t = 0:.01:1;

% origin signal
x = sin(2*pi*90*t) + sin(2*pi*400*t);

% add noise to signal
y = x + 2*randn(size(t));

figure;
plot(y);
xlabel('Time');
% discrete fourier transform
Y = fft(y);
Pyy = Y.*conj(Y);

figure;
plot(Pyy);
title('Power spectral density');
xlabel('Frequency');

可以感覺到Time domain的樣子好像不太對，似乎frequency變慢的感覺

從frequency來看，的確frequency被誤判為10Hz喏，且會發現400Hz好像不見喏！這種現象就是Aliasing

後記:當然更多MATLAB關於FFT的information可以去官網看看，花了不少時間去看Discrete Fourier Transform總算有點了解，很多人說Discrete Fourier Transform為Continuity Fourier Transform演變而來的，只要Continuity通喏Discrete也跟著通喏，個人認為不然，Discrete Fourier Transform牽扯到了不少關於Sampling的知識，雖然沒有很難，但是有不少小細節是與Continuity是不一樣的，這裡需要花點時間來思考才能釐清兩者的不同！