- 一、实验目的和要求
- 理解傅里叶变换的思想与功能。
- 掌握傅里叶变换的Python实现方法：包括OpenCV工具包相关函数和NumPy工具包相关函数。
- 理解并掌握图像频域滤波的用途与实现方法。
二、实验内容和步骤
1. 编程分别使用OpenCV和NumPy工具包中的函数对图像peppers.bmp进行傅里叶变换，要求显示原图像、变换后图像的幅值谱、变换后图像的对数变换幅值谱、变换后图像的相位谱，以及重构后的图像。
2. 图像alphabet.jpg分别使用截止频率15,30,80,220进行理想低通滤波，观察不同的滤波效果。
3. 图像alphabet.jpg分别使用截止频率15,30,80,220进行高斯低通滤波，观察不同的滤波效果。
4. 图像lena.jpg分别使用截止频率15,30,80,220，阶数分别取2和5进行巴特沃斯高通滤波，观察不同的滤波效果。
1. 编程分别使用OpenCV和NumPy工具包中的函数对图像peppers.bmp进行傅里叶变换，要求显示原图像、变换后图像的幅值谱、变换后图像的对数变换幅值谱、变换后图像的相位谱，以及重构后的图像。
import numpy as np import cv2 from matplotlib import pyplot as plt plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签 plt.rcParams['axes.unicode_minus']=False #用来正常显示负 #Numpy函数实现傅里叶变换 img = cv2.imread(r"C:\Users\Administrator\Desktop\photo\peppers.bmp",0)#读取bmp文件 dft = np.fft.fft2(img)#对图像进行二维FFT magnitude0= np.abs(dft)#计算幅值谱 magnitude1= 20*np.log(1+magnitude0)# 可视化幅值谱（取对数并显示） phase_angle0=np.angle(dft) print(phase_angle0) img_back=np.fft.ifft2(dft) img_back1 = np.abs(img_back) plt.figure(figsize=(10,6)) plt.subplot(151),plt.imshow(img, cmap = 'gray') plt.title('原图像'), plt.axis('off') plt.subplot(152),plt.imshow(magnitude0, cmap = 'hsv') plt.title('幅值谱'), plt.axis('off') plt.subplot(153),plt.imshow(magnitude1, cmap = 'hsv') plt.title('对数变换后的幅值谱'), plt.axis('off') plt.subplot(154),plt.imshow(phase_angle0, cmap = 'hsv') plt.title('相位谱'), plt.axis('off') plt.subplot(155),plt.imshow(img_back1,cmap = 'hsv') plt.title('重构图像'),plt.axis('off') plt.show() plt.savefig(r"C:\Users\Administrator\Desktop\photo\sys4-1-2.jpg")
- dft = np.fft.fft2(img)：对输入的二维数组（通常是一个图像）img进行二维FFT，并将结果存储在dft中。
- img_back=np.fft.ifft2(dft)：dft 是原始图像经过二维快速傅里叶变换（FFT）得到的结果。为了从 dft 重构（或反变换）回原始图像，需要执行逆FFT（Inverse FFT）。在NumPy中，可以使用 np.fft.ifft2 函数来实现这一点。
- img_back1 = np.abs(img_back)：img_back 是通过逆FFT得到的复数数组。由于傅里叶变换涉及复数运算，因此逆FFT的结果也是复数。然而，对于图像来说，我们通常只对幅度感兴趣，因此通过取复数的绝对值（np.abs(img_back)）得到 img_back1，代表了重构图像的幅度信息。请注意，由于在进行FFT和IFFT时可能会引入一些数值误差，因此重构的图像可能不会与原始图像完全相同，但它们应该非常接近。
- phase_angle0=np.angle(dft)：np.angle(dft) 函数用于计算复数数组 dft 中每个元素的相位角（phase angle）。当 dft 是通过二维快速傅里叶变换（FFT）得到的结果时，np.angle(dft) 返回的是每个频率分量的相位信息。
运行结果：
2. 图像alphabet.jpg分别使用截止频率15,30,80,220进行理想低通滤波，观察不同的滤波效果。
import cv2 import numpy as np from math import * import random from matplotlib import pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负 img = cv2.imread(r"C:\Users\Administrator\Desktop\photo\alphabet.jpg", 0) f = np.fft.fft2(img)#对图像进行二维FFT fshift = np.fft.fftshift(f) #实现傅里叶频谱的低频移中 magnitude_spectrum0 = 20 * np.log(1 + np.abs(fshift)) plt.figure(figsize=(20, 8)) # 进行理想低通滤波 dd=[15,30,80,220] #截止频率 for w in range(4): r =dd[w] # 截止频率的设置 [m, n] = fshift.shape H = np.zeros((m, n), dtype=np.complex_)#初始化一个与fshift相同大小的复数矩阵H，用于存储滤波器的值。 # 接下来的双重循环创建了一个理想低通滤波器H。 # 如果某个频率分量（由i和j索引表示）的距离小于截止频率r， # 则对应的H值设置为1，否则为0。 for i in range(m): for j in range(n): d = sqrt((i - m / 2) * (i - m / 2) + (j - n / 2) * (j - n / 2)) if d < r: H[i, j] = 1 G = H * fshift#将滤波器H应用于频谱fshift，得到滤波后的频谱。 magnitude_spectrum1 = 20 * np.log(1 + np.abs(G)) # 理想低通滤波后的幅值谱 f1 = np.fft.ifftshift(G) img1 = abs(np.fft.ifft2(f1)) # 重构图像 plt.subplot(2,4,w+1), plt.imshow(magnitude_spectrum1, cmap='hsv') # 显示滤波后幅值谱 plt.title('截止频率'+str(dd[w])+'ILPF滤波后幅值谱'), plt.axis("off") plt.subplot(2,4,w+5), plt.imshow(img1, cmap='hsv') # 显示重构图像 plt.title('截止频率'+str(dd[w])+'ILPF滤波后重构图像'), plt.axis("off") plt.show() plt.savefig(r"C:\Users\Administrator\Desktop\photo\sys4-2.jpg")- fshift = np.fft.fftshift(f) ：
- 在NumPy中，np.fft.fftshift函数用于将零频率分量（DC分量）移动到频谱的中心。在未经移动的FFT结果中，零频率分量通常位于频谱的第一个位置，而在使用fftshift之后，它会移动到频谱的中间位置。这对于许多信号和图像处理任务来说是非常有用的，因为它允许我们更容易地分析和可视化频谱。
- 在图像处理中，通常希望低频（接近DC）分量位于频谱的中心，因为低频分量通常对应于图像的整体亮度或颜色，而高频分量对应于图像的细节和边缘。通过将低频分量移动到频谱中心，我们可以更容易地观察到这些特征。
- magnitude_spectrum0 = 20 * np.log(1 + np.abs(fshift))：
这里做了几件事：
- np.abs(fshift)：计算傅里叶变换结果 fshift 的绝对值，这给出了每个频率分量的幅值。
- 1 + np.abs(fshift)：在取对数之前，加1是为了避免对0取对数，因为对0取对数在数学上是没有定义的。这样做可以确保即使幅值为0，计算也不会出错。
- np.log(1 + np.abs(fshift))：取对数，通常是自然对数（np.log 默认计算的是自然对数，即以 e 为底的对数）。对数变换将幅值谱的范围压缩到一个更易于可视化的区间。
- 20 * ...：这是一个常见的缩放因子，用于将对数变换后的结果调整到一个合适的显示范围。这个因子是基于经验选择的，以便在图像显示时能够清楚地看到高频和低频成分
最终，magnitude_spectrum0 将包含对数变换后的幅值谱，它更适合于可视化，因为高频和低频成分在图像上都会有一个合理的亮度表示，而不是被高动态范围所掩盖。
运行结果：
3. 图像alphabet.jpg分别使用截止频率15,30,80,220进行高斯低通滤波，观察不同的滤波效果。
import cv2 as cv import numpy as np from matplotlib import pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负 img = cv.imread(r"C:\Users\Administrator\Desktop\photo\alphabet.jpg", 0) f = np.fft.fft2(img) fshift0 = np.fft.fftshift(f) dd = [15, 30, 80, 220] rows, cols = img.shape crow, ccol = rows // 2, cols // 2 plt.figure(figsize=(20, 8)) #以下这段循环代码的主要目的是对给定的截止频率列表dd中的每一个值， # 执行二维高通滤波，并显示滤波后的幅值谱和重构的图像。 for w in range(4): d0 = dd[w] ** 2 H = np.zeros((rows, cols)) for u in range(rows): for v in range(cols): d = (u - crow) ** 2 + (v - ccol) ** 2#d是像素到图像中心的距离的平方 H[u, v] = np.exp(-d / (2 * d0)) fshift = fshift0 * H#对原始频域数据fshift0应用高通滤波器H。 f_ishift = np.fft.ifftshift(fshift)#将滤波后的频域数据fshift进行逆移位，以准备进行逆傅里叶变换。 temp = np.fft.ifft2(f_ishift)#对逆移位后的频域数据进行二维逆傅里叶变换，得到空间域的图像。 img_back= np.abs(temp)#取逆傅里叶变换结果的绝对值，得到重构的图像。 plt.subplot(2,4,w+1), plt.imshow(20 * np.log(1 + np.abs(fshift)), cmap='hsv') plt.title('截止频率'+str(dd[w])+'GLPF滤波后幅值谱'), plt.axis("off") plt.subplot(2,4,w+5), plt.imshow(img_back, cmap='hsv') plt.title('截止频率'+str(dd[w])+'GLPF滤波后重构图像'), plt.axis("off") plt.show() plt.savefig(r"C:\Users\Administrator\Desktop\photo\sys4-3.jpg")
- H[u, v] = np.exp(-d / (2 * d0))：
这个表达式可以用于计算一个二维高斯核（或者称为高斯滤波器），该核可以用于图像处理中的平滑或模糊操作。在这种情况下，u 和 v 通常表示图像中的像素坐标，而 d 可以是像素之间的某种距离度量（例如欧几里得距离）。d0 则控制了高斯核的宽度或平滑程度。
运行结果：
4. 图像lena.jpg分别使用截止频率15,30,80,220，阶数分别取2和5进行巴特沃斯高通滤波，观察不同的滤波效果。
import cv2 as cv import numpy as np import math from matplotlib import pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负 img = cv.imread(r"C:\Users\Administrator\Desktop\photo\lena.jpg", 0) f = np.fft.fft2(img) fshift0 = np.fft.fftshift(f) dd = [15, 30, 80, 220] rows, cols = img.shape crow, ccol = rows // 2, cols // 2 plt.figure(figsize=(20, 8)) for w in range(4): d0 = dd[w] ** 2 H = np.zeros((rows, cols)) for u in range(rows): for v in range(cols): d = (u - crow) ** 2 + (v - ccol) ** 2 H[u, v] =1.0/(1.0+(d0/(d+0.01))**5) fshift = fshift0 * H f_ishift = np.fft.ifftshift(fshift) temp = np.fft.ifft2(f_ishift) img_back= np.abs(temp) plt.subplot(2,4,w+1), plt.imshow(20 * np.log(1 + np.abs(fshift)), cmap='hsv') plt.title('截止频率'+str(dd[w])+'BHPF滤波后幅值谱'), plt.axis("off") plt.subplot(2,4,w+5), plt.imshow(img_back, cmap='hsv') plt.title('截止频率'+str(dd[w])+'BHPF滤波后重构图像'), plt.axis("off") plt.show() plt.savefig(r"C:\Users\Administrator\Desktop\photo\sys4-4.jpg")
- H[u, v] =1.0/(1.0+(d0/(d+0.01))**5)：
这个表达式可以看作是一个平滑过渡函数或者称为窗口函数，用于控制某种影响的范围或强度，根据 d 和 d0 的关系来调整 H[u, v] 的值。具体来说：
- 当 d 远大于 d0 时，d0 / (d + 0.01) 将接近于 0，因此 H[u, v] 将接近于 1.0。
- 当 d 远小于 d0 时，d0 / (d + 0.01) 将非常大，因此 H[u, v] 将接近于 0.0。
- 当 d 接近于 d0 时，H[u, v] 的值将在 0.0 和 1.0 之间变化，具体取决于 d 和 d0 的具体数值。
这里的 +0.01 是一个小的常数，用于防止分母为 0 的情况，确保表达式在所有情况下都是定义良好的。
这个函数可以用于多种应用，例如图像处理中的权重分配、机器学习中的核函数、空间插值等。根据上下文，u 和 v 可以代表像素坐标、数据点坐标或其他类型的空间坐标。d 可以是这些坐标之间的距离或相似度度量，而 d0 是一个控制函数形状的参数。
例如，在图像处理中，H[u, v] 可以用于定义一个平滑的权重函数，用于在像素 u, v 周围进行加权平均。这样，距离中心较远的像素将有较小的权重，而距离较近的像素将有较大的权重。这种权重函数可以用于实现图像模糊、特征提取或其他类型的空间滤波操作。
运行结果：
三、补充知识点
- 相位角（Phase Angle）是指复数在复平面上与正实轴之间的角度，它描述了信号在不同频率上的相位延迟。
- 相位谱（Phase Spectrum）则是频率域中每个频率分量的相位角的表示。在信号处理中，相位谱与幅值谱（Magnitude Spectrum）一样重要，因为它们共同决定了信号在频率域中的完整表示。
- 快速傅里叶变换（Fast Fourier Transform，FFT）是离散傅里叶变换（DFT）的高效算法，用于计算DFT及其逆变换。FFT算法通过将DFT分解为若干个小DFT的方式，显著降低了计算复杂度，使得大规模数据处理变得可行。FFT在计算机科学、数字信号处理、图像处理、通信等领域有广泛应用。
- cmap全称为"colormap parameter"，中文为颜色映射参数，用于定义颜色映射的方式。
使用HSV色彩映射来显示相位谱，因为相位通常表示为一个角度，HSV色彩空间中的色调（Hue）可以很好地表示这种周期性变化。
扫一扫
22万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
