傅立葉變換是物理學家、數學家、工程師和計算機科學家常用的最有用的工具之一。本篇文章我們將使用Python來實現一個連續函數的傅立葉變換。

我們使用以下定義來表示傅立葉變換及其逆變換。

設 f: ℝ → ℂ 是一個既可積又可平方積分的複值函數。那麽它的傅立葉變換，記爲 f̂，是由以下複值函數給出：

同樣地，對于一個複值函數 ĝ，我們定義其逆傅立葉變換（記爲 g）爲

這些積分進行數值計算是可行的，但通常是棘手的——特別是在更高維度上。所以必須采用某種離散化的方法。

在Numpy文檔中關于傅立葉變換如下，實現這一點的關鍵是離散傅立葉變換（DFT）：

當函數及其傅立葉變換都被離散化的對應物所取代時，這被稱爲離散傅立葉變換（DFT）。離散傅立葉變換由于計算它的一種非常快速的算法而成爲數值計算的重要工具，這個算法被稱爲快速傅立葉變換（FFT），這個算法最早由高斯（1805年）發現，我們現在使用的形式是由Cooley和Tukey公開的

根據Numpy文檔，一個具有 n 個元素的序列 a₀, …, aₙ₋₁ 的 DFT 計算如下：

我們將積分分解爲黎曼和。在 n 個不同且均勻間隔的點 xₘ = x₀ + m Δx 處對 x 進行采樣，其中 m 的範圍從 0 到 n-1，x₀ 是任意選擇的最左側點。然後就可以近似表示積分爲

現在對變量 k 進行離散化，在 n 個均勻間隔的點 kₗ = l Δk 處對其進行采樣。然後積分變爲:

這使得我們可以用類似于 DFT 的形式來計算函數的傅立葉變換。這與DFT的計算形式非常相似，這讓我們可以使用FFT算法來高效計算傅立葉變換的近似值。

最後一點是將Δx和Δk聯系起來，以便指數項變爲-2π I ml/n，這是Numpy的實現方法;

這就是不確定性原理，所以我們得到了最終的方程

我們可以對逆變換做同樣的處理。在Numpy中，它被定義爲

1/n是歸一化因子：

概念和公式我們已經通過Numpy的文檔進行了解了，下面開始我們自己的Python實現

import numpy as npimport matplotlib.pyplot as pltdef 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, gdef 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

我們來通過一些例子看看我們自己實現是否正確。

第一個例子：階躍函數

函數在-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的采樣值和整個連續體上計算的解析解:

k, g = fourier_transform_1d(f, x, sort_results=True) # make it easier to plotkk = 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();

看起來是沒問題的，然後我們把它轉換回來:

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();

我們可以清楚地看到不連續邊緣處的 Gibbs 現象——這是傅裏葉變換的一個預期特征。

第二個例子：高斯PDF

傅裏葉變換

下面，我們繪制數值傅裏葉變換和解析值:

以及傅裏葉逆變換與原函數的對比

可以看到，我們的實現沒有任何問題

最後，如果你對機器學習的基礎計算和算法比較感興趣，可以多多關注Numpy和SK-learn的文檔（還有scipy但是這個更複雜），這兩個庫不僅有很多方法的實現，還有這些方法的詳細解釋，這對于我們學習是非常有幫助的。

例如本文的一些數學的公式和概念就是來自于Numpy的文檔，有興趣的可以直接看看

https://avoid.overfit.cn/post/546692942b9144a5a56d734c5a007099

作者：Alessandro Morita Gagliardi