来源：Deephub Imba
        
       
      
     
    
   
   

    

     

      

       

        
本文约2500字，建议阅读10分钟
        
本文我们将使用Python来实现一个连续函数的傅立叶变换。
       
      
     
    
   

傅立叶变换是物理学家、数学家、工程师和计算机科学家常用的最有用的工具之一。

 当函数及其傅立叶变换都被离散化的对应物所取代时，这被称为离散傅立叶变换（DFT）。离散傅立叶变换由于计算它的一种非常快速的算法而成为数值计算的重要工具，这个算法被称为快速傅立叶变换（FFT），这个算法最早由高斯（1805年）发现，我们现在使用的形式是由Cooley和Tukey公开的

   

    

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

    
    
 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
   


   

    我们来通过一些例子看看我们自己实现是否正确。
    

   

   

    

     

     

     

     

     

     

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