python - python 中离散傅立叶变换 (FT) 的教程、技巧和香蕉皮

Question

我对在Python中进行傅立叶变换时有用的技巧和香蕉皮感兴趣。

请提供介绍性代码示例以进入该主题，并就高级主题提供更多建议，例如:频率滤波器、连续 FT、高维 FT。

Answer 1

Python 中的二维离散傅立叶变换入门(numpy):

1。频率排序

应该小心 numpy.fft 库的频率排序。fft 例程首先对零进行排序，然后是正频率，最后是负频率(图 1b、1c)。函数 fftshift 可用于建立预期的排序，即从负值到正值需要 fftshift 例程(图 1d、1e)。

代码示例:

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2,3, figsize=(18.3,12.5))


#Spectrum and FT
spec=np.ndarray(shape=(12,9))
for ii in range(0,spec.shape[0]):
    for jj in range(0,spec.shape[1]):
        spec[ii,jj]=np.cos(2*np.pi*ii/spec.shape[0])*np.cos(4*np.pi*jj/spec.shape[1])
        
spec_ft = np.fft.fft2(spec)
spec_ftShifted = np.fft.fftshift(np.fft.fft2(spec))
   
#Frequencies / labels axes
xFreq=np.around(np.fft.fftfreq(spec.shape[0]),decimals=2)
yFreq=np.around(np.fft.fftfreq(spec.shape[1]),decimals=2)
xFreqShifted=np.around(np.fft.fftshift(np.fft.fftfreq(spec.shape[0])),decimals=2)
yFreqShifted=np.around(np.fft.fftshift(np.fft.fftfreq(spec.shape[1])),decimals=2)


#Plotting
ax=ax1
ax.set_title('1a) Spectrum')
im=ax.imshow(spec.T,origin='lower',cmap='bwr',interpolation='none')
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax2
ax.set_title('1b) FT Spectrum (not shifted): real part')
im=ax.imshow(spec_ft.real.T,origin='lower',cmap='afmhot',interpolation='none')
ax.set_xticks([0,1,int(spec.shape[0]/2),int(spec.shape[0]/2+1),spec.shape[0]-1])
ax.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(spec.shape[0]/2)],xFreq[int(spec.shape[0]/2+1)],xFreq[spec.shape[0]-1]])
ax.set_yticks([0,1,int(spec.shape[1]/2),int(spec.shape[1]/2+1),spec.shape[1]-1])
ax.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(spec.shape[1]/2)],yFreq[int(spec.shape[1]/2+1)],yFreq[spec.shape[1]-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax3
ax.set_title('1c) FT Spectrum (not shifted): imag part')
im=ax.imshow(spec_ft.imag.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,1,int(spec.shape[0]/2),int(spec.shape[0]/2+1),spec.shape[0]-1])
ax.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(spec.shape[0]/2)],xFreq[int(spec.shape[0]/2+1)],xFreq[spec.shape[0]-1]])
ax.set_yticks([0,1,int(spec.shape[1]/2),int(spec.shape[1]/2+1),spec.shape[1]-1])
ax.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(spec.shape[1]/2)],yFreq[int(spec.shape[1]/2+1)],yFreq[spec.shape[1]-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax5
ax.set_title('1d) FT Spectrum (shifted): real part')
im=ax.imshow(spec_ftShifted.real.T,origin='lower',cmap='afmhot',interpolation='none')
ax.set_xticks([0,1,int(spec.shape[0]/2),int(spec.shape[0]/2+1),spec.shape[0]-1])
ax.set_xticklabels([xFreqShifted[0],xFreqShifted[1],xFreqShifted[int(spec.shape[0]/2)],xFreqShifted[int(spec.shape[0]/2+1)],xFreqShifted[spec.shape[0]-1]])
ax.set_yticks([0,1,int(spec.shape[1]/2),int(spec.shape[1]/2+1),spec.shape[1]-1])
ax.set_yticklabels([yFreqShifted[0],yFreqShifted[1],yFreqShifted[int(spec.shape[1]/2)],yFreqShifted[int(spec.shape[1]/2+1)],yFreqShifted[spec.shape[1]-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax6
ax.set_title('1e) FT Spectrum (shifted): imag part')
im=ax.imshow(spec_ftShifted.imag.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,1,int(spec.shape[0]/2),int(spec.shape[0]/2+1),spec.shape[0]-1])
ax.set_xticklabels([xFreqShifted[0],xFreqShifted[1],xFreqShifted[int(spec.shape[0]/2)],xFreqShifted[int(spec.shape[0]/2+1)],xFreqShifted[spec.shape[0]-1]])
ax.set_yticks([0,1,int(spec.shape[1]/2),int(spec.shape[1]/2+1),spec.shape[1]-1])
ax.set_yticklabels([yFreqShifted[0],yFreqShifted[1],yFreqShifted[int(spec.shape[1]/2)],yFreqShifted[int(spec.shape[1]/2+1)],yFreqShifted[spec.shape[1]-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)

plt.show()

输出:

2。实值 FFT

实值输入数组 (rfft) 的 FT 利用了 FT 必须是 Hermitian 的事实。函数 F 是 Hermitian 函数，当:F(f)=F(-f)*
因此只需要存储一半大小的数组。小心:我们会在3.中看到使用rfft时有一些香蕉皮。

代码示例:

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig, ((ax1,ax2,ax3)) = plt.subplots(1,3, figsize=(18.3,12.5))


#Spectrum and FT
spec=np.ndarray(shape=(10,8))
for ii in range(0,spec.shape[0]):
    for jj in range(0,spec.shape[1]):
        spec[ii,jj]=np.cos(2*np.pi*ii/spec.shape[0])*np.cos(4*np.pi*jj/spec.shape[1])
        #spec[ii,jj]=(-1)**(ii+jj)
        #spec[:,jj]=np.cos(2*np.pi*jj/spec.shape[1])
        #spec[ii,:]=np.cos(2*np.pi*ii/spec.shape[0])
        
spec_ft=np.fft.fftshift(np.fft.rfft2(spec),axes=0)

#Frequencies / labels axes
xFreq=np.around(np.fft.fftshift(np.fft.fftfreq(spec.shape[0])),decimals=2)
yFreq=np.around(np.fft.rfftfreq(spec.shape[1]),decimals=2)


#Plotting   
ax=ax1
ax.set_title('2a) Real spectrum')
im=ax.imshow(spec.T,origin='lower',cmap='bwr',interpolation='none')
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax2
ax.set_title('2b) FT: real part')
im=ax.imshow(spec_ft.real.T,origin='lower',cmap='afmhot',interpolation='none')
ax.set_xticks([0,1,int(len(xFreq)/2),int(len(xFreq)/2+1),len(xFreq)-1])
ax.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(len(xFreq)/2)],xFreq[int(len(xFreq)/2+1)],xFreq[len(xFreq)-1]])
ax.set_yticks([0,1,int(len(yFreq)/2),int(len(yFreq)/2+1),len(yFreq)-1])
ax.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(len(yFreq)/2)],yFreq[int(len(yFreq)/2+1)],yFreq[len(yFreq)-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax3.set_title('2c) FT: imag part')
im3=ax3.imshow(spec_ft.imag.T,origin='lower',cmap='bwr',interpolation='none')
ax3.set_xticks([0,1,int(len(xFreq)/2),int(len(xFreq)/2+1),len(xFreq)-1])
ax3.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(len(xFreq)/2)],xFreq[int(len(xFreq)/2+1)],xFreq[len(xFreq)-1]])
ax3.set_yticks([0,1,int(len(yFreq)/2),int(len(yFreq)/2+1),len(yFreq)-1])
ax3.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(len(yFreq)/2)],yFreq[int(len(yFreq)/2+1)],yFreq[len(yFreq)-1]])
divider0 = make_axes_locatable(ax3)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im3,cax=cax)
    
plt.show()

输出:

3。实值逆FFT

rfft 的问题是，它没有歧义。这意味着，从 FT ndarry 的像素数 n_ft 中，无法推断出实际空间中的数组 n_rs 有多少个像素。有2 个选项:

n_rs= (n_ft-1)*2或者n_rs=(n_ft-1)*2+1

图 3a 和 3d 显示了 2 个实值二维数组。我们看到他们的 FT 看起来完全一样(3b，3e)。虚部(3c 和 3f)几乎都为 0。因此，如果我们不知道 FT 谱的维数/形状，我们就无法从 FT 谱得出真实空间谱是什么样子的结论。 因此，使用标准 fft 而不是 rfft 可能更通用。

警告:irfft 不明确!!!

代码示例:

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2,3,figsize=(20,10))

#Spectrum and FT
#generating spectrum: this is just used to generate real spectrum option 1/2 and is not plotted
#it is equivalent to: spec_ft1 and spec_ft2
spec_ft=np.zeros(shape=(10,5))
spec_ft[4,2]=20
spec_ft[6,2]=20
#
spec_r1=np.fft.irfft2(np.fft.ifftshift(spec_ft,axes=0),s=(spec_ft.shape[0],(spec_ft.shape[1]-1)*2))
spec_r2=np.fft.irfft2(np.fft.ifftshift(spec_ft,axes=0),s=(spec_ft.shape[0],(spec_ft.shape[1]-1)*2+1))
#
spec_ft1=np.fft.fftshift(np.fft.rfft2(spec_r1),axes=0)
spec_ft2=np.fft.fftshift(np.fft.rfft2(spec_r2),axes=0)

#Frequencies / labels axes
xFreq1=np.around(np.fft.fftshift(np.fft.fftfreq(spec_r1.shape[0])),decimals=2)
yFreq1=np.around(np.fft.rfftfreq(spec_r1.shape[1]),decimals=2)
xFreq2=np.around(np.fft.fftshift(np.fft.fftfreq(spec_r2.shape[0])),decimals=2)
yFreq2=np.around(np.fft.rfftfreq(spec_r2.shape[1]),decimals=2)

#Plotting  
ax=ax1
ax.set_title('3a) Spectrum 1')
im=ax.imshow(spec_r1.real.T,origin='lower',cmap='bwr',interpolation='none')
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax2
ax.set_title('3b) FT Spectrum 1: real part')
im=ax.imshow(spec_ft1.real.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,int(len(xFreq1)/2),len(xFreq1)-1])
ax.set_xticklabels([xFreq1[0],xFreq1[int(len(xFreq1)/2)],xFreq1[-1]])
ax.set_yticks([0,len(yFreq1)-1])
ax.set_yticklabels([yFreq1[0],yFreq1[-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax3
ax.set_title('3c) FT Spectrum 1: imag part')
im=ax.imshow(spec_ft1.imag.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,int(len(xFreq1)/2),len(xFreq1)-1])
ax.set_xticklabels([xFreq1[0],xFreq1[int(len(xFreq1)/2)],xFreq1[-1]])
ax.set_yticks([0,len(yFreq1)-1])
ax.set_yticklabels([yFreq1[0],yFreq1[-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax4
ax.set_title('3d) Spectrum 2')
im=ax.imshow(spec_r2.real.T,origin='lower',cmap='bwr',interpolation='none')
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax5
ax.set_title('3e) FT Spectrum 2: real part')
im=ax.imshow(spec_ft2.real.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,int(len(xFreq2)/2),len(xFreq2)-1])
ax.set_xticklabels([xFreq2[0],xFreq2[int(len(xFreq2)/2)],xFreq2[-1]])
ax.set_yticks([0,len(yFreq2)-1])
ax.set_yticklabels([yFreq2[0],yFreq2[-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax6
ax.set_title('3f) FT Spectrum 2: imag part')
im=ax.imshow(spec_ft2.imag.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,int(len(xFreq2)/2),len(xFreq2)-1])
ax.set_xticklabels([xFreq2[0],xFreq2[int(len(xFreq2)/2)],xFreq2[-1]])
ax.set_yticks([0,len(yFreq2)-1])
ax.set_yticklabels([yFreq2[0],yFreq2[-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
    
plt.show()

输出:

4。快速傅里叶变换

这是一个如何在不了解真实频谱知识的情况下使用 FT 的示例。我们看到FT阵列的形状与真实空间阵列的形状相同。

代码示例:

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig, ((ax1,ax2,ax3)) = plt.subplots(1,3, figsize=(18.3,12.5))

#Spectrum and FT
spec=np.ndarray(shape=(11,8))
for ii in range(0,spec.shape[0]):
    for jj in range(0,spec.shape[1]):
        #spec[ii,jj]=(-1)**(ii+jj)
        #spec[ii,jj]=np.sin(2*np.pi*ii/spec.shape[0])*np.sin(4*np.pi*jj/spec.shape[1])
        spec[ii,jj]=np.cos(2*np.pi*ii/spec.shape[0])*np.cos(4*np.pi*jj/spec.shape[1])
        #spec[:,jj]=np.cos(2*np.pi*jj/spec.shape[1])
        #spec[ii,:]=np.cos(2*np.pi*ii/spec.shape[0])

spec_ft=np.fft.fftshift(np.fft.fft2(spec))

#Frequencies / labels axes
xFreq=np.around(np.fft.fftshift(np.fft.fftfreq(spec.shape[0])),decimals=2)
yFreq=np.around(np.fft.fftshift(np.fft.fftfreq(spec.shape[1])),decimals=2)

#Plotting
ax=ax1
ax.set_title('4a) Real spectrum')
im=ax.imshow(spec.T,origin='lower',cmap='bwr',interpolation='none')
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax2
ax.set_title('4b) FT: real part')
im=ax.imshow(spec_ft.real.T,origin='lower',cmap='afmhot',interpolation='none')
ax.set_xticks([0,1,int(len(xFreq)/2),int(len(xFreq)/2+1),len(xFreq)-1])
ax.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(len(xFreq)/2)],xFreq[int(len(xFreq)/2+1)],xFreq[len(xFreq)-1]])
ax.set_yticks([0,1,int(len(yFreq)/2),int(len(yFreq)/2+1),len(yFreq)-1])
ax.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(len(yFreq)/2)],yFreq[int(len(yFreq)/2+1)],yFreq[len(yFreq)-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax3
ax.set_title('4c) FT: imag part')
im=ax.imshow(spec_ft.imag.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,1,int(len(xFreq)/2),int(len(xFreq)/2+1),len(xFreq)-1])
ax.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(len(xFreq)/2)],xFreq[int(len(xFreq)/2+1)],xFreq[len(xFreq)-1]])
ax.set_yticks([0,1,int(len(yFreq)/2),int(len(yFreq)/2+1),len(yFreq)-1])
ax.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(len(yFreq)/2)],yFreq[int(len(yFreq)/2+1)],yFreq[len(yFreq)-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
    
plt.show()

输出:

5。逆FFT

当我们执行逆 FT 时，得到的实空间 ndarray 将是一个复数。由于数值精度的原因，这些值将具有一些无穷小的虚部。

代码示例:

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig, ((ax1,ax2,ax3)) = plt.subplots(1,3, figsize=(18.3,12.5))

#Spectrum and FT
spec_ft=np.zeros(shape=(10,8))
spec_ft[4,2]=20
spec_ft[4,6]=20
spec_ft[6,2]=20
spec_ft[6,6]=20

spec_ift=np.fft.ifft2(np.fft.ifftshift(spec_ft))
        

#Frequencies / labels axes
xFreq=np.around(np.fft.fftshift(np.fft.fftfreq(spec_ift.shape[0])),decimals=2)
yFreq=np.around(np.fft.fftshift(np.fft.fftfreq(spec_ift.shape[1])),decimals=2)

#Plotting
ax=ax1
ax.set_title('FT: real part')
im=ax.imshow(spec_ft.real.T,origin='lower',cmap='afmhot',interpolation='none')
ax.set_xticks([0,1,int(len(xFreq)/2),int(len(xFreq)/2+1),len(xFreq)-1])
ax.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(len(xFreq)/2)],xFreq[int(len(xFreq)/2+1)],xFreq[len(xFreq)-1]])
ax.set_yticks([0,1,int(len(yFreq)/2),int(len(yFreq)/2+1),len(yFreq)-1])
ax.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(len(yFreq)/2)],yFreq[int(len(yFreq)/2+1)],yFreq[len(yFreq)-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax2
ax.set_title('FT: imag part')
im=ax.imshow(spec_ft.imag.T,origin='lower',cmap='bwr',interpolation='none')
ax.set_xticks([0,1,int(len(xFreq)/2),int(len(xFreq)/2+1),len(xFreq)-1])
ax.set_xticklabels([xFreq[0],xFreq[1],xFreq[int(len(xFreq)/2)],xFreq[int(len(xFreq)/2+1)],xFreq[len(xFreq)-1]])
ax.set_yticks([0,1,int(len(yFreq)/2),int(len(yFreq)/2+1),len(yFreq)-1])
ax.set_yticklabels([yFreq[0],yFreq[1],yFreq[int(len(yFreq)/2)],yFreq[int(len(yFreq)/2+1)],yFreq[len(yFreq)-1]])
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
#
ax=ax3
ax.set_title('real cosine spectrum')
im=ax.imshow(spec_ift.real.T,origin='lower',cmap='bwr',interpolation='none')
divider0 = make_axes_locatable(ax)
cax=divider0.append_axes("right", size="5%", pad=0.05)
cbar=plt.colorbar(im,cax=cax)
        
plt.show()

输出: