《数字图像处理》第 4 章 - 频域滤波

pip install numpy opencv-python matplotlib scipy

import numpy as np

# 定义复数（三种正确方式，任选其一即可）
# 方式1：直接用Python原生复数（推荐，最简单）
z1 = 3 + 4j  # 3+4j
z2 = 1 + 2j  # 1+2j

# 方式2：np.complex_ 传入字符串（兼容旧版本NumPy）
# z1 = np.complex_("3+4j")
# z2 = np.complex_("1+2j")

# 方式3：np.array 构造复数数组（适合批量定义）
# z1 = np.array([3, 4], dtype=np.complex128)
# z2 = np.array([1, 2], dtype=np.complex128)

# 基本运算
z_add = z1 + z2       # 加法
z_mul = z1 * z2       # 乘法
z_abs = np.abs(z1)    # 模长
z_angle = np.angle(z1)# 辐角（弧度）

# 输出结果
print(f"z1 = {z1}, z2 = {z2}")
print(f"z1 + z2 = {z_add}")
print(f"z1 * z2 = {z_mul}")
print(f"|z1| = {z_abs}, 辐角 = {z_angle:.2f} rad")

z1 = (3+4j), z2 = (1+2j)
z1 + z2 = (4+6j)
z1 * z2 = (-5+10j)
|z1| = 5.0, 辐角 = 0.93 rad

import numpy as np
import matplotlib.pyplot as plt

# 设置中文字体（解决matplotlib中文显示问题）
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义方波函数
def square_wave(t, T=2*np.pi):
    return np.where(np.mod(t, T) < T/2, 1, -1)

# 傅里叶级数逼近
def fourier_series(t, n_max, T=2*np.pi):
    omega0 = 2*np.pi / T
    result = 0
    for n in range(1, n_max+1, 2):  # 只保留奇次谐波
        an = 0
        bn = 4/(n*np.pi)
        result += an*np.cos(n*omega0*t) + bn*np.sin(n*omega0*t)
    return result

# 生成数据
t = np.linspace(-2*np.pi, 2*np.pi, 1000)
y_original = square_wave(t)

# 不同项数的逼近
y_1 = fourier_series(t, 1)
y_5 = fourier_series(t, 5)
y_20 = fourier_series(t, 20)

# 可视化对比
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.plot(t, y_original, label='原始方波')
plt.title('原始方波')
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(t, y_1, label='1项逼近', color='orange')
plt.title('1项傅里叶级数逼近')
plt.grid(True)

plt.subplot(2, 2, 3)
plt.plot(t, y_5, label='5项逼近', color='red')
plt.title('5项傅里叶级数逼近')
plt.grid(True)

plt.subplot(2, 2, 4)
plt.plot(t, y_20, label='20项逼近', color='green')
plt.title('20项傅里叶级数逼近')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 模拟冲激函数（离散近似）
def delta_func(t, t0=0, eps=0.01):
    return np.where(np.abs(t - t0) < eps, 1/eps, 0)

# 定义测试函数
def f(t):
    return np.sin(t) + 0.5*t

# 生成数据
t = np.linspace(-5, 5, 1000)
delta = delta_func(t, t0=1)  # t0=1处的冲激
f_t = f(t)

# 验证筛分性质：积分f(t)*delta(t-1) ≈ f(1)
integral = np.trapz(f_t * delta, t)
f_1 = f(1)

# 可视化
plt.figure(figsize=(10, 6))
plt.plot(t, f_t, label='f(t) = sin(t) + 0.5t', color='blue')
plt.plot(t, delta, label='δ(t-1)', color='red', alpha=0.7)
plt.scatter(1, f_1, color='green', s=100, label=f'f(1) = {f_1:.2f}')
plt.title(f'冲激函数的筛分性质：积分结果 = {integral:.2f} ≈ f(1)')
plt.xlabel('t')
plt.ylabel('幅值')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义两个函数
t = np.linspace(-5, 5, 1000)
f = np.exp(-t**2)  # 高斯函数
g = np.where((t >= -1) & (t <= 1), 1, 0)  # 矩形窗

# 计算卷积
conv = convolve(f, g, mode='same') / len(t)  # 归一化

# 可视化
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.plot(t, f, label='f(t) = e^(-t²)')
plt.title('函数f(t)')
plt.grid(True)

plt.subplot(1, 3, 2)
plt.plot(t, g, label='g(t) 矩形窗')
plt.title('函数g(t)')
plt.grid(True)

plt.subplot(1, 3, 3)
plt.plot(t, conv, label='f*g', color='red')
plt.title('卷积结果 (f*g)(t)')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义原始信号：10Hz + 20Hz 正弦波
def signal(t):
    return np.sin(2*np.pi*10*t) + 0.5*np.sin(2*np.pi*20*t)

# 生成连续信号
t_continuous = np.linspace(0, 1, 1000)
y_continuous = signal(t_continuous)

# 采样1：满足奈奎斯特（采样频率50Hz > 2*20Hz）
fs1 = 50
t1 = np.linspace(0, 1, fs1)
y1 = signal(t1)

# 采样2：不满足奈奎斯特（采样频率30Hz < 2*20Hz）
fs2 = 30
t2 = np.linspace(0, 1, fs2)
y2 = signal(t2)

# 可视化
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(t_continuous, y_continuous, label='原始连续信号', color='blue')
plt.scatter(t1, y1, color='red', label=f'采样频率{fs1}Hz（无混叠）')
plt.title('满足奈奎斯特定理：无混叠')
plt.xlabel('时间 (s)')
plt.ylabel('幅值')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(t_continuous, y_continuous, label='原始连续信号', color='blue')
plt.scatter(t2, y2, color='orange', label=f'采样频率{fs2}Hz（混叠）')
plt.title('不满足奈奎斯特定理：混叠失真')
plt.xlabel('时间 (s)')
plt.ylabel('幅值')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 生成离散信号：5Hz正弦波，采样频率50Hz，长度100
fs = 50
N = 100
t = np.arange(N) / fs
f_n = np.sin(2*np.pi*5*t)

# 计算DFT
F_k = np.fft.fft(f_n)
# 频率轴
freq = np.fft.fftfreq(N, 1/fs)
# 幅值谱（归一化）
amp = np.abs(F_k) / N

# 可视化
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t, f_n)
plt.title('原始离散信号（5Hz正弦波）')
plt.xlabel('时间 (s)')
plt.ylabel('幅值')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.stem(freq, amp, basefmt='b-')
plt.title('DFT幅值谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅值')
plt.xlim(0, fs/2)  # 只显示正频率
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import cv2
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 读取图像（灰度图）
img = cv2.imread('lena.jpg', cv2.IMREAD_GRAYSCALE)
if img is None:
    # 若读取失败，生成测试图像
    img = np.zeros((256, 256), dtype=np.uint8)
    img[100:150, 100:150] = 255  # 白色方块

# 计算二维DFT
dft = np.fft.fft2(img)
# 将低频移到中心
dft_shift = np.fft.fftshift(dft)
# 幅值谱（对数缩放，便于可视化）
magnitude_spectrum = 20 * np.log(np.abs(dft_shift))

# 逆DFT
dft_shift_back = np.fft.ifftshift(dft_shift)
img_back = np.fft.ifft2(dft_shift_back)
img_back = np.abs(img_back).astype(np.uint8)

# 可视化
plt.figure(figsize=(12, 6))
plt.subplot(1, 3, 1)
plt.imshow(img, cmap='gray')
plt.title('原始图像')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('DFT幅值谱（中心对齐）')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.imshow(img_back, cmap='gray')
plt.title('逆DFT恢复图像')
plt.axis('off')

plt.tight_layout()
plt.show()

import numpy as np
import cv2
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 读取图像
img = cv2.imread('lena.jpg', cv2.IMREAD_GRAYSCALE)
if img is None:
    img = np.zeros((256, 256), dtype=np.uint8)
    img[100:150, 100:150] = 255

# 图像平移（右移50，下移50）
M, N = img.shape
dx, dy = 50, 50
img_shift = np.roll(np.roll(img, dx, axis=1), dy, axis=0)

# 计算DFT
dft_original = np.fft.fftshift(np.fft.fft2(img))
dft_shift = np.fft.fftshift(np.fft.fft2(img_shift))

# 幅值谱
mag_original = 20 * np.log(np.abs(dft_original))
mag_shift = 20 * np.log(np.abs(dft_shift))

# 可视化
plt.figure(figsize=(12, 6))
plt.subplot(2, 2, 1)
plt.imshow(img, cmap='gray')
plt.title('原始图像')
plt.axis('off')

plt.subplot(2, 2, 2)
plt.imshow(img_shift, cmap='gray')
plt.title('平移后图像')
plt.axis('off')

plt.subplot(2, 2, 3)
plt.imshow(mag_original, cmap='gray')
plt.title('原始图像DFT幅值谱')
plt.axis('off')

plt.subplot(2, 2, 4)
plt.imshow(mag_shift, cmap='gray')
plt.title('平移图像DFT幅值谱')
plt.axis('off')

plt.tight_layout()
plt.show()

import numpy as np
import cv2
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

img = cv2.imread('lena.jpg', cv2.IMREAD_GRAYSCALE)
if img is None:
    img = np.zeros((256, 256), dtype=np.uint8)
    img[100:150, 100:150] = 255

# 计算DFT并移到中心
dft = np.fft.fft2(img)
dft_shift = np.fft.fftshift(dft)

# 幅值谱和相位谱
magnitude = 20 * np.log(np.abs(dft_shift))
phase = np.angle(dft_shift)

# 可视化
plt.figure(figsize=(12, 6))
plt.subplot(1, 3, 1)
plt.imshow(img, cmap='gray')
plt.title('原始图像')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.imshow(magnitude, cmap='gray')
plt.title('幅值谱')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.imshow(phase, cmap='gray')
plt.title('相位谱')
plt.axis('off')

plt.tight_layout()
plt.show()

import numpy as np
import cv2
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 1. 读取并预处理图像
img = cv2.imread('../picture/AALi.jpg', cv2.IMREAD_GRAYSCALE)
if img is None:
    # 生成测试图像（添加噪声）
    img = np.zeros((256, 256), dtype=np.uint8)
    img[80:180, 80:180] = 255
    # 添加高斯噪声
    noise = np.random.normal(0, 30, img.shape).astype(np.float32)
    img = np.clip(img + noise, 0, 255).astype(np.uint8)

M, N = img.shape  # M=高度(行), N=宽度(列)
print(f"图像形状：M={M}, N={N}")  # 打印维度，方便排查

# 2. 计算DFT并中心对齐
dft = np.fft.fft2(img)
dft_shift = np.fft.fftshift(dft)

# 3. 生成频域网格（关键修正：先v后u，匹配图像行列）
# 注意：meshgrid第一个参数对应列(v)，第二个对应行(u)
v = np.arange(N)  # 列方向（宽度）
u = np.arange(M)  # 行方向（高度）
u_grid, v_grid = np.meshgrid(u, v, indexing='ij')  # indexing='ij' 确保(u,v)对应(行,列)

# 中心坐标
u0, v0 = M // 2, N // 2
# 距离D(u,v)：计算每个频域点到中心的欧式距离
D = np.sqrt((u_grid - u0) ** 2 + (v_grid - v0) ** 2)
print(f"距离矩阵形状：{D.shape}，DFT移位后形状：{dft_shift.shape}")  # 验证维度匹配

# 4. 设计滤波器（截止频率D0=30）
D0 = 30
n = 2  # 巴特沃斯阶数

# 理想低通
H_ilpf = np.where(D <= D0, 1, 0)
# 巴特沃斯低通
H_blpf = 1 / (1 + (D / D0) ** (2 * n))
# 高斯低通
H_glpf = np.exp(-(D ** 2) / (2 * D0 ** 2))

# 5. 频域滤波（现在维度匹配，可正常相乘）
g_shift_ilpf = dft_shift * H_ilpf
g_shift_blpf = dft_shift * H_blpf
g_shift_glpf = dft_shift * H_glpf

# 6. 逆DFT（添加小技巧：取实部后归一化，避免数值误差）
# 理想低通
g_ilpf = np.fft.ifftshift(g_shift_ilpf)
g_ilpf = np.fft.ifft2(g_ilpf)
g_ilpf = np.real(g_ilpf)  # 优先取实部，避免虚部残留
g_ilpf = np.clip(g_ilpf, 0, 255).astype(np.uint8)  # 归一化到0-255

# 巴特沃斯低通
g_blpf = np.fft.ifftshift(g_shift_blpf)
g_blpf = np.fft.ifft2(g_blpf)
g_blpf = np.real(g_blpf)
g_blpf = np.clip(g_blpf, 0, 255).astype(np.uint8)

# 高斯低通
g_glpf = np.fft.ifftshift(g_shift_glpf)
g_glpf = np.fft.ifft2(g_glpf)
g_glpf = np.real(g_glpf)
g_glpf = np.clip(g_glpf, 0, 255).astype(np.uint8)

# 7. 可视化（补充第8个子图，避免布局警告）
plt.figure(figsize=(15, 10))
# 原始图像
plt.subplot(2, 4, 1)
plt.imshow(img, cmap='gray')
plt.title('原始图像（含噪声）')
plt.axis('off')
# 理想低通
plt.subplot(2, 4, 2)
plt.imshow(g_ilpf, cmap='gray')
plt.title('理想低通滤波（振铃效应明显）')
plt.axis('off')
# 巴特沃斯低通
plt.subplot(2, 4, 3)
plt.imshow(g_blpf, cmap='gray')
plt.title('巴特沃斯低通（轻微振铃）')
plt.axis('off')
# 高斯低通
plt.subplot(2, 4, 4)
plt.imshow(g_glpf, cmap='gray')
plt.title('高斯低通（无振铃）')
plt.axis('off')

# 滤波器可视化
plt.subplot(2, 4, 5)
plt.imshow(H_ilpf, cmap='gray')
plt.title('理想低通滤波器')
plt.axis('off')
plt.subplot(2, 4, 6)
plt.imshow(H_blpf, cmap='gray')
plt.title('巴特沃斯低通滤波器')
plt.axis('off')
plt.subplot(2, 4, 7)
plt.imshow(H_glpf, cmap='gray')
plt.title('高斯低通滤波器')
plt.axis('off')
# 补充第8个子图（避免布局警告）
plt.subplot(2, 4, 8)
plt.text(0.5, 0.5, '低通滤波对比', ha='center', va='center', fontsize=12)
plt.axis('off')

plt.tight_layout()
plt.show()

import numpy as np
import cv2
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 1. 读取图像
img = cv2.imread('../picture/TianHuoSanXuanBian.jpg', cv2.IMREAD_GRAYSCALE)
if img is None:
    # 生成测试图像（低对比度）
    img = np.zeros((256, 256), dtype=np.uint8)
    img[80:180, 80:180] = 255
    # 降低对比度
    img = (img * 0.5).astype(np.uint8)

M, N = img.shape  # M=高度(行), N=宽度(列)
print(f"图像形状：M={M}, N={N}")  # 打印维度，方便排查

# 2. 计算DFT并中心对齐
dft = np.fft.fft2(img.astype(np.float32))
dft_shift = np.fft.fftshift(dft)

# 3. 生成频域网格（关键修正：匹配图像行列维度）
# 注意：先定义列(v)、后定义行(u)，并使用indexing='ij'确保维度匹配
v = np.arange(N)  # 列方向（宽度）
u = np.arange(M)  # 行方向（高度）
u_grid, v_grid = np.meshgrid(u, v, indexing='ij')  # 生成(M,N)形状的网格

# 中心坐标
u0, v0 = M // 2, N // 2
# 距离D(u,v)：添加极小值eps避免除零
eps = 1e-8  # 防止D=0导致除零
D = np.sqrt((u_grid - u0) ** 2 + (v_grid - v0) ** 2) + eps
print(f"距离矩阵形状：{D.shape}，DFT移位后形状：{dft_shift.shape}")  # 验证维度匹配

# 4. 设计高通滤波器（截止频率D0=30）
D0 = 30
n = 2    # 巴特沃斯阶数

# 理想高通
H_ihpf = np.where(D <= D0, 0, 1)

# 巴特沃斯高通（修正除零问题）
H_bhpf = 1 / (1 + (D0 / D) ** (2 * n))

# 高斯高通
H_ghpf = 1 - np.exp(-(D ** 2) / (2 * D0 ** 2))

# 高频强调滤波（加入直流分量，提升亮度）
H_hfe = 0.5 + 0.7 * H_ghpf

# 5. 同态滤波（增强对比度）
gamma_h = 2.0  # 高频增益（增强细节）
gamma_l = 0.5  # 低频增益（抑制背景）
c = 1
D1 = 30
H_homomorphic = (gamma_h - gamma_l) * (1 - np.exp(-c * (D ** 2) / (D1 ** 2))) + gamma_l

# 6. 滤波计算（统一优化：取实部+归一化，避免虚部误差）
# 理想高通
g_ihpf = np.fft.ifftshift(dft_shift * H_ihpf)
g_ihpf = np.fft.ifft2(g_ihpf)
g_ihpf = np.real(g_ihpf)  # 优先取实部，避免虚部残留
g_ihpf = np.clip(g_ihpf, 0, 255).astype(np.uint8)

# 巴特沃斯高通
g_bhpf = np.fft.ifftshift(dft_shift * H_bhpf)
g_bhpf = np.fft.ifft2(g_bhpf)
g_bhpf = np.real(g_bhpf)
g_bhpf = np.clip(g_bhpf, 0, 255).astype(np.uint8)

# 高频强调
g_hfe = np.fft.ifftshift(dft_shift * H_hfe)
g_hfe = np.fft.ifft2(g_hfe)
g_hfe = np.real(g_hfe)
g_hfe = np.clip(g_hfe, 0, 255).astype(np.uint8)

# 同态滤波（优化数值稳定性）
img_log = np.log1p(img.astype(np.float32))  # log(1+x)避免0值
dft_log = np.fft.fft2(img_log)
dft_log_shift = np.fft.fftshift(dft_log)
g_log_shift = dft_log_shift * H_homomorphic
g_log = np.fft.ifftshift(g_log_shift)
g_log = np.fft.ifft2(g_log)
g_homomorphic = np.expm1(np.real(g_log))  # exp(x)-1还原
g_homomorphic = np.clip(g_homomorphic, 0, 255).astype(np.uint8)

# 7. 可视化（补充第6个子图，避免布局警告）
plt.figure(figsize=(15, 10))
plt.subplot(2, 3, 1)
plt.imshow(img, cmap='gray')
plt.title('原始图像（低对比度）')
plt.axis('off')

plt.subplot(2, 3, 2)
plt.imshow(g_ihpf, cmap='gray')
plt.title('理想高通滤波（振铃+过暗）')
plt.axis('off')

plt.subplot(2, 3, 3)
plt.imshow(g_bhpf, cmap='gray')
plt.title('巴特沃斯高通滤波')
plt.axis('off')

plt.subplot(2, 3, 4)
plt.imshow(g_hfe, cmap='gray')
plt.title('高频强调滤波（亮度提升）')
plt.axis('off')

plt.subplot(2, 3, 5)
plt.imshow(g_homomorphic, cmap='gray')
plt.title('同态滤波（对比度增强）')
plt.axis('off')

# 补充第6个子图（显示滤波器示意图，提升实用性）
plt.subplot(2, 3, 6)
plt.imshow(H_hfe, cmap='gray')
plt.title('高频强调滤波器')
plt.axis('off')

plt.tight_layout()
plt.show()

import numpy as np
import cv2
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 1. 读取/生成图像
img = cv2.imread('../picture/GaoDa.png', cv2.IMREAD_GRAYSCALE)
if img is None:
    # 备用测试图像（确保维度是行×列）
    img = np.ones((256, 256), dtype=np.uint8) * 128
    img[80:180, 80:180] = 255

# 核心：明确图像的行(M)和列(N)
M, N = img.shape  # M=高度(行), N=宽度(列)
print(f"✅ 图像形状（行×列）：{img.shape}")

# 2. 生成与图像完全匹配的网格（行×列 = M×N）
# 正确逻辑：先列(N)、后行(M)，indexing='xy' → 最终网格形状M×N
x = np.arange(N)  # 列坐标（0~N-1）
y = np.arange(M)  # 行坐标（0~M-1）
x_grid, y_grid = np.meshgrid(x, y, indexing='xy')  # 结果：x_grid(M×N), y_grid(M×N)
print(f"✅ 网格形状（行×列）：x_grid={x_grid.shape}, y_grid={y_grid.shape}")

# 3. 添加周期性噪声（维度100%匹配）
freq = 8  # 噪声频率（越大条纹越密）
noise_amp = 30  # 噪声幅值（越大条纹越明显）
# 生成与图像同维度的噪声
noise = noise_amp * (
    np.sin(2 * np.pi * freq * x_grid / N) +  # 水平条纹（沿列方向）
    np.sin(2 * np.pi * freq * y_grid / M)    # 垂直条纹（沿行方向）
)
print(f"✅ 噪声形状（行×列）：{noise.shape}")
# 叠加噪声并归一化
img_noisy = np.clip(img + noise, 0, 255).astype(np.uint8)

# 4. 计算DFT并中心化
dft = np.fft.fft2(img_noisy.astype(np.float32))
dft_shift = np.fft.fftshift(dft)
# 幅值谱（避免log(0)）
magnitude = 20 * np.log(np.abs(dft_shift) + 1e-8)

# 5. 设计陷波滤波器（精准抑制噪声频率）
# 频域中心
cy, cx = M // 2, N // 2
# 噪声频率对应的频域偏移（根据实际效果调整）
offset = freq
# 陷波中心（4个方向的噪声频率点）
notch_centers = [
    (cy, cx + offset), (cy, cx - offset),
    (cy + offset, cx), (cy - offset, cx)
]
print(f"✅ 陷波中心坐标：{notch_centers}")

# 巴特沃斯陷波滤波器参数
D0 = 15  # 陷波半径（大图像适当增大）
n = 4    # 阶数
eps = 1e-8  # 防除零

# 初始化全通滤波器
H_notch = np.ones((M, N), dtype=np.float32)
for (yc, xc) in notch_centers:
    # 计算每个点到陷波中心的距离
    D = np.sqrt((y_grid - yc)**2 + (x_grid - xc)**2) + eps
    # 陷波阻带（抑制该中心周围的频率）
    H_notch *= 1 / (1 + (D0 / D)**(2*n))

# 6. 频域滤波 + 逆变换
g_shift = dft_shift * H_notch
g = np.fft.ifftshift(g_shift)
g = np.fft.ifft2(g)
# 取实部并归一化（避免虚部残留）
g_filtered = np.clip(np.real(g), 0, 255).astype(np.uint8)

# 7. 可视化（清晰对比）
plt.figure(figsize=(18, 10))

# 原始图像
plt.subplot(2, 3, 1)
plt.imshow(img, cmap='gray')
plt.title('1. 原始无噪声图像', fontsize=12)
plt.axis('off')

# 含噪声图像
plt.subplot(2, 3, 2)
plt.imshow(img_noisy, cmap='gray')
plt.title('2. 含周期性条纹噪声的图像', fontsize=12)
plt.axis('off')

# DFT幅值谱
plt.subplot(2, 3, 3)
plt.imshow(magnitude, cmap='gray')
plt.title('3. 噪声图像DFT幅值谱\n（亮点=噪声频率）', fontsize=12)
plt.axis('off')

# 陷波滤波器
plt.subplot(2, 3, 4)
plt.imshow(H_notch, cmap='gray')
plt.title('4. 陷波滤波器\n（黑色=噪声抑制区）', fontsize=12)
plt.axis('off')

# 滤波后图像
plt.subplot(2, 3, 5)
plt.imshow(g_filtered, cmap='gray')
plt.title('5. 陷波滤波后图像\n（噪声已去除）', fontsize=12)
plt.axis('off')

# 噪声残留
plt.subplot(2, 3, 6)
residual = np.abs(img_noisy - g_filtered).astype(np.uint8)
plt.imshow(residual, cmap='gray')
plt.title('6. 噪声残留（越黑残留越少）', fontsize=12)
plt.axis('off')

plt.tight_layout()
plt.show()