傅里葉變換算法和Python代碼實現(xiàn)
傅立葉變換是物理學家、數(shù)學家、工程師和計算機科學家常用的最有用的工具之一。本篇文章我們將使用Python來實現(xiàn)一個連續(xù)函數(shù)的傅立葉變換。
我們使用以下定義來表示傅立葉變換及其逆變換。
設 f: ? → ? 是一個既可積又可平方積分的復值函數(shù)。那么它的傅立葉變換，記為 f?，是由以下復值函數(shù)給出：
同樣地，對于一個復值函數(shù) g?，我們定義其逆傅立葉變換（記為 g）為
這些積分進行數(shù)值計算是可行的，但通常是棘手的——特別是在更高維度上。所以必須采用某種離散化的方法。
在Numpy文檔中關(guān)于傅立葉變換如下，實現(xiàn)這一點的關(guān)鍵是離散傅立葉變換（DFT）：
當函數(shù)及其傅立葉變換都被離散化的對應物所取代時，這被稱為離散傅立葉變換（DFT）。離散傅立葉變換由于計算它的一種非?？焖俚乃惴ǘ蔀閿?shù)值計算的重要工具，這個算法被稱為快速傅立葉變換（FFT），這個算法最早由高斯（1805年）發(fā)現(xiàn)，我們現(xiàn)在使用的形式是由Cooley和Tukey公開的根據(jù)Numpy文檔，一個具有 n 個元素的序列 a?, …, a??? 的 DFT 計算如下：
我們將積分分解為黎曼和。在 n 個不同且均勻間隔的點 x? = x? + m Δx 處對 x 進行采樣，其中 m 的范圍從 0 到 n-1，x? 是任意選擇的最左側(cè)點。然后就可以近似表示積分為
現(xiàn)在對變量 k 進行離散化，在 n 個均勻間隔的點 k? = l Δk 處對其進行采樣。然后積分變?yōu)?
這使得我們可以用類似于 DFT 的形式來計算函數(shù)的傅立葉變換。這與DFT的計算形式非常相似，這讓我們可以使用FFT算法來高效計算傅立葉變換的近似值。
最后一點是將Δx和Δk聯(lián)系起來，以便指數(shù)項變?yōu)?2π I ml/n，這是Numpy的實現(xiàn)方法;
這就是不確定性原理，所以我們得到了最終的方程
我們可以對逆變換做同樣的處理。在Numpy中，它被定義為
1/n是歸一化因子：
概念和公式我們已經(jīng)通過Numpy的文檔進行了解了，下面開始我們自己的Python實現(xiàn)
import numpy as np
import matplotlib.pyplot as plt
def fourier_transform_1d(func, x, sort_results=False):
"""
Computes the continuous Fourier transform of function `func`, following the physicist's convention
Grid x must be evenly spaced.
Parameters
----------
- func (callable): function of one argument to be Fourier transformed
- x (numpy array) evenly spaced points to sample the function
- sort_results (bool): reorders the final results so that the x-axis vector is sorted in a natural order.
Warning: setting it to True makes the output not transformable back via Inverse Fourier transform
Returns
--------
- k (numpy array): evenly spaced x-axis on Fourier domain. Not sorted from low to high, unless `sort_results` is set to True
- g (numpy array): Fourier transform values calculated at coordinate k
"""
x0, dx = x[0], x[1] - x[0]
f = func(x)
g = np.fft.fft(f) # DFT calculation
# frequency normalization factor is 2*np.pi/dt
w = np.fft.fftfreq(f.size)*2*np.pi/dx
# Multiply by external factor
g *= dx*np.exp(-complex(0,1)*w*x0)
if sort_results:
zipped_lists = zip(w, g)
sorted_pairs = sorted(zipped_lists)
sorted_list1, sorted_list2 = zip(*sorted_pairs)
w = np.array(list(sorted_list1))
g = np.array(list(sorted_list2))
return w, g
def inverse_fourier_transform_1d(func, k, sort_results=False):
"""
Computes the inverse Fourier transform of function `func`, following the physicist's convention
Grid x must be evenly spaced.
Parameters
----------
- func (callable): function of one argument to be inverse Fourier transformed
- k (numpy array) evenly spaced points in Fourier space to sample the function
- sort_results (bool): reorders the final results so that the x-axis vector is sorted in a natural order.
Warning: setting it to True makes the output not transformable back via Fourier transform
Returns
--------
- y (numpy array): evenly spaced x-axis. Not sorted from low to high, unless `sort_results` is set to True
- h (numpy array): inverse Fourier transform values calculated at coordinate x
"""
dk = k[1] - k[0]
f = np.fft.ifft(func) * len(k) * dk /(2*np.pi)
x = np.fft.fftfreq(f.size)*2*np.pi/dk
if sort_results:
zipped_lists = zip(x, f)
sorted_pairs = sorted(zipped_lists)
sorted_list1, sorted_list2 = zip(*sorted_pairs)
x = np.array(list(sorted_list1))
f = np.array(list(sorted_list2))
return x, f我們來通過一些例子看看我們自己實現(xiàn)是否正確。
第一個例子：階躍函數(shù)
函數(shù)在-1/2和1/2之間是1，在其他地方是0。它的傅里葉變換是
N = 2048
# Define the function f(x)
f = lambda x: np.where((x >= -0.5) & (x <= 0.5), 1, 0)
x = np.linspace(-1, 1, N)
plt.plot(x, f(x));
畫出傅里葉變換，以及在k的采樣值和整個連續(xù)體上計算的解析解:
k, g = fourier_transform_1d(f, x, sort_results=True) # make it easier to plot
kk = np.linspace(-30,30, 100)
plt.plot(k, np.real(g), label='Numerical');
plt.plot(k, np.sin(k/2)/(k/2), linestyle='-.', label='Analytic (samples)')
plt.plot(kk, np.sin(kk/2)/(kk/2), linestyle='--', label='Analytic (full)')
plt.xlim(-30, 30)
plt.legend();
看起來是沒問題的，然后我們把它轉(zhuǎn)換回來:
k, g = fourier_transform_1d(f, x)
y, h = inverse_fourier_transform_1d(g, k, sort_results=True)
plt.plot(y, np.real(h), label='Numerical transform')
plt.plot(x, f(x), linestyle='--', label='Analytical')
plt.legend();
我們可以清楚地看到不連續(xù)邊緣處的 Gibbs 現(xiàn)象——這是傅里葉變換的一個預期特征。
第二個例子：高斯PDF
傅里葉變換
下面，我們繪制數(shù)值傅里葉變換和解析值:
以及傅里葉逆變換與原函數(shù)的對比
可以看到，我們的實現(xiàn)沒有任何問題
最后，如果你對機器學習的基礎計算和算法比較感興趣，可以多多關(guān)注Numpy和SK-learn的文檔（還有scipy但是這個更復雜），這兩個庫不僅有很多方法的實現(xiàn)，還有這些方法的詳細解釋，這對于我們學習是非常有幫助的。
例如本文的一些數(shù)學的公式和概念就是來自于Numpy的文檔，有興趣的可以直接看看
https://numpy.org/doc/stable/reference/routines.fft.html
