傅里葉變換算法和Python代碼實現(xiàn)

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

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

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

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

同樣地，對于一個復(fù)值函數(shù) g?，我們定義其逆傅立葉變換（記為 g）為

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

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

當(dāng)函數(shù)及其傅立葉變換都被離散化的對應(yīng)物所取代時，這被稱為離散傅立葉變換（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();