企业官网建设流程全解析

Python实现信号3倍抽取与2倍内插：从时域波形到频谱验证

引言

数字信号处理中的采样率转换技术是现代通信系统和音频处理的核心组件之一。想象一下，你正在开发一个智能音箱系统，需要处理来自不同来源的音频信号——有些是48kHz的高质量音乐，有些是16kHz的语音指令。如何高效地在不同采样率之间转换而不损失关键信息？这正是抽取和内插技术大显身手的场景。

本文将带你用Python从零实现信号的3倍抽取和2倍内插操作，并通过时域波形和频域频谱的双重视角验证这些操作的效果。不同于纯理论讲解，我们将聚焦于工程实践，使用NumPy、SciPy和Matplotlib构建完整的信号处理流程，包括：

• 信号生成与可视化
• 抗混叠滤波器设计
• 抽取与内插的Python实现
• 频谱分析与验证

1. 基础准备：构建测试信号

1.1 生成多频测试信号

我们首先创建一个包含多个频率成分的复合信号，这样可以更全面地观察采样率转换对频谱的影响：

import numpy as np import matplotlib.pyplot as plt from scipy import signal # 参数设置 fs_original = 6000 # 原始采样率(Hz) duration = 0.1 # 信号时长(s) t = np.arange(0, duration, 1/fs_original) # 生成包含三个频率成分的信号 freqs = [300, 800, 1500] # 信号频率成分(Hz) signal_components = [np.sin(2*np.pi*f*t) for f in freqs] x = np.sum(signal_components, axis=0) # 合成信号 # 可视化 plt.figure(figsize=(12, 4)) plt.plot(t[:200], x[:200]) # 显示前200个采样点 plt.title("原始信号时域波形 (部分)") plt.xlabel("时间(s)") plt.ylabel("幅度") plt.grid() plt.show()

1.2 信号频谱分析

在开始处理前，先观察原始信号的频谱特性：

def plot_spectrum(x, fs, title): n = len(x) freq = np.fft.fftfreq(n, 1/fs)[:n//2] spectrum = np.abs(np.fft.fft(x))[:n//2] * 2/n plt.figure(figsize=(12, 4)) plt.stem(freq, spectrum, 'b', markerfmt=" ", basefmt="-b") plt.title(title) plt.xlabel("频率(Hz)") plt.ylabel("幅度") plt.grid() plt.xlim(0, fs/2) plot_spectrum(x, fs_original, "原始信号频谱")

2. 信号抽取：理论与实践

2.1 抽取的数学原理

抽取（Decimation）是指按整数因子D降低信号采样率的过程，包含两个关键步骤：

1. 抗混叠滤波：防止抽取后出现频谱混叠
2. 采样率降低：每隔D-1个样本保留一个样本

数学表达式为： $$ y[n] = x[nD] $$

其中D为抽取因子。在频域，抽取会导致频谱扩展和周期性重复：

$$ Y(e^{j\omega}) = \frac{1}{D}\sum_{k=0}^{D-1}X\left(e^{j(\omega-2\pi k)/D}\right) $$

2.2 Python实现3倍抽取

D = 3 # 抽取因子 fs_decimated = fs_original // D # 抽取后采样率 # 设计抗混叠滤波器 nyq = fs_original / 2 cutoff = fs_decimated / 2 # 新奈奎斯特频率 ripple = 0.5 # 通带波纹(dB) attenuation = 60 # 阻带衰减(dB) # 计算滤波器阶数和系数 order, normal_cutoff = signal.buttord(cutoff/nyq, 1.2*cutoff/nyq, ripple, attenuation) b, a = signal.butter(order, normal_cutoff, btype='low') # 应用滤波器 x_filtered = signal.lfilter(b, a, x) # 执行抽取 x_decimated = x_filtered[::D] # 绘制结果 t_decimated = t[::D] plt.figure(figsize=(12, 8)) plt.subplot(3, 1, 1) plt.plot(t[:200], x[:200], label='原始信号') plt.plot(t[:200], x_filtered[:200], 'r', label='滤波后信号') plt.title("抗混叠滤波效果对比") plt.legend() plt.subplot(3, 1, 2) plt.plot(t[:200], x[:200], 'b', label='原始信号') plt.stem(t_decimated[:70], x_decimated[:70], 'r', markerfmt="ro", basefmt=" ", label='抽取后信号') plt.title("3倍抽取时域效果") plt.legend() plt.subplot(3, 1, 3) plot_spectrum(x_decimated, fs_decimated, "抽取后信号频谱") plt.tight_layout() plt.show()

2.3 关键参数对比

3. 信号内插：从理论到实现

3.1 内插的数学原理

内插（Interpolation）是按整数因子I增加采样率的过程，包含三个主要步骤：

1. 零值插入：在原始样本间插入I-1个零
2. 镜像滤波：去除由零插入引入的高频镜像
3. 采样率提升：最终采样率为I×fs

数学表达式为： $$ v[n] = \begin{cases} x[n/I] & \text{当 } n/I \text{为整数} \ 0 & \text{其他情况} \end{cases} $$

频域效应表现为频谱压缩： $$ V(e^{j\omega}) = X(e^{j\omega I}) $$

3.2 Python实现2倍内插

I = 2 # 内插因子 fs_interpolated = fs_original * I # 内插后采样率 # 零值插入 x_zeros = np.zeros(len(x) * I) x_zeros[::I] = x # 每隔I-1个零插入一个原始样本 # 设计镜像滤波器 cutoff_interp = fs_original / 2 # 需要保留的原始带宽 nyq_interp = fs_interpolated / 2 order_interp, normal_cutoff_interp = signal.buttord(cutoff_interp/nyq_interp, 1.2*cutoff_interp/nyq_interp, ripple, attenuation) b_interp, a_interp = signal.butter(order_interp, normal_cutoff_interp, btype='low') # 应用滤波器 x_interpolated = signal.lfilter(b_interp, a_interp, x_zeros) # 绘制结果 t_interpolated = np.arange(0, duration, 1/fs_interpolated) plt.figure(figsize=(12, 8)) plt.subplot(3, 1, 1) plt.plot(t[:100], x[:100], 'bo-', label='原始信号') plt.stem(t_interpolated[:200], x_zeros[:200], 'r', markerfmt="ro", basefmt=" ", label='零插入信号') plt.title("2倍内插零插入效果") plt.legend() plt.subplot(3, 1, 2) plt.plot(t[:100], x[:100], 'bo-', label='原始信号') plt.plot(t_interpolated[:200], x_interpolated[:200], 'g.-', label='内插后信号') plt.title("2倍内插最终效果") plt.legend() plt.subplot(3, 1, 3) plot_spectrum(x_interpolated, fs_interpolated, "内插后信号频谱") plt.tight_layout() plt.show()

3.3 内插过程中的频谱变化

1. 零插入阶段：在原始频谱基础上，会在fs_original/2到fs_interpolated/2之间产生镜像频谱
2. 滤波阶段：低通滤波器保留原始频带(0到fs_original/2)，消除镜像成分
3. 最终频谱：与原始信号频谱形状相同，但横轴范围扩展到fs_interpolated/2

4. 综合应用：抽取与内插联合操作

在实际系统中，经常需要组合使用抽取和内插来实现非整数倍的采样率转换。让我们实现一个先3倍抽取再2倍内插的流程：

# 先3倍抽取 x_decimated = x_filtered[::D] # 对抽取后信号进行2倍内插 x_zeros_combined = np.zeros(len(x_decimated) * I) x_zeros_combined[::I] = x_decimated # 使用相同的镜像滤波器 x_combined = signal.lfilter(b_interp, a_interp, x_zeros_combined) # 最终采样率 fs_combined = fs_decimated * I # 2000 * 2 = 4000 Hz # 结果可视化 plt.figure(figsize=(12, 8)) plt.subplot(2, 1, 1) plt.plot(t[:150], x[:150], 'b', label='原始信号') plt.plot(np.arange(0, duration, 1/fs_combined)[:300], x_combined[:300], 'g', label='处理后信号') plt.title("3倍抽取+2倍内插时域对比") plt.legend() plt.subplot(2, 1, 2) plot_spectrum(x, fs_original, "原始信号频谱") plot_spectrum(x_combined, fs_combined, "处理后信号频谱") plt.legend(['原始', '处理后']) plt.tight_layout() plt.show()

4.1 性能指标对比

通过计算信噪比(SNR)来量化处理过程中的信号质量变化：

def calculate_snr(original, processed, fs_orig, fs_proc): # 对齐信号长度 min_len = min(len(original), int(len(processed) * fs_orig / fs_proc)) original = original[:min_len] processed = processed[:int(min_len * fs_proc / fs_orig):int(fs_proc/fs_orig)] # 计算信号功率和噪声功率 signal_power = np.sum(original**2) noise_power = np.sum((original - processed)**2) return 10 * np.log10(signal_power / noise_power) snr = calculate_snr(x, x_combined, fs_original, fs_combined) print(f"处理前后信号SNR: {snr:.2f} dB")

5. 工程实践中的注意事项

5.1 滤波器设计考量

在实际应用中，滤波器设计需要权衡三个关键参数：

1. 过渡带宽：越窄则滤波器阶数越高，计算量越大
2. 阻带衰减：要求越高则滤波器阶数越高
3. 通带波纹：越小则信号失真越小，但实现难度增加

推荐使用scipy.signal.remez设计等波纹滤波器，可以获得最优的幅频特性：

# 等波纹滤波器设计示例 bands = [0, 0.4, 0.6, 1] # 归一化频率 desired = [1, 0] # 期望增益 weight = [1, 100] # 通带和阻带权重 numtaps = 45 # 滤波器阶数 b_remez = signal.remez(numtaps, bands, desired, weight=weight)

5.2 多级处理策略

对于大的采样率转换因子（如D=10），采用多级处理可以显著降低计算复杂度：

1. 多级抽取：将D分解为D1×D2×...×Dn，每级使用较小的抽取因子
2. 多级内插：同理分解I因子

例如，10倍抽取可以分解为2×5：

# 两级抽取示例 D1, D2 = 2, 5 # 第一级：2倍抽取 cutoff1 = fs_original / (2 * D1) b1, a1 = signal.butter(8, cutoff1/nyq) x_dec1 = signal.lfilter(b1, a1, x)[::D1] # 第二级：5倍抽取 fs_dec1 = fs_original // D1 cutoff2 = fs_dec1 / (2 * D2) nyq_dec1 = fs_dec1 / 2 b2, a2 = signal.butter(8, cutoff2/nyq_dec1) x_dec2 = signal.lfilter(b2, a2, x_dec1)[::D2]

5.3 实时处理优化

在实时系统中，可以采用以下优化策略：

1. 多相滤波：减少计算量，特别适合内插操作
2. 环形缓冲区：高效处理连续数据流
3. 定点运算：在嵌入式系统中提高效率

多相滤波实现示例：

def polyphase_interp(x, I, b): """多相内插实现""" n_poly = I poly_filters = [b[i::n_poly] for i in range(n_poly)] # 上采样 x_up = np.zeros(len(x) * I) x_up[::I] = x # 多相滤波 y = np.zeros(len(x_up)) for i in range(n_poly): y[i::I] = signal.lfilter(poly_filters[i], 1, x_up[i::I]) return y

6. 验证与调试技巧

6.1 时频联合分析

开发过程中，建议同时观察时域和频域变化：

def analyze_signal(x, fs, title): plt.figure(figsize=(12, 5)) # 时域 plt.subplot(1, 2, 1) plt.plot(np.arange(len(x))/fs, x) plt.title(f"{title} - 时域") plt.xlabel("时间(s)") # 频域 plt.subplot(1, 2, 2) plot_spectrum(x, fs, f"{title} - 频域") plt.tight_layout() plt.show() analyze_signal(x, fs_original, "原始信号") analyze_signal(x_decimated, fs_decimated, "抽取后信号") analyze_signal(x_interpolated, fs_interpolated, "内插后信号")

6.2 常见问题排查

1. 频谱混叠：检查抗混叠滤波器的截止频率是否设置正确
2. 幅度失真：确认滤波器通带波纹是否过大
3. 相位失真：考虑使用零相位滤波(filtfilt代替lfilter)
4. 计算延迟：多级处理时注意各阶段的延迟累积

零相位滤波示例：

# 使用零相位滤波减少相位失真 x_filtered_zero = signal.filtfilt(b, a, x)

7. 扩展应用与性能评估

7.1 音频重采样案例

将44.1kHz音频转换为48kHz的典型应用：

def audio_resample(audio, fs_in, fs_out): """音频重采样函数""" # 计算有理数因子 L/M gcd_val = np.gcd(fs_in, fs_out) L = fs_out // gcd_val M = fs_in // gcd_val # 先L倍内插 x_up = np.zeros(len(audio) * L) x_up[::L] = audio # 设计抗镜像/抗混叠组合滤波器 cutoff = min(fs_in, fs_out) / 2 nyq = max(fs_in, fs_out) / 2 b = signal.remez(101, [0, 0.9*cutoff, 1.1*cutoff, nyq], [1, 0], fs=nyq*2) # 滤波 x_filt = signal.lfilter(b, 1, x_up) # M倍抽取 return x_filt[::M] # 示例使用 fs_audio_in = 44100 fs_audio_out = 48000 t_audio = np.arange(0, 1, 1/fs_audio_in) audio_in = np.sin(2*np.pi*1000*t_audio) + 0.5*np.sin(2*np.pi*5000*t_audio) audio_out = audio_resample(audio_in, fs_audio_in, fs_audio_out)

7.2 计算复杂度分析

不同实现方式的复杂度对比：

实际项目中，可以使用Python的timeit模块测量不同实现的运行时间：

import timeit setup = """ import numpy as np from scipy import signal x = np.random.randn(10000) b = signal.firwin(101, 0.1) """ direct_time = timeit.timeit('signal.lfilter(b, 1, x)[::3]', setup=setup, number=100) poly_time = timeit.timeit('signal.lfilter(b[::3], 1, x[::3])', setup=setup, number=100) print(f"直接方法: {direct_time:.4f}s, 多相方法: {poly_time:.4f}s")

8. 高级话题：抗混叠滤波器的优化设计

8.1 滤波器类型选择

不同滤波器类型的特性对比：

8.2 半带滤波器优化

对于2倍抽取/内插，半带滤波器可以显著减少计算量：

def halfband_filter_design(): """设计高效的半带滤波器""" # 半带滤波器的特点：约一半系数为零 taps = signal.remez(32, [0, 0.4, 0.6, 1], [1, 0], weight=[1, 100]) # 验证半带特性 print("非零系数数量:", np.sum(np.abs(taps) > 1e-6)) return taps hb_filter = halfband_filter_design()

8.3 自适应滤波技术

在非平稳信号处理中，可以考虑使用自适应滤波器动态调整截止频率：

def adaptive_decimation(x, fs, target_fs): """自适应抽取示例""" D = fs // target_fs if D < 1: return x, fs # 估计信号最高频率成分 spectrum = np.abs(np.fft.fft(x)) max_freq = np.argmax(spectrum[:len(spectrum)//2]) * fs / len(x) # 动态设置滤波器截止频率 cutoff = min(target_fs/2, max_freq*1.2) nyq = fs/2 b = signal.remez(101, [0, 0.9*cutoff/nyq, 1.1*cutoff/nyq, 1], [1, 0]) # 滤波和抽取 x_filt = signal.lfilter(b, 1, x) return x_filt[::D], target_fs

9. 可视化工具开发

为了更直观地理解采样率转换过程，我们可以开发一个交互式可视化工具：

from ipywidgets import interact, FloatSlider def interactive_resampling(f1=300, f2=1500, D=3, I=2): """交互式采样率转换演示""" fs = 6000 t = np.arange(0, 0.01, 1/fs) x = np.sin(2*np.pi*f1*t) + 0.5*np.sin(2*np.pi*f2*t) # 抽取 cutoff_dec = fs/(2*D) b_dec, a_dec = signal.butter(6, cutoff_dec/fs*2) x_dec = signal.lfilter(b_dec, a_dec, x)[::D] # 内插 x_zeros = np.zeros(len(x)*I) x_zeros[::I] = x cutoff_int = fs/2 b_int, a_int = signal.butter(6, cutoff_int/(fs*I)*2) x_int = signal.lfilter(b_int, a_int, x_zeros) # 绘图 plt.figure(figsize=(15, 10)) # 原始信号 plt.subplot(3, 2, 1) plt.plot(t, x) plt.title(f"原始信号 (fs={fs}Hz)") plt.subplot(3, 2, 2) plot_spectrum(x, fs, "原始频谱") # 抽取后 plt.subplot(3, 2, 3) plt.plot(np.arange(len(x_dec))/(fs/D), x_dec) plt.title(f"{D}倍抽取后 (fs={fs/D}Hz)") plt.subplot(3, 2, 4) plot_spectrum(x_dec, fs/D, "抽取后频谱") # 内插后 plt.subplot(3, 2, 5) plt.plot(np.arange(len(x_int))/(fs*I), x_int) plt.title(f"{I}倍内插后 (fs={fs*I}Hz)") plt.subplot(3, 2, 6) plot_spectrum(x_int, fs*I, "内插后频谱") plt.tight_layout() plt.show() # 创建交互界面 interact(interactive_resampling, f1=FloatSlider(min=100, max=2000, step=100, value=300), f2=FloatSlider(min=100, max=2000, step=100, value=1500), D=FloatSlider(min=2, max=5, step=1, value=3), I=FloatSlider(min=2, max=4, step=1, value=2))

10. 实际工程中的挑战与解决方案

10.1 有限字长效应

在FPGA或DSP硬件实现时，需要考虑有限精度带来的影响：

1. 系数量化误差：滤波器系数舍入误差
2. 算术溢出：定点运算中的溢出风险
3. 舍入噪声：累积计算误差

解决方案：

• 使用系数缩放和饱和算术
• 增加内部字长
• 采用噪声整形技术

10.2 实时性要求

对于实时系统，需要优化处理延迟：

1. 分段处理：将长信号分块处理
2. 流水线设计：重叠滤波与抽取/内插操作
3. 并行计算：利用多核处理器或GPU加速

10.3 非理想信号处理

真实世界信号的挑战：

• 非平稳信号（时变频谱）
• 带外噪声
• 时钟抖动

应对策略：

• 自适应滤波
• 过采样技术
• 抖动抑制算法

11. 性能优化技巧

11.1 内存高效实现

避免不必要的数组拷贝：

def memory_efficient_decimate(x, D, b): """内存高效的抽取实现""" # 原位操作减少内存使用 x = np.ascontiguousarray(x) x_filtered = signal.lfilter(b, 1, x, overwrite_x=True) return x_filtered[::D]

11.2 多线程处理

利用多核CPU并行处理：

from multiprocessing import Pool def parallel_resample(chunk): """并行处理的块处理函数""" x_chunk, b, D = chunk return signal.lfilter(b, 1, x_chunk)[::D] def parallel_decimation(x, D, b, n_workers=4): """并行抽取实现""" chunk_size = len(x) // n_workers chunks = [(x[i*chunk_size:(i+1)*chunk_size], b, D) for i in range(n_workers)] with Pool(n_workers) as p: results = p.map(parallel_resample, chunks) return np.concatenate(results)

11.3 GPU加速

使用CuPy库实现GPU加速：

try: import cupy as cp def gpu_resample(x, D, I, b): """GPU加速的重采样""" x_gpu = cp.asarray(x) b_gpu = cp.asarray(b) # 抽取 if D > 1: x_filtered = cp.zeros_like(x_gpu) for i in range(len(b)): x_filtered[i:] += b_gpu[i] * x_gpu[:len(x_gpu)-i] x_dec = x_filtered[::D] else: x_dec = x_gpu # 内插 if I > 1: x_up = cp.zeros(len(x_dec)*I) x_up[::I] = x_dec x_filtered = cp.zeros_like(x_up) for i in range(len(b)): x_filtered[i:] += b_gpu[i] * x_up[:len(x_up)-i] return cp.asnumpy(x_filtered) else: return cp.asnumpy(x_dec) except ImportError: print("CuPy not available, using CPU version")

12. 测试与验证框架

12.1 单元测试设计

确保各处理阶段的正确性：

import unittest class TestResampling(unittest.TestCase): def test_decimation(self): fs = 1000 t = np.arange(0, 1, 1/fs) x = np.sin(2*np.pi*100*t) D = 2 b = signal.firwin(31, 0.4) x_dec = signal.lfilter(b, 1, x)[::D] # 验证抽取后频率 spectrum = np.abs(np.fft.fft(x_dec)) max_freq = np.argmax(spectrum[:len(spectrum)//2]) * (fs/D) / len(x_dec) self.assertAlmostEqual(max_freq, 100, delta=5) def test_interpolation(self): fs = 1000 t = np.arange(0, 0.1, 1/fs) x = np.sin(2*np.pi*200*t) I = 3 x_up = np.zeros(len(x)*I) x_up[::I] = x b = signal.firwin(31, 1/I) x_int = signal.lfilter(b, 1, x_up) # 验证内插后频率 spectrum = np.abs(np.fft.fft(x_int)) max_freq = np.argmax(spectrum[:len(spectrum)//2]) * (fs*I) / len(x_int) self.assertAlmostEqual(max_freq, 200, delta=5) if __name__ == '__main__': unittest.main()

12.2 性能基准测试

比较不同实现的处理速度：

def benchmark(): """性能基准测试""" fs = 48000 x = np.random.randn(10*fs) # 10秒音频 # 测试不同抽取方法 print("=== 抽取性能测试 ===") b = signal.firwin(101, 0.4) # 直接方法 t1 = %timeit -o -r 3 signal.lfilter(b, 1, x)[::2] # 多相方法 b_poly = b.reshape(-1, 2).sum(axis=1) # 多相分解 t2 = %timeit -o -r 3 np.convolve(x[::2], b_poly, mode='same') print(f"直接方法: {t1.average:.4f}s, 多相方法: {t2.average:.4f}s") print(f"加速比: {t1.average/t2.average:.2f}x")

13. 完整Jupyter Notebook实现

将所有功能整合到一个可执行的Jupyter Notebook中：

# %% [markdown] # # 完整的信号抽取与内插实验 # # 这个笔记本实现了信号的3倍抽取和2倍内插操作，包含完整的时域和频域分析。 # %% import numpy as np import matplotlib.pyplot as plt from scipy import signal from IPython.display import display, Markdown # %% [markdown] # ## 1. 信号生成 # %% def generate_test_signal(fs, duration, freqs): """生成多频测试信号""" t = np.arange(0, duration, 1/fs) components = [np.sin(2*np.pi*f*t) for f in freqs] return t, np.sum(components, axis=0) fs = 6000 duration = 0.1 freqs = [300, 800, 1500] t, x = generate_test_signal(fs, duration, freqs) # %% [markdown] # ## 2. 抽取与内插函数定义 # %% def decimate_signal(x, D, fs): """带抗混叠滤波的抽取""" nyq = fs / 2 cutoff = fs / (2 * D) b = signal.remez(101, [0, 0.9*cutoff, 1.1*cutoff, nyq], [1, 0], fs=fs) x_filt = signal.filtfilt(b, 1, x) return x_filt[::D], fs // D def interpolate_signal(x, I, fs): """带镜像滤波的内插""" x_up = np.zeros(len(x) * I) x_up[::I] = x cutoff = fs / 2 nyq = fs * I / 2 b = signal.remez(101, [0, 0.9*cutoff, 1.1*cutoff, nyq], [1, 0], fs=nyq*2) return signal.filtfilt(b, 1, x_up), fs * I # %% [markdown] # ## 3. 执行处理并可视化 # %% # 3倍抽取 x_dec, fs_dec = decimate_signal(x, 3, fs) # 2倍内插 x_int, fs_int = interpolate_signal(x, 2, fs) # 组合操作：先3倍抽取再2倍内插 x_dec, _ = decimate_signal(x, 3, fs) x_combined, fs_combined = interpolate_signal(x_dec, 2, fs//3) # %% [markdown] # ## 4. 结果分析 # %% def plot_comparison(): fig, axs = plt.subplots(3, 2, figsize=(15, 12)) # 原始信号 axs[0,0].plot(t[:200], x[:200]) axs[0,0].set_title("Original Signal") axs[0,1].magnitude_spectrum(x, Fs=fs, scale='dB') axs[0,1].set_xlim(0, fs/2) # 抽取后 t_dec = t[::3] axs[1,0].plot(t_dec[:70], x_dec[:70], 'ro-') axs[1,0].set_title("Decimated by 3") axs[1,1].magnitude_spectrum(x_dec, Fs=fs_dec, scale='dB') axs[1,1].set_xlim(0, fs_dec/2) # 内插后 t_int = np.arange(0, duration, 1/fs_int) axs[2,0].plot(t_int[:400], x_int[:400], 'g.-') axs[2,0].set_title("Interpolated by 2") axs[2,1].magnitude_spectrum(x_int, Fs=fs_int, scale='dB') axs[2,1].set_xlim(0, fs_int/2) plt.tight_layout() plt.show() plot_comparison() # %% [markdown] # ## 5. SNR计算 # %% def calculate_snr(original, processed, fs_orig, fs_proc): """计算信号处理前后的信噪比""" # 对齐信号 min_len = min(len(original), int(len(processed) * fs_orig / fs_proc)) original = original[:min_len] processed = processed[:int(min_len * fs_proc / fs_orig)] # 计算功率 signal_power = np.sum(original**2) noise_power = np.sum((original - processed[:len(original)])**2) return 10 * np.log10(signal_power / noise_power) snr_dec = calculate_snr(x, x_dec, fs, fs_dec) snr_int = calculate_snr(x, x_int, fs, fs_int) snr_combined = calculate_snr(x, x_combined, fs, fs_combined) display(Markdown(f""" ### 信噪比结果: - 3倍抽取: {snr_dec:.2f} dB - 2倍内插: {snr_int:.2f} dB - 3倍抽取+2倍内插: {snr_combined:.2f} dB """))