手動(dòng)實(shí)現(xiàn)離散傅里葉正變換與逆變換(程序+例子)

前言

實(shí)際野外采集到的信號(hào)都是有一定時(shí)間間隔的，因?yàn)楸疚闹饕P(guān)注于"離散傅里葉變換(DFT)"。另外，本文主要關(guān)注的是"如何使用"DFT這樣一個(gè)工具(公式的使用+程序編寫)，并非介紹其背后的原理。最后，本文將自己編寫的程序與Matlab自帶的相關(guān)語(yǔ)句的結(jié)果做對(duì)比，驗(yàn)證方法與程序的正確性。

DFT相關(guān)公式及使用

我們假設(shè)原始時(shí)域信號(hào)為x(n)，經(jīng)DFT后頻域信號(hào)為X(k)，其中n是采樣點(diǎn)的序號(hào)，k是對(duì)應(yīng)到頻率域中的頻率值，N是采樣點(diǎn)總數(shù)。傅里葉正變換公式為：

公式1：DFT原始公式

注意1：在等號(hào)右邊的最后一項(xiàng)中，nk乘積是的冪指數(shù)！
注意2：是一個(gè)常數(shù)！其公式為：

公式1怎么用？先把它展開看看：

公式2：把公式1展開成多項(xiàng)來(lái)看

公式2看上去還是不夠直觀，我們?cè)侔阉镁仃嚤硎境桑?div id="u0z1t8os" class="image-package">

公式3：公式2的矩陣形式

是不是清晰多了呢？
已經(jīng)知道是一個(gè)常數(shù)！
那公式3"等號(hào)右邊第一項(xiàng)"可進(jìn)一步簡(jiǎn)化：

公式4：對(duì)公式3中右邊第一項(xiàng)的簡(jiǎn)化

到此，DFT公式可變?yōu)椋?img class="math-inline" src="https://math.jianshu.com/math?formula=X%20%3D%20W_%7BN%7D%5E%7BE%7D%5C%20x" alt="X = W_{N}^{E}\ x" mathimg="1"> 這樣的矩陣相乘。

有了上面的公式的基礎(chǔ)，我們可以手寫Matlab程序：

data = [2 3 8 9 12 4 8 10];   % 原始離散時(shí)域信號(hào)
N = length(data);

% DFT需要的相關(guān)矩陣:
WN = exp(-i*2*pi/N);  % 常數(shù)
WN_nk = zeros(N)+WN;  % WN_kn矩陣(初始)
xk = data';     % 時(shí)域信號(hào)振幅(列矩陣)
E = zeros(N);   % 輔助的E(WN_kn的冪,單獨(dú)拿出來(lái)算)

%%% 傅里葉正變換即結(jié)果 %%%
for row = 0:N-1
    for cow = 0:N-1
        E(row+1,cow+1) = row*cow;
        WN_nk(row+1,cow+1) = WN_nk(row+1,cow+1)^E(row+1,cow+1);
    end
end

Xk = WN_nk * xk;   % 頻域結(jié)果
fft(data);         % Matlab自帶語(yǔ)句,兩者結(jié)果對(duì)比一樣

IDFT相關(guān)公式及使用

很DFT的公式非常相似，只不過(guò)多了負(fù)號(hào)：

公式5：IDFT原始公式

將公式轉(zhuǎn)換成矩陣，再寫出矩陣E和對(duì)應(yīng)的即可。最后不要忘記除以總采樣點(diǎn)數(shù)N。

對(duì)應(yīng)的手寫Matlab程序?yàn)椋?/p>

% 繼DFT代碼后跟著寫
% IDFT需要的相關(guān)矩陣：
WN_nk_n = zeros(N)+WN;  % WN_kn矩陣(初始)
E_n = zeros(N);                   % 輔助的E(WN_kn的冪,單獨(dú)拿出來(lái)算)

for row = 0:N-1
    for cow = 0:N-1
        E_n(row+1,cow+1) = -row*cow;     % 注意負(fù)號(hào)
        WN_nk_n(row+1,cow+1) = WN_nk_n(row+1,cow+1)^E_n(row+1,cow+1);
    end
end

xk_n = real((WN_nk_n * Xk)/N)'     % IDFT手寫結(jié)果
ifft(Xk)                           % Matlab自帶語(yǔ)句，上下倆結(jié)果對(duì)比一樣。

更多完整程序與相關(guān)實(shí)例程序參考：傅里葉變換程序及實(shí)例.

補(bǔ)充說(shuō)明：

頻率域的優(yōu)越性：時(shí)間域信號(hào)雖然可以反映原始信號(hào)不少的信息，但是這種圖像不易修改！因?yàn)榈厍蛭锢硪巴饨邮盏降男盘?hào)，不論是有效信號(hào)還是特定噪聲，它們都各有自己所屬的頻帶寬度！也就是說(shuō)有效信號(hào)與噪聲它們?cè)?strong>頻率域是很好識(shí)別和區(qū)分開來(lái)的！雖然頻率信息也能少量反應(yīng)在時(shí)域圖像中(圖像的緊密或舒張程度)，但是要想從時(shí)域圖像中將二者分離開來(lái)是非常非常困難的！這也體現(xiàn)出頻域的優(yōu)越性。不只是信號(hào)提取與消除，頻率域或者類頻率域(后面的小波域、τp域等)還能做很多更高級(jí)的操作，這些都是時(shí)域圖像所不能進(jìn)行的?？偨Y(jié)之：時(shí)域信號(hào)好采集不好處理，頻域信號(hào)好處理但直接采集不了。