np.fft
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
numpy.rad2deg(array, out)：将弧度转换为度。参数：
array：要计算其度值的角度。deg：输出数组的形状。返回值：
- 返回一个数组，包含弧度对应的角度。
import numpy as np import math arr = [0, math.pi/2, math.pi/4, math.pi/6 ] print ("弧度值 : \n", arr) degval = np.rad2deg(arr) print ("角度值 : \n", degval)弧度值 : [0, 1.5707963267948966, 0.7853981633974483, 0.5235987755982988] 角度值 : [ 0. 90. 45. 30.]
2. np.fft.fft()
numpy.fft.fft(a, n=None, axis=-1, norm=None)：计算一维离散傅里叶变换。
参数：
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\)）。例如，样本间隔为秒，则频率单位为周期/秒。
给定窗口长度n和d：
示例：
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]
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.]]
