python中离散傅里叶变换(FT)的教程，技巧和香蕉皮 [英] Tutorial, tricks and banana skins for discrete Fourier transformation (FT) in python

问题描述

在python中进行傅立叶转换时，我对有用的技巧和香蕉皮感兴趣.

I am interested in usefull tricks and and banana skins when doing Fourier transformation in python.

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

Please give introductive code examples to get into the topic, as well as more advise in advanced topics like: frequency filters, continuous FT, high dimensional FT.

推荐答案

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

1.频率排序

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

Getting started with 2 dimensional discrete Fourier transformation in python (numpy):

1. Frequency ordering

One should be carefull with the ordering of frequencies by the numpy.fft library. The fft routine orders first zero, then positive frequencies and then negative ones (Figure 1b, 1c). The function fftshift can be used to establish ordering as expected i.e. from negative to positive values one needs the fftshift routine (Figure 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()

输出:

实值输入数组(rfft)的FT利用了这样一个事实，即FT必须是埃尔米特式的. 在以下情况下，函数F是厄米函数:F(f)= F(-f)*
因此，仅需要存储一半大小的数组. 小心::我们将在3.中看到使用rfft时有一些香蕉皮.

The FT for real valued input arrays (rfft) makes use of the fact, that the FT has to be hermitian. A function F is hermitian, when: F(f)=F(-f)*
Hence only an array of half the size needs to be stored. Carefull: We will see in 3. that there are some banana skins, when using 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()

输出:

rfft的问题在于，它不是模棱两可的.这意味着，从FT数据的像素数n_ft，我们不能得出实数空间n_rs中的阵列有多少像素.有 2个选项:

The problem with rfft is, that it not ambiguous. This means, that from the number of pixels n_ft of a FT ndarry, one can not conclude how many pixels the array in real space n_rs has. There are 2 options:

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

n_rs= (n_ft-1)*2 or n_rs= (n_ft-1)*2+1

图3a和3d显示了2个实值2d数组.我们看到它们的FT看起来完全一样(3b，3e).虚部(3c和3f)都几乎为0. 因此，从FT光谱我们无法得出结论，如果我们不知道其空间尺寸/形状，那么真实的空间光谱将是什么样. 因此，使用标准fft代替rfft可能更普遍.

Figures 3a and 3d show 2 real valued 2d-arrays. We see that their FT looks exactly the same (3b, 3e). The imaginary parts (3c and 3f) are both almost 0. So from the FT spectrum we can not conclude, what the real space spectrum looks like, if we do not know its dimension / shape. For this reason it might be more universal to use the standard fft instead of rfft.

代码示例:

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

输出:

这是一个如何使用FT而不包含真实频谱知识的示例.我们看到FT阵列的形状与真实空间阵列的形状相同.

Here is an example how to use FT without including the knowledge of the spectrum beeing real. We see taht the shape of the FT array is the same like the shape of the real space array.

代码示例:

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

输出:

当我们执行逆FT时，实际空间ndarray将是一个复数. 由于数值的准确性，这些值将具有一些无穷小的虚部.

When we perform the inverse FT, the rsulting real space ndarray will be a complex number. Because of the numeric accuracy the values will have some infinitesimal small imaginary part.

代码示例:

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

输出:

这篇关于python中离散傅里叶变换(FT)的教程，技巧和香蕉皮的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！