傅里葉變換算法和Python代碼實現

作者：Alessandro Morita Ga 2024-03-12 12:49:17

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

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

設 f: ? → ? 是一個既可積又可平方積分的復值函數。那么它的傅立葉變換，記為 f?，是由以下復值函數給出：

同樣地，對于一個復值函數 g?，我們定義其逆傅立葉變換（記為 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 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

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

第一個例子：階躍函數

函數在-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 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();

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

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