1. FFT 知识

傅里叶变换( \(Fourier\ Transform,FT\) ) 是一种线性积分变换,用于信号在时域(或空域)到频域之间的变换。

\(FFT\) 变换( \(Fast\ Fourier\ Transform,快速傅里叶变换\) ) 是针对一组数值进行运算,这组数的长度 \(N\) 必须是 \(2\) 的整数次幂,例如: \(64、128、256\) 等。

数值可以是实数也可以复数,通常我们的时域信号都是实数。

对这 \(N\) 个实数值进行 \(FFT\) 变换,得到一个有 \(N\) 个复数的数组,称此 复数数组为频域信号。 复数数组有如下规律:

  • 下标为 \(0\) \(N/2\) 的两个复数的虚数部分为 \(0\)
  • 下标为 \(i\) \(N-i\) 的两个复数共轭,也就是虚数部分值相同、符号相反。
  • 示例: \(rand\) 随机生成实数数组 \(x\) ,对其 \(fft\) 变换,结果为 \(8\) 个复数。

    import numpy as np
    x = np.random.rand(8)
    print(x)
    xf = np.fft.fft(x)
    print(xf)
    
    [0.66622462 0.49234934 0.40298128 0.48832192 0.76572852 0.18502897
     0.49050222 0.36262891]
    [ 3.85376579+0.j          0.02892604-0.21866575j  0.53846964+0.17357253j
     -0.22793384-0.39370763j  0.79710752+0.j         -0.22793384+0.39370763j
      0.53846964-0.17357253j  0.02892604+0.21866575j]
    

    通过 \(ifft\) 变换(逆 \(fft\) 变换)还原:

    print(np.fft.ifft(xf))
    
    [0.66622462+0.j 0.49234934+0.j 0.40298128+0.j 0.48832192+0.j
     0.76572852+0.j 0.18502897+0.j 0.49050222+0.j 0.36262891+0.j]
    

    \(ifft\) 运算结果实际与 \(x\) 相同,由于浮点数运算误差,出现非常小的虚部。

    \(FFT\) 变换后复数代表什么意思?

  • 下标为 \(0\) 的实数表示时域信号中的直流成分的多少。
  • 下标为 \(i\) 的复数 \(a+b*j\) 表示时域信号中周期为 \(N/i\) 个取样值的正弦波和余弦波的成分的多少,其中 \(a\) 表示 \(cos\) 波形的成分,\(b\) 表示 \(sin\) 波形的成分。
  • 示例:对直流信号进行 \(FFT\) 变换。

    x = np.ones(8)
    print(x)
    print(np.fft.fft(x) / len(x))
    
    [1. 1. 1. 1. 1. 1. 1. 1.]
    [1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]
    

    所谓直流信号,就是其值不随时间变化。

    为了计算各个成分的能量多少,将 \(FFT\) 的结果除以长度。

    ① 对周期为 \(8\) 个取样的正、余弦波进行 \(FFT\) 变换。

    x = np.arange(0, 2*np.pi, 2*np.pi/8)
    y = np.sin(x)
    print(np.fft.fft(y)/len(y))
    z = np.cos(x)
    print(np.fft.fft(z)/len(z))
    
    [ 1.53080850e-17+0.00000000e+00j -4.30636606e-17-5.00000000e-01j
      1.53080850e-17-2.77555756e-17j  1.24474906e-17+0.00000000e+00j
      1.53080850e-17+0.00000000e+00j  1.24474906e-17+0.00000000e+00j
      1.53080850e-17+2.77555756e-17j -4.30636606e-17+5.00000000e-01j]
    [-4.30636606e-17+0.00000000e+00j  5.00000000e-01-5.83717456e-17j
      1.53080850e-17+0.00000000e+00j  0.00000000e+00+2.86059436e-18j
      1.24474906e-17+0.00000000e+00j  0.00000000e+00-2.86059436e-18j
      1.53080850e-17+0.00000000e+00j  5.00000000e-01+5.83717456e-17j]
    

    正弦波的 \(FFT\) 结果:下标为 \(1\) 的复数的虚部为 \(-0.5\),产生的正弦波的放大系数(振幅)为 \(1\),它们之间的关系是 \(-0.5*(-2)=1\)

    余弦波的 \(FFT\) 结果:只有下标为 \(1\) 的复数的实数部分有有效数值 \(0.5\),和余弦波的放大系数(振幅) \(1\) 之间的关系是 \(0.5*2=1\)

    ② 对周期为 \(4\) 个取样的正、余弦波进行 \(FFT\) 变换。

    print(np.fft.fft(2*np.sin(2*x))/len(x))
    print(np.fft.fft(0.8*np.cos(2*x))/len(x))
    
    [ 6.1232340e-17+0.000000e+00j  6.1232340e-17+6.123234e-17j
     -1.8369702e-16-1.000000e+00j  6.1232340e-17-6.123234e-17j
      6.1232340e-17+0.000000e+00j  6.1232340e-17+6.123234e-17j
     -1.8369702e-16+1.000000e+00j  6.1232340e-17-6.123234e-17j]
    [-2.44929360e-17+0.00000000e+00j -3.46382422e-17-3.08148791e-33j
      4.00000000e-01-9.79717439e-17j  3.46382422e-17-3.08148791e-33j
      2.44929360e-17+0.00000000e+00j  3.46382422e-17+3.08148791e-33j
      4.00000000e-01+9.79717439e-17j -3.46382422e-17+3.08148791e-33j]
    

    \(FFT\) 的有效成分在下标为 \(2\) 的复数中。

    正弦波的 \(FFT\) 结果:正弦波的放大系数(振幅)为 \(2\),因此频域虚数部分的值为 \(-1\)

    余弦波的 \(FFT\) 结果:余弦波的放大系数(振幅)为 \(0.8\),因此其对应的值为 \(0.4\)

    同频率的正弦波和余弦波通过不同的系数叠加,可以产生同频率的各种相位的余弦波,因此我们可以这样来理解频域中的复数:

  • 复数的模(绝对值)代表此频率的余弦波的振幅。
  • 复数的辐角代表了此频率的余弦波的相位。
  • ③ 最后一个例子。

    x = np.arange(0, 2*np.pi, 2*np.pi/128)
    y = 0.3*np.cos(x) + 0.5*np.cos(2*x+np.pi/4) + 0.8*np.cos(3*x-np.pi/3)
    yf = np.fft.fft(y)/len(y)
    print(yf[:4])
    print(np.angle(yf[1]))
    print(np.abs(yf[1]), np.rad2deg(np.angle(yf[1])))
    print(np.abs(yf[2]), np.rad2deg(np.angle(yf[2])))
    print(np.abs(yf[3]), np.rad2deg(np.angle(yf[3])))
    
    [1.22514845e-17+0.00000000e+00j
     1.50000000e-01-5.20417043e-18j
     1.76776695e-01+1.76776695e-01j
     2.00000000e-01-3.46410162e-01j]
    -3.469446951953613e-17
    0.15000000000000005 -1.987846675914697e-15
    0.2500000000000001 45.000000000000014
    0.39999999999999997 -60.00000000000002
    

    上例产生了三个不同频率的余弦波,并且给他们不同的振幅和相位:

  • 周期为 \(128/1\) 点的余弦波的相位为 \(0\), 振幅为 \(0.3\)
  • 周期为 \(64/2\) 点的余弦波的相位为 \(45\) 度, 振幅为 \(0.5\)
  • 周期为 \(128/3\) 点的余弦波的相位为 \(-60\) 度,振幅为 \(0.8\)
  • numpy.angle(z, deg=0):计算复数的辐角主值。

  • z:复数或复数组成的列表。
  • deg:若为 \(False\)(默认),返回值是弧度值。若为 \(True\),返回值是角度值。
  • \(n\) 维数组或标量。复平面上从正实半轴出发沿逆时针方向到复数所在点走过的角度。
  • import numpy as np
    print('结果用弧度制表示:{}'.format(np.angle(1+1j)))
    print('结果用角度制表示:{}'.format(np.angle(1+1j, deg=True)))
    
    结果用弧度制表示:0.7853981633974483
    结果用角度制表示:45.0
    
  • a:输入数组。
  • n:输出的变换轴的长度。如果 \(n\) 小于输入的长度,则裁剪输入。如果 \(n\) 小于输入的长度,则用 \(0\) 填充。未给出 \(n\),则使用沿轴指定的轴的输入长度。
  • axis:计算 \(FFT\) 的轴。如果未给出,则使用最后一个轴。
  • norm:标准化模式。取值为{“backward”, “ortho”, “forward”}。默认设置是backward。指示缩放前向/后向转换对应的哪个方向以及使用什么标准化因子。
  • 截断或补零的输入,沿轴指示的轴进行转换,如果未指定轴,则为最后一个轴。
  • import numpy as np
    t = np.arange(8) / 8
    a = np.fft.fft(np.exp(2j * np.pi * t))
    print(a)
    
    [-3.44509285e-16+1.22464680e-16j  8.00000000e+00-8.11483250e-16j
      3.44509285e-16+1.22464680e-16j  0.00000000e+00+1.22464680e-16j
      9.95799250e-17+1.22464680e-16j  0.00000000e+00+7.66951701e-17j
     -9.95799250e-17+1.22464680e-16j  0.00000000e+00+1.22464680e-16j]
    

    3. np.fft.fft2()

    numpy.fft.fft2(a, s=None, axes=(- 2, - 1), norm=None)计算二维离散傅里叶变换。

    此函数通过快速傅里叶变换 (\(FFT\)) 计算 \(M\) 维数组中任意轴上的 \(2\) 维离散傅里叶变换。默认情况下,变换是在输入数组的最后两个轴上计算的,即二维 \(FFT\)

  • a:输入数组。
  • s:整数序列,可选。
  • ​ 输出的形状(每个变换轴的长度)(\(s[0]\) 指轴 \(0\)\(s[1]\) 指轴 \(1\))。

    ​ 如果 \(n\) 小于输入的长度,则裁剪输入。如果 \(n\) 小于输入的长度,则用 \(0\) 填充。未给出 \(n\),则使用沿轴指定的轴的输入长度。

  • axes:整数序列,可选。
  • ​ 计算 \(FFT\) 的轴。如果未给出,则使用最后两个轴。

  • norm:标准化模式。取值为{“backward”, “ortho”, “forward”}。默认设置是backward。指示缩放前向/后向转换对应的哪个方向以及使用什么标准化因子。
  • 截断或补零的输入,沿轴指示的轴进行转换,如果未指定轴,则为最后两个轴转换。
  • import numpy as np
    a = np.mgrid[:5, :5][0]
    print(a)
    print(np.fft.fft2(a))
    
    [[0 0 0 0 0]
     [1 1 1 1 1]
     [2 2 2 2 2]
     [3 3 3 3 3]
     [4 4 4 4 4]]
    [[ 50.  +0.j           0.  +0.j           0.  +0.j
        0.  +0.j           0.  +0.j        ]
     [-12.5+17.20477401j   0.  +0.j           0.  +0.j
        0.  +0.j           0.  +0.j        ]
     [-12.5 +4.0614962j    0.  +0.j           0.  +0.j
        0.  +0.j           0.  +0.j        ]
     [-12.5 -4.0614962j    0.  +0.j           0.  +0.j
        0.  +0.j           0.  +0.j        ]
     [-12.5-17.20477401j   0.  +0.j           0.  +0.j
        0.  +0.j           0.  +0.j        ]]
    

    4. np.fft.fftfreq

    numpy.fft.fftfreq(n, d=1.0)返回离散傅里叶变换采样频率

  • n:窗口长度。
  • d:采样间隔。默认为 \(1\)
  • 长度为 \(n\) 的数组,包含采样频率。
  • 返回浮点数组f包含频率单元中心,该频率单元中心以每单位采样间隔的周期为周期(开始为 \(0\))。例如,样本间隔为秒,则频率单位为周期/秒。

    给定窗口长度nd

    \[\begin{aligned} & f = [0, 1, ...,\frac{n}{2}-1,-\frac{n}{2},...,-1] / (d*n) \ \ \ \ 如果n是偶数 \\ & f = [0,1,...,\frac{n-1}{2},-\frac{n-1}{2},...,-1] / (d*n) \ \ \ \ 如果n是奇数 \end{aligned} \]

    示例:

    import numpy as np
    from numpy.fft import *
    signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
    fourier = np.fft.fft(signal)
    n, timestep = signal.size, 0.1            # 窗口长度n=8,采样间隔timestep=0.1
    freq = np.fft.fftfreq(n, d=timestep)
    print(freq)
    
    [ 0.    1.25  2.5   3.75 -5.   -3.75 -2.5  -1.25]
    

    \[\begin{aligned} & 0 = 0 \\ & 1.25 = 1 / (0.1*8) \\ & 2.5 = 2 / (0.1*8) \\ & 3.75 = 3 / (0.1*8) \\ & -5 = -4 / (0.1*8) \\ & -3.75 = -3 / (0.1*8) \\ & -2.5 = -2 / (0.1*8) \\ & -1.25 = -1 / (0.1*8) \end{aligned}

    5. np.fft.fftshift

    numpy.fft.fftshift(x, axes=None):将零频率分量移动到频谱中心。即中心化。

  • x:输入数组。
  • d:要移动的轴。默认为无,表示移动所有轴。
  • 移动的数组。
  • freqs = np.fft.fftfreq(10, 0.1)
    print(freqs)
    print(np.fft.fftshift(freqs))
    
    [ 0.  1.  2.  3.  4. -5. -4. -3. -2. -1.]
    [-5. -4. -3. -2. -1.  0.  1.  2.  3.  4.]
    

    仅沿第二轴移动零频分量。

    freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
    print(freqs)
    print(np.fft.fftshift(freqs))				# 沿所有轴移动零频分量
    print(np.fft.fftshift(freqs, axes=(1,)))	        # 沿第二轴移动零频分量
    
    [[ 0.  1.  2.]
     [ 3.  4. -4.]
     [-3. -2. -1.]]
    [[-1. -3. -2.]
     [ 2.  0.  1.]
     [-4.  3.  4.]]
    [[ 2.  0.  1.]
     [-4.  3.  4.]
     [-1. -3. -2.]]
    

    6. np.fft.ifftshift

    numpy.fft.ifftshift(x, axes=None)去中心化。

    freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
    print(freqs)
    print(np.fft.ifftshift(np.fft.fftshift(freqs)))
    
    [[ 0.  1.  2.]
     [ 3.  4. -4.]
     [-3. -2. -1.]]
    [[ 0.  1.  2.]
     [ 3.  4. -4.]
     [-3. -2. -1.]]