# fft
傅里叶变换 (一维、二维、多维)
函数库: TyMath
# 语法
Y = fft(X)
Y = fft(X,dim)
# 说明
Y = fft(X) 用快速傅里叶变换 (FFT) 算法计算 X 的离散傅里叶变换 (DFT)。
- 如果 X 是向量，则 fft(X) 返回该向量的傅里叶变换；示例
- 如果 X 是矩阵，则 fft(X) 返回该矩阵的二维傅里叶变换；示例
- 如果 X 是一个多维数组，返回该多维数组的 N 维傅里叶变换。示例
Y = fft(X,dim) 返回对数组 X 沿维度 dim 进行计算傅里叶变换。示例
# 示例
含噪信号
使用傅里叶变换求噪声中隐藏的信号的频率分量。
指定信号的参数，采样频率为 1 kHz，信号持续时间为 1.5 秒。
using TyMath
using TyPlot
Fs = 1000
T = 1/Fs
L = 1500
t = (0:L-1)*T
构造一个信号，其中包含幅值为 0.7 的 50 Hz 正弦量和幅值为 1 的 120 Hz 正弦量。
S = 0.7*sin.(2*pi*50*t) + sin.(2*pi*120*t)
用均值为 0、方差为 4 的白噪声扰乱该信号。
X = S + 2*randn(size(t))
在时域中绘制含噪信号。通过查看信号 X(t) 很难确定频率分量。
plot(1000*t[1:50],X[1:50])
title("均值为0的随机噪声信号")
xlabel("t(milliseconds)")
ylabel("X(t)")
计算信号的傅里叶变换。
Y = fft(X)
计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs.(Y/L)
P1 = P2[1:Int(L/2)+1]
P1[2:end-1] = 2*P1[2:end-1]
定义频域 f 并绘制单侧幅值频谱 P1。与预期相符，由于增加了噪声，幅值并不精确等于 0.7 和 1。一般情况下，较长的信号会产生更好的频率近似值。
f = Fs*(0:(L/2))/L
plot(f,P1)
title("X(t) 的单侧振幅图")
xlabel("f (Hz)")
ylabel("|P1(f)|")
现在，采用原始的、未破坏信号的傅里叶变换并检索精确幅值 0.7 和 1.0。
Y = fft(S)
P2 = abs.(Y/L)
P1 = P2[1:Int(L/2)+1]
P1[2:end-1] = 2*P1[2:end-1]
plot(f,P1)
title("S(t)的单侧振幅谱")
xlabel("f (Hz)")
ylabel("|P1(f)|")
高斯脉冲
将高斯脉冲从时域转换为频域。
定义信号参数和高斯脉冲 X。
using TyMath
using TyPlot
Fs = 100
t = -0.5:1/Fs:0.5
L = length(t)
X = 1/(4*sqrt.(2*pi*0.01))*(exp.(-t.^2/(2*0.01)))
在时域中绘制脉冲。
plot(t,X)
title("时域高斯脉冲")
xlabel("Time (t)")
ylabel("X(t)")
要使用 fft 将信号转换为频域，首先从原始信号长度确定下一个 2 次幂的新输入长度。这将用尾随零填充信号 X 以改善 fft 的性能。
n = 2^nextpow2(L)
将高斯脉冲转换为频域。
X = append!(X,zeros(n-L))
Y = fft(X)
定义频域并绘制唯一频率。
f = Fs*(0:Int(n/2))/n
P = abs.(Y/n) .^ 2
plot(f,P[1:Int(n/2)+1])
title("频域上的高斯脉冲")
xlabel("Frequency (f)")
ylabel("|P(f)|^2")
余弦波
比较时域和频域中的余弦波。
指定信号的参数，采样频率为 1kHz，信号持续时间为 1 秒。
using TyMath
using TyPlot
Fs = 1000
T = 1/Fs
L = 1000
t = (0:L-1)*T
创建一个矩阵，其中每一列代表一个频率经过缩放的余弦波。结果 X 为 1000×3 矩阵。第一列的波频为 50，第二列的波频为 150，第三列的波频为 300。
x1 = cos.(2*pi*50*t)
x2 = cos.(2*pi*150*t)
x3 = cos.(2*pi*300*t)
X = hcat(x1,x2,x3)
在单个图窗中按顺序绘制 X 的每列的前 100 个项，并比较其频率。
for i = 1:3
ax = subplot(3,1,i)
plot(ax,t[1:100],X[1:100,i])
title("Row $i in the Time Domain")
end
出于算法性能的考虑，fft 允许您用尾随零填充输入。在这种情况下，用零填充 X 的每一列，以使每列的长度为比当前长度大的下一个最小的 2 的次幂值。使用 nextpow2 函数定义新长度。
n = 2^nextpow2(L)
指定 dim 参数沿 X 的列（即对每个信号）使用 fft。
dim = 1
计算信号的傅里叶变换。
X = vcat(X,zeros(n-L,3))
Y = fft(X,dim)
计算每个信号的双侧频谱和单侧频谱。
P2 = abs.(Y/L)
P1 = P2[1:Int(n/2)+1,:]
P1[2:end-1,:] = 2*P1[2:end-1,:]
在频域内，为单个图窗中的每一行绘制单侧幅值频谱。
for i = 1:3
ax = subplot(3,1,i)
plot(ax,0:(Fs/n):(Fs/2-Fs/n),P1[1:Int(n/2),i])
title("Row $i in the Time Domain")
end
二维变换
二维傅里叶变换对处理二维信号和其他二维数据（如图像）很有用。
创建并绘制具有重复块的二维数据。
using TyMath
using TyPlot
P = peaks(20)[3]
X = repeat(P,5,10)
imagesc(X)
计算数据的二维傅里叶变换。将零频分量移动到输出的中心，并绘制生成的 100×200 矩阵，它与 X 的大小相同。
Y = fft(X)
imagesc(abs.(fftshift(Y)))
用零填充 X 以计算 128×256 变换。
n1 = 2^nextpow2(100)
n2 = 2^nextpow2(200)
X = hcat(X,zeros(100,n2-200))
X = vcat(X,zeros(n1-100,n2))
Y = fft(X)
imagesc(abs.(fftshift(Y)))
三维变换
您可以使用 fft 函数计算多维数组中每个维度的一维快速傅里叶变换。
创建一个三维信号 X。X 的大小为 20×20×20。
using TyMath
x = [1:20;]
y = [1:20;]'
z = reshape(1:20,1,1,20)
X = cos.(2*pi*0.01*x) .+ sin.(2*pi*0.02*y) .+ cos.(2*pi*0.03*z)
计算该信号的三维傅里叶变换，这也是一个 20×20×20 数组。
Y = fft(X)
# 输入参数
X - 输入数组向量 | 矩阵 | 多维数组
数据类型： Int | Float |Bool
复数支持： 是
dim - 沿其运算的维度正整数标量
沿其运算的维度，指定为正整数标量。如果未指定值，则默认是对数组的 ndims(X) 维傅里叶变换。
fft(X,1) 沿 X 的第一个维度进行傅里叶变换。
数据类型： Int |Bool
# 输出参数
Y - 频域表示向量 | 矩阵 | 多维数组
频域表示，以向量、矩阵或多维数组形式返回。
数据类型： Float | Int
复数支持： 是
# 详细信息
向量的离散傅里叶变换
Y = fft(X) 和 X = ifft(Y) 分别实现傅里叶变换和傅里叶逆变换。对于长度为 n 的 X 和 Y，这些变换定义如下：
其中，
二维傅里叶变换
以下公式定义 m×n 矩阵 X 的离散傅里叶变换 Y：
i 是虚数单位。p 和 j 是值范围从 0 到 m–1 的索引，q 和 k 是值范围从 0 到 n–1 的索引。此公式将 X 和 Y 的索引平移 1 位。
N 维傅里叶变换
N 维数组 X 的离散傅里叶变换 Y 定义为：
每个维度的长度为
# 提示
- fft 的执行时间取决于变换的长度。仅具有小质因数的变换长度的 fft 执行时间明显快于本身是质数或具有较大质因数的变换长度的 fft 执行时间；
- 对于大多数 n 值，实数输入的 DFT 需要的计算时间大致是复数输入的 DFT 计算时间的一半。但是，当 n 有较大的质因数时，速度很少有差别或没有差别。
# 算法
FFT 函数（fft、ifft）基于一个称为 FFTW [1]的库。
# 参考文献
[1] Frigo, M., and S. G. Johnson. “FFTW: An Adaptive Software Architecture for the FFT.” Proceedings of the International Conference on Acoustics, Speech, and Signal Processing. Vol. 3, 1998, pp. 1381-1384.