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

📅 2026/7/11 7:50:05
Python实现信号3倍抽取与2倍内插:从时域波形到频谱验证
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, axis0) # 合成信号 # 可视化 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降低信号采样率的过程包含两个关键步骤抗混叠滤波防止抽取后出现频谱混叠采样率降低每隔D-1个样本保留一个样本数学表达式为 $$ y[n] x[nD] $$其中D为抽取因子。在频域抽取会导致频谱扩展和周期性重复$$ Y(e^{j\omega}) \frac{1}{D}\sum_{k0}^{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, btypelow) # 应用滤波器 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, markerfmtro, 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 关键参数对比参数原始信号抽取后信号采样率6000 Hz2000 Hz奈奎斯特频率3000 Hz1000 Hz数据点数量600200最高可处理频率3000 Hz1000 Hz3. 信号内插从理论到实现3.1 内插的数学原理内插Interpolation是按整数因子I增加采样率的过程包含三个主要步骤零值插入在原始样本间插入I-1个零镜像滤波去除由零插入引入的高频镜像采样率提升最终采样率为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, btypelow) # 应用滤波器 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, markerfmtro, 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 内插过程中的频谱变化零插入阶段在原始频谱基础上会在fs_original/2到fs_interpolated/2之间产生镜像频谱滤波阶段低通滤波器保留原始频带(0到fs_original/2)消除镜像成分最终频谱与原始信号频谱形状相同但横轴范围扩展到fs_interpolated/24. 综合应用抽取与内插联合操作在实际系统中经常需要组合使用抽取和内插来实现非整数倍的采样率转换。让我们实现一个先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 滤波器设计考量在实际应用中滤波器设计需要权衡三个关键参数过渡带宽越窄则滤波器阶数越高计算量越大阻带衰减要求越高则滤波器阶数越高通带波纹越小则信号失真越小但实现难度增加推荐使用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, weightweight)5.2 多级处理策略对于大的采样率转换因子如D10采用多级处理可以显著降低计算复杂度多级抽取将D分解为D1×D2×...×Dn每级使用较小的抽取因子多级内插同理分解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 实时处理优化在实时系统中可以采用以下优化策略多相滤波减少计算量特别适合内插操作环形缓冲区高效处理连续数据流定点运算在嵌入式系统中提高效率多相滤波实现示例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 y6. 验证与调试技巧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 常见问题排查频谱混叠检查抗混叠滤波器的截止频率是否设置正确幅度失真确认滤波器通带波纹是否过大相位失真考虑使用零相位滤波(filtfilt代替lfilter)计算延迟多级处理时注意各阶段的延迟累积零相位滤波示例# 使用零相位滤波减少相位失真 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], fsnyq*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 计算复杂度分析不同实现方式的复杂度对比方法每输出样本乘法次数适用场景直接滤波抽取N (滤波器阶数)简单系统多相抽取N/D大D值抽取多相内插N/I大I值内插FFT卷积log2(L) (L为FFT长度)长滤波器实际项目中可以使用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], setupsetup, number100) poly_time timeit.timeit(signal.lfilter(b[::3], 1, x[::3]), setupsetup, number100) print(f直接方法: {direct_time:.4f}s, 多相方法: {poly_time:.4f}s)8. 高级话题抗混叠滤波器的优化设计8.1 滤波器类型选择不同滤波器类型的特性对比滤波器类型优点缺点适用场景Butterworth最大平坦通带过渡带较宽要求平坦通带Chebyshev I陡峭过渡带通带波纹需要锐截止Chebyshev II平坦通带阻带波纹需要良好阻带Elliptic最陡过渡带通阻带波纹严格资源限制FIR线性相位高阶数实现锐截止零相位需求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_fs9. 可视化工具开发为了更直观地理解采样率转换过程我们可以开发一个交互式可视化工具from ipywidgets import interact, FloatSlider def interactive_resampling(f1300, f21500, D3, I2): 交互式采样率转换演示 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, f1FloatSlider(min100, max2000, step100, value300), f2FloatSlider(min100, max2000, step100, value1500), DFloatSlider(min2, max5, step1, value3), IFloatSlider(min2, max4, step1, value2))10. 实际工程中的挑战与解决方案10.1 有限字长效应在FPGA或DSP硬件实现时需要考虑有限精度带来的影响系数量化误差滤波器系数舍入误差算术溢出定点运算中的溢出风险舍入噪声累积计算误差解决方案使用系数缩放和饱和算术增加内部字长采用噪声整形技术10.2 实时性要求对于实时系统需要优化处理延迟分段处理将长信号分块处理流水线设计重叠滤波与抽取/内插操作并行计算利用多核处理器或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_xTrue) 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_workers4): 并行抽取实现 chunk_size len(x) // n_workers chunks [(x[i*chunk_size:(i1)*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, delta5) 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, delta5) 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(axis1) # 多相分解 t2 %timeit -o -r 3 np.convolve(x[::2], b_poly, modesame) 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, axis0) 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], fsfs) 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], fsnyq*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, Fsfs, scaledB) 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, Fsfs_dec, scaledB) 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, Fsfs_int, scaledB) 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 ))