数字图像处理-13-图像频域变换数学基础之快速傅里叶变换

前面已经介绍过快速傅里叶变换的实现原理,这里用python代码实现下

1.一维快速傅里叶变换 (1D-FFT)

一维快速傅里叶变换 (1D-FFT)，Cooley-Tukey 时域抽取(DIT)算法。

数学定义（DFT）： X[k] = Σ_{n=0}^{N-1} x[n] · e^{-j2πkn/N}, k = 0,1,...,N-1

• 算法原理：

  将长度 N 的 DFT 分解为两个长度 N/2 的 DFT（分治策略）： X_偶[k] = Σ x[2n] · e^{-j2πkn/(N/2)} X_奇[k] = Σ x[2n+1] · e^{-j2πkn/(N/2)} 合并： X[k] = X_偶[k] + W_N^k · X_奇[k] X[k + N/2] = X_偶[k] - W_N^k · X_奇[k] 其中旋转因子 W_N^k = e^{-j2πk/N} 递归展开后得到 log₂N 级蝴蝶运算，总复杂度 O(N log N)。

• 实现步骤：
  1. 位倒序排列输入（bit-reversal permutation）
  2. 执行 log₂N 级蝴蝶运算（每级使用 numpy reshape 向量化所有组）

• python代码

# ===========================================================================# 辅助函数# ===========================================================================def_next_pow2(n:int)->int:"""返回大于等于 n 的最小 2 的幂次。用于对图像尺寸做补零填充。"""p=1whilep<n:p<<=1returnpdef_bit_rev_indices(N:int,power:int)->np.ndarray:""" 计算长度 N = 2^power 的位倒序排列索引数组。 FFT 蝴蝶运算前需将输入按位倒序重排（bit-reversal permutation）。 例如 power=3： 原序 0(000)→0, 1(001)→4, 2(010)→2, 3(011)→6, 4(100)→1, 5(101)→5, 6(110)→3, 7(111)→7 """idx=np.empty(N,dtype=np.int32)foriinrange(N):x,r=i,0for_inrange(power):r=(r<<1)|(x&1)x>>=1idx[i]=rreturnidx# ===========================================================================# 快速傅里叶变换 (FFT) — Cooley-Tukey DIT 算法# ===========================================================================deffft1d(x:np.ndarray)->np.ndarray:""" 参数： x: 长度必须是 2 的幂次的一维 numpy 数组（实数或复数） 返回： 频域复数数组，与输入等长 """N=len(x)power=int(round(math.log2(N)))assert(1<<power)==N,f"输入长度{N}必须是 2 的幂次"# 步骤 1：位倒序排列rev=_bit_rev_indices(N,power)X=np.asarray(x,dtype=np.complex128)[rev]# 步骤 2：逐级蝴蝶运算# half: 当前半块长度；m=2*half: 当前整块长度,这里做的就是将元素递归2分half=1whilehalf<N:m=half<<1# 整块长度k=np.arange(half)# 旋转因子 W[k] = e^{-j·2π·k/m} = cos(-2πk/m) + j·sin(-2πk/m)angle=-2.0*math.pi*k/m W=np.cos(angle)+1j*np.sin(angle)# 将 X 重塑为 (N/m, m)，对每一行同时做蝴蝶运算# 这里会循环将数字序列按2,4,8 2^n的方式重新费城多列,比如下面如果m=2的话,则说明把# X中元素分成2列,-1这里代表不知道元素的个数,系统会自动推断长度.X=X.reshape(-1,m)u=X[:,:half].copy()# 上半部分t=X[:,half:]*W# 下半部分乘旋转因子X[:,:half]=u+t# 蝴蝶和X[:,half:]=u-t# 蝴蝶差X=X.reshape(-1)#将数组 X 展平（Flatten）为一个一维数组。half=mreturnX

• 效果图
  

2. 一维快速傅里叶逆变换 (1D-IFFT)

• 原理
  利用共轭对称性，逆变换可由 FFT 导出：

  IDFT(X)[n] = (1/N) · conj( FFT( conj(X) ) )[n]

  只需调用正向 FFT，无需重新实现。

defifft1d(X:np.ndarray)->np.ndarray:""" 只需调用正向 FFT，无需重新实现。 """N=len(X)returnnp.conj(fft1d(np.conj(X)))/N

3. 二维快速傅里叶变换 (2D-FFT)

• 原理（可分离性）：

  二维 DFT 可分解为行方向和列方向两次一维 DFT： F(u,v) = Σ_y [ e^{-j2πvy/N} · Σ_x f(x,y)·e^{-j2πux/M} ] 先对每行做 1D-FFT，再对每列做 1D-FFT，结果相同。

  非 2 的幂次尺寸自动向上补零到 2 的幂次。 即这里我们是按N^2的元素个数序列进行变换

deffft2d(image:np.ndarray)->tuple:""" 参数： image: 二维 numpy 数组（灰度图像，uint8 或 float） 返回： (F, pad_rows, pad_cols) F : 复数频谱矩阵，尺寸 (pad_rows, pad_cols) pad_rows : 行方向填充后尺寸 pad_cols : 列方向填充后尺寸 """rows,cols=image.shape pr=_next_pow2(rows)pc=_next_pow2(cols)# 零填充到 2 的幂次F=np.zeros((pr,pc),dtype=np.complex128)F[:rows,:cols]=image.astype(np.float64)# 行方向 FFTforiinrange(pr):F[i,:]=fft1d(F[i,:])# 列方向 FFT（转置后按行处理，再转置回来）F=F.T.copy()foriinrange(pc):F[i,:]=fft1d(F[i,:])F=F.TreturnF,pr,pc

4. 二维快速傅里叶逆变换 (2D-IFFT)

defifft2d(F:np.ndarray)->np.ndarray:""" 二维快速傅里叶逆变换 (2D-IFFT)。 先对每行做 1D-IFFT，再对每列做 1D-IFFT。 """rows,cols=F.shape G=F.copy()# 行方向 IFFTforiinrange(rows):G[i,:]=ifft1d(G[i,:])# 列方向 IFFTG=G.T.copy()foriinrange(cols):G[i,:]=ifft1d(G[i,:])G=G.TreturnG

5. 一维离散余弦变换 (1D-DCT-II)

一维离散余弦变换（1D DCT）是数字信号处理中一种至关重要的技术。它能够将一个长度为N 的离散信号序列从空间域（或时域）转换到频域。
在数学原理上，DCT 的核心思想是将有限长度的信号展开为不同频率的余弦基函数的加权和。由于自然界的实数信号通常具有偶对称特性，因此其傅里叶变换只包含余弦项，这使得 DCT 相比于涉及复数的离散傅里叶变换（DFT），仅需使用实数进行计算，极大地简化了运算复杂度。

这里由于sin函数是奇函数,所以其在一个周期内和为0,根据欧拉公式可知

e^iΘ = cos(Θ) + j sin(Θ)

根据对偶性,所以上面的旋转因子, W[k] = e^{-j·2π·k/m} = cos(-2πk/m) + j·sin(-2πk/m)
可以去掉奇函数部分.

• 数学定义（DCT-II）：

  C[0] = (1/√N) · Σ_{n=0}^{N-1} x[n] C[k] = √(2/N) · Σ_{n=0}^{N-1} x[n] · cos( π·k·(2n+1) / (2N) ), k≥1

即:

• FFT 加速原理:

  构造长度 2N 的序列 y： y[n] = x[n], n = 0,...,N-1 y[n] = 0, n = N,...,2N-1 计算 Y = FFT(y)，则 DCT 系数可由 Y 导出： C[0] = Re(Y[0]) / √N C[k] = (Re(Y[k])·cos(πk/2N) + Im(Y[k])·sin(πk/2N)) · √(2/N), k≥1 这利用了 DCT 与 DFT 之间的代数关系，避免重新实现。

• DCT 相比 DFT 的优势：

  • 纯实数变换（输入输出均为实数，无虚部）
  • 更强的能量集中性（低频系数携带大部分信息）
  • 不存在 Gibbs 振铃现象（余弦基函数与边界自然兼容）
  • JPEG 图像压缩标准的核心变换

• python代码

# ===========================================================================# 离散余弦变换 (DCT-II) — 基于 FFT 的快速算法# ===========================================================================defdct1d(x:np.ndarray)->np.ndarray:""" 一维离散余弦变换 (1D-DCT-II)，基于 FFT 的快速实现。 参数： x: 长度必须是 2 的幂次的一维实数数组 返回： 实数 DCT 系数数组，与输入等长 """N=len(x)power=int(round(math.log2(N)))assert(1<<power)==N,f"输入长度{N}必须是 2 的幂次"# 构造长度 2N 的序列（后半填零）y=np.zeros(2*N,dtype=np.complex128)y[:N]=x.astype(np.float64)# 2N 点 FFTY=fft1d(y)# 从 FFT 结果提取 DCT 系数k=np.arange(N,dtype=np.float64)angle=math.pi*k/(2.0*N)cos_a=np.cos(angle)sin_a=np.sin(angle)C=np.empty(N)C[0]=Y[0].real/math.sqrt(N)# 只取 FFT 输出的前 N 项（2N 点 FFT 只需前 N 项来重建 N 个 DCT 系数）C[1:]=(Y[1:N].real*cos_a[1:]+Y[1:N].imag*sin_a[1:])*math.sqrt(2.0/N)returnC

6. 二维离散余弦变换 (2D-DCT-II)

• 原理

  利用可分离性： 先对每行做 1D-DCT，再对每列做 1D-DCT， 等价于二维 DCT： C[u,v] = a[u]·a[v] · ΣΣ f(x,y)·cos(πu(2x+1)/2M)·cos(πv(2y+1)/2N) 非 2 的幂次尺寸自动向上补零。

• python代码

defdct2d(image:np.ndarray)->np.ndarray:""" 返回： 实数 DCT 系数矩阵 """rows,cols=image.shape pr=_next_pow2(rows)pc=_next_pow2(cols)F=np.zeros((pr,pc),dtype=np.float64)F[:rows,:cols]=image.astype(np.float64)# 行方向 DCTforiinrange(pr):F[i,:]=dct1d(F[i,:])# 列方向 DCTF=F.T.copy()foriinrange(pc):F[i,:]=dct1d(F[i,:])F=F.TreturnFdefdct_spectrum_to_image(C:np.ndarray)->np.ndarray:"""将 DCT 系数矩阵转换为可显示的 uint8 灰度图像（对数幅度）。"""mag=np.abs(C)mag=np.log1p(mag)lo,hi=mag.min(),mag.max()ifhi>lo:mag=(mag-lo)/(hi-lo)*255.0else:mag=np.zeros_like(mag)returnmag.astype(np.uint8)

• 效果图