什么是快速傅里叶变换？

最新推荐文章于 2025-04-11 11:48:50 发布

CM莫问

最新推荐文章于 2025-04-11 11:48:50 发布

阅读量4.5k

点赞数 24

分类专栏：机器学习模型人工智能算法常见概念文章标签：算法人工智能机器学习 python 傅里叶傅立叶图像处理

本文链接：https://blog.youkuaiyun.com/ChaneMo/article/details/144176789

版权

人工智能算法常见概念同时被 2 个专栏收录

53 篇文章

订阅专栏

机器学习模型

19 篇文章

订阅专栏

一、FFT概念

快速傅里叶变换（Fast Fourier Transform，FFT）是一种高效计算离散傅里叶变换（Discrete Fourier Transform，简称DFT）及其逆变换的算法。DFT和FFT都是音频处理、图像分析、振动分析、无线通信等许多领域中数字信号处理的关键技术。

1、离散傅里叶变换

首先，我们需要知道什么是DFT。DFT是一种将复数序列（通常为实数信号的离散样本）转换为另一个复数序列的数学运算。给定一个长度为 N 的复数序列 x[n]，其离散傅里叶变换 X[k] 为：

$eq?X%5Bk%5D%20%3D%20%5Csum_%7Bn%3D0%7D%5E%7BN-1%7Dx%5Bn%5D%20%5Ccdot%20e%5E%7B-j%20%5Cfrac%7B2%20%5Cpi%20%7D%7BN%7Dnk%7D$

其中，X[k]是一个复数序列，表示原始信号在频域的表示，每个X[k]对应于信号在第k（取值范围0到N-1）个频率分量上的复数幅度和相位；N为序列的长度，每个x[n]表示信号在时间域的第n个样本；j为虚数单位，满足j的平方等于-1； $eq?e%5E%7B-j%20%5Cfrac%7B2%20%5Cpi%20%7D%7BN%7Dnk%7D$ 为复数指数权重，表示第n个样本在第k个频率分量上的权重，其模为1，相位为 $eq?-%5Cfrac%7B2%20%5Cpi%20%7D%7BN%7Dnk$ 。

DFT公式计算了每个频率分量 k 上的复数和，这个和是所有时间域样本 x[n] 乘以相应的复数指数权重的结果。通过这种方式，DFT将时间域信号转换为频域信号，使得我们可以分析信号的频率成分。

2、快速傅里叶变换

FFT是DFT的快速算法，它利用了DFT的对称性、周期性和冗余性，将DFT的计算复杂度从 O(N²) 降低到 O(Nlog⁡N)。FFT的实现有多种算法，其中最常用的是Cooley-Tukey算法。Cooley-Tukey算法基于分而治之的思想，将DFT分解为更小的DFT问题：

分解：将长度为 N 的序列 x[n]分解为偶数索引部分 x[n]（n 为偶数）和奇数索引部分 x[n]（n 为奇数）。
递归：对偶数索引部分和奇数索引部分分别进行DFT，得到 $eq?X_%7Be%7D%5Bk%5D$ 和 $eq?X_%7Bo%7D%5Bk%5D$ 。
合并：利用DFT的性质，将 $eq?X_%7Be%7D%5Bk%5D$ 和 $eq?X_%7Bo%7D%5Bk%5D$ 合并为原始序列的X[k]。

3、FFT示例

这里，我们简单构建一个带噪音的正弦波图像。作为一种周期性波形，正弦波的频率决定了波形的周期性重复速率。

import numpy as np
import matplotlib.pyplot as plt

# 创建一个简单的正弦波信号
t = np.linspace(0, 1, 400, endpoint=False)  # 时间序列
f = 5  # 信号频率（Hz）
signal = np.sin(2 * np.pi * f * t)  # 生成正弦波

# 添加一些噪声
noise = np.random.normal(0, 0.1, signal.shape)
noisy_signal = signal + noise

# 计算FFT
fft_result = np.fft.fft(noisy_signal)
fft_magnitude = np.abs(fft_result)  # 计算幅度
fft_phase = np.angle(fft_result)  # 计算相位
fft_frequency = np.fft.fftfreq(len(t), t[1] - t[0])  # 计算频率轴

# 绘制原始信号
plt.figure(figsize=(12, 9))
plt.subplot(3, 1, 1)
plt.plot(t, noisy_signal, label='Noisy Signal')
plt.title('Original Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()

# 绘制FFT幅度结果
plt.subplot(3, 1, 2)
plt.stem(fft_frequency[:len(fft_frequency)//2], fft_magnitude[:len(fft_magnitude)//2], 'b', markerfmt=" ", basefmt="-b")
plt.title('FFT Magnitude')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim([0, 10])  # 限制x轴范围以显示正频率部分

# 绘制FFT相位结果
plt.subplot(3, 1, 3)
plt.stem(fft_frequency[:len(fft_frequency)//2], fft_phase[:len(fft_phase)//2], 'r', markerfmt=" ", basefmt="-r")
plt.title('FFT Phase')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Phase (Radians)')
plt.xlim([0, 10])  # 限制x轴范围以显示正频率部分

plt.tight_layout()
plt.show()

如下图所示，FFT振幅结果显示频率的峰值出现在X轴上的5Hz位置，表明信号中存在一个频率为5Hz的周期性成分，即每秒钟波动5次。我们通过观察所构建的正弦波图像可知，该波形一秒内的确波动了5次。在FFT的幅度谱图中，峰值表示信号中特定频率成分的强度或能量。明显的峰值意味着在该频率点，信号具有较大的能量集中，即信号中该频率成分较为显著。

二、FFT应用

现在，我们知道了FFT可以将原始数据分离成振幅＋相位的形式，那么如何应用于信号分析?这里，我们做一个FFT进行图像体积压缩的实践。首先，导入必要的库：

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

其次，我们直接使用numpy封装好的FFT模块，为了便于理解，这里为每一个模块都构建一个函数：

def fft2(x):
    '''
    函数np.fft.fft2计算二维离散傅里叶变换，
    它适用于二维输入数组，比如图像数据（通常是灰度图或彩色图像的每个通道），
    这个函数将空间域（或时间域）的信号转换到频率域，揭示了信号的频率成分，
    对于图像而言，fft2可以帮助识别图像中的频率信息，比如边缘和纹理等。
    '''
    return np.fft.fft2(x)

def ifft2(x):
    '''
    函数np.fft.ifft2是fft2的逆变换，它计算二维离散逆傅里叶变换，
    这个函数将频率域的信号转换回空间域，用于从频率信息重建原始信号，
    在图像处理中，ifft2可以用来从其频域表示中恢复图像。
    '''
    return np.fft.ifft2(x)

def fftshift(x):
    '''
    函数np.fft.fftshift将零频率分量移动到数组的中心，
    在FFT变换后，零频率分量（直流分量）通常位于数组的开始位置，
    fftshift通过对数组进行重新排序，将这个零频率分量移动到中心位置，
    这样做的目的是为了更直观地观察频谱的分布并处理，特别是在进行频域滤波时。
    '''
    return np.fft.fftshift(x)

def ifftshift(x):
    '''
    函数np.fft.ifftshift是fftshift的逆操作，它将经过fftshift处理的数组恢复到原始的顺序，
    在进行了频域处理之后，比如滤波，使用ifftshift将数组恢复到零频率分量位于数组开始位置的状态，
    以便进行逆FFT变换。
    '''
    return np.fft.ifftshift(x)

接着，我们定义图像压缩的主函数，下面举的是一个很简单的策略，代码中已经给出了详细的注释。

def compress_image(img_array, compression_ratio):
    # 应用二维快速傅里叶变换
    f = fft2(img_array)
    fshift = fftshift(f)

    # 定义一个简单的低通滤波器
    rows, cols = img_array.shape[:2]  # 获取维度
    # 计算图像的中心行和列的坐标
    crow, ccol = rows//2 , cols//2
    # 创建一个与输入图像大小相同的遮罩数组 mask，初始值为零
    mask = np.zeros((rows, cols), np.uint8)
    # 根据压缩比率计算滤波器半径，这里取图像的最小维度的一半乘以压缩比率，以确保滤波器半径不会超出图像边界
    r = int((min(rows, cols) / 2) * compression_ratio)  
    center = [crow, ccol]
    # 使用 np.ogrid 创建两个开放网格 x 和 y，用于生成遮罩规则数组
    x, y = np.ogrid[:rows, :cols]
    # 计算一个圆形区域，该区域内的点到中心点的距离小于或等于半径 r，结果存储在布尔数组 mask_area 中
    mask_area = (x - center[0]) ** 2 + (y - center[1]) ** 2 <= r*r
    # 将圆形区域内的遮罩值设置为1，表示这些区域将被保留
    mask[mask_area] = 1

    # 应用低通滤波器，将 FFT 结果 fshift 与遮罩 mask 相乘，实现低通滤波效果
    fshift = fshift * mask

    # 应用逆FFT
    ishift = ifftshift(fshift)
    img_back = ifft2(ishift)
    # 取 IFFT 结果 img_back 的绝对值，因为原始图像是复数。
    img_back = np.abs(img_back)

    return img_back

最后，在目标图片上应用我们定义好的压缩函数：

# 读取图像
# 这里会将图像转为灰度图
img = Image.open('1.jpg').convert('L')
img_array = np.array(img)
compressed_img_array = compress_image(img_array, 0.5)  # 设置压缩比率

# 显示结果
plt.figure(figsize=(12, 6))

plt.subplot(121), plt.imshow(img_array)
plt.title('Original Image'), plt.xticks([]), plt.yticks([])

plt.subplot(122), plt.imshow(compressed_img_array)
plt.title('Compressed Image'), plt.xticks([]), plt.yticks([])
plt.show()

Image.fromarray((compressed_img_array).astype(np.uint8)).save('compressed_image.jpg')

随便拿一张图片跑一下，查看保存后的图片大小可知，新图片的大小为原图的2/3左右。