用Python和Matplotlib动手画窄带与宽带干扰的频谱/时频图附完整代码在信号处理领域理解干扰信号的特性至关重要。无论是通信系统设计还是雷达信号分析工程师们都需要直观掌握窄带干扰(NBI)和宽带干扰(WBI)的时频特性。本文将带你用Python从零开始构建这些干扰模型并通过Matplotlib可视化它们的频谱和时频特征。1. 环境准备与基础设置在开始之前确保你的Python环境已安装以下库import numpy as np import matplotlib.pyplot as plt from scipy import signal from scipy.fft import fft, fftshift我们将使用这些工具来生成信号并进行分析。首先设置一些基础参数# 通用参数设置 sample_rate 1000 # 采样率(Hz) duration 1.0 # 信号持续时间(s) n_samples int(sample_rate * duration) # 总采样点数 t np.linspace(0, duration, n_samples, endpointFalse) # 时间轴注意采样率应至少是信号最高频率的两倍奈奎斯特准则这里设为1000Hz足够我们演示大多数情况。2. 窄带干扰(NBI)的建模与可视化窄带干扰通常表现为频谱中的尖峰我们可以用单频信号或多个单频信号的叠加来模拟。2.1 单频干扰生成def generate_nbi(freq, amplitude1.0, phase0.0): 生成单频窄带干扰信号 return amplitude * np.exp(1j * (2 * np.pi * freq * t phase)) # 示例生成1个250Hz的单频干扰 nbi_single generate_nbi(freq250)2.2 多频窄带干扰现实中的NBI往往是多个单频信号的组合# 生成包含3个频率分量的NBI frequencies [100, 250, 400] # 各分量频率(Hz) amplitudes [0.8, 1.0, 0.6] # 各分量幅度 phases [0, np.pi/4, np.pi/2] # 各分量相位 nbi_multi sum(generate_nbi(f, a, p) for f, a, p in zip(frequencies, amplitudes, phases))2.3 频谱分析使用FFT分析NBI的频谱特性def plot_spectrum(signal, title): n len(signal) freq np.fft.fftfreq(n, d1/sample_rate) spectrum np.abs(fftshift(fft(signal))) plt.figure(figsize(10, 4)) plt.plot(freq, spectrum) plt.title(title) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.grid(True) plt.xlim(-500, 500) plt.show() plot_spectrum(nbi_multi, Multi-component NBI Spectrum)关键观察点频谱在100Hz、250Hz和400Hz处出现明显尖峰各峰高度对应信号分量幅度其他频率成分几乎为零体现窄带特性2.4 时频分析(STFT)短时傅里叶变换能展示频率随时间的变化def plot_stft(signal, title): f, t, Zxx signal.stft(signal, fssample_rate, nperseg64) plt.figure(figsize(10, 4)) plt.pcolormesh(t, f, np.abs(Zxx), shadinggouraud) plt.title(title) plt.ylabel(Frequency [Hz]) plt.xlabel(Time [sec]) plt.colorbar(labelMagnitude) plt.ylim(0, 500) plt.show() plot_stft(nbi_multi, NBI STFT (Time-Frequency Representation))NBI的时频图显示三条平行于时间轴的直线印证了其频率不随时间变化的特性。3. 宽带干扰(WBI)的建模与可视化宽带干扰占据较宽的频带我们重点分析两种常见类型线性调频和正弦调频干扰。3.1 线性调频宽带干扰线性调频信号的频率随时间线性变化def generate_lfm(start_freq, end_freq, amplitude1.0): 生成线性调频宽带干扰 rate (end_freq - start_freq) / duration # 调频斜率 phase 2 * np.pi * (start_freq * t 0.5 * rate * t**2) return amplitude * np.exp(1j * phase) # 示例生成从50Hz到450Hz的线性调频信号 wbi_lfm generate_lfm(50, 450)3.1.1 频谱特性分析plot_spectrum(wbi_lfm, LFM WBI Spectrum)LFM信号的频谱呈现宽带的平台特征能量分布在50-450Hz范围内。3.1.2 时频特性分析plot_stft(wbi_lfm, LFM WBI STFT)时频图上清晰显示频率随时间线性增长的特征斜率为(450-50)/1400 Hz/s。3.2 正弦调频宽带干扰正弦调频信号的频率按正弦规律变化def generate_sfm(carrier_freq, mod_freq, mod_index, amplitude1.0): 生成正弦调频宽带干扰 phase 2 * np.pi * carrier_freq * t mod_index * np.sin(2 * np.pi * mod_freq * t) return amplitude * np.exp(1j * phase) # 示例载波300Hz调制频率5Hz调制指数10 wbi_sfm generate_sfm(300, 5, 10)3.2.1 频谱特性分析plot_spectrum(wbi_sfm, SFM WBI Spectrum)SFM信号的频谱呈现对称的宽带特征中心在载波频率300Hz处带宽与调制指数相关。3.2.2 时频特性分析plot_stft(wbi_sfm, SFM WBI STFT)时频图显示频率按正弦规律波动直观反映了调制过程。4. 干扰信号参数影响分析不同参数会显著影响干扰信号的时频特性理解这些影响对干扰识别和抑制至关重要。4.1 调制系数对SFM信号的影响我们固定载波频率(300Hz)和调制频率(5Hz)改变调制指数调制指数频谱宽度时频图波动幅度2窄小5中等中等10宽大# 对比不同调制指数的SFM信号 for beta in [2, 5, 10]: sfm_signal generate_sfm(300, 5, beta) plot_stft(sfm_signal, fSFM with β{beta})4.2 调频率对LFM信号的影响保持起止频率不变(50-450Hz)改变信号持续时间持续时间(s)调频率(Hz/s)频谱形状0.5800更陡峭1.0400中等2.0200更平缓# 不同持续时间的LFM信号比较 for dur in [0.5, 1.0, 2.0]: t_local np.linspace(0, dur, int(sample_rate * dur), endpointFalse) rate 400 / dur # 保持总频率跨度相同 phase 2 * np.pi * (50 * t_local 0.5 * rate * t_local**2) lfm_signal np.exp(1j * phase) plot_stft(lfm_signal, fLFM with duration{dur}s)5. 实际应用中的注意事项在真实场景中分析干扰信号时有几个关键点需要考虑频谱泄露问题使用合适的窗函数(如Hamming窗)确保采样时长包含完整周期# 加窗减少频谱泄露 window np.hamming(len(nbi_multi)) plot_spectrum(nbi_multi * window, Windowed NBI Spectrum)时频分辨率权衡STFT中时间分辨率与频率分辨率相互制约需要根据具体需求选择适当的窗口长度# 不同窗口长度的STFT比较 for nperseg in [32, 64, 128]: f, t, Zxx signal.stft(nbi_multi, fssample_rate, npersegnperseg) plt.pcolormesh(t, f, np.abs(Zxx), shadinggouraud) plt.title(fSTFT with {nperseg} points window) plt.show()计算效率优化对于长信号考虑使用重叠分段处理合理选择FFT点数(通常是2的幂次)# 优化FFT计算 def efficient_spectrum(signal, n_fft1024): spectrum np.abs(fftshift(fft(signal, nn_ff4))) freq np.fft.fftshift(np.fft.fftfreq(n_fft, d1/sample_rate)) return freq, spectrum6. 完整代码示例以下是生成和可视化NBI与WBI的完整Python脚本import numpy as np import matplotlib.pyplot as plt from scipy import signal from scipy.fft import fft, fftshift # 参设置 sample_rate 1000 # Hz duration 1.0 # seconds t np.linspace(0, duration, int(sample_rate * duration), endpointFalse) def generate_signals(): 生成各种干扰信号 # NBI - 多频 nbi_freqs [100, 250, 400] nbi_amps [0.8, 1.0, 0.6] nbi_phases [0, np.pi/4, np.pi/2] nbi sum(a * np.exp(1j*(2*np.pi*f*t p)) for f, a, p in zip(nbi_freqs, nbi_amps, nbi_phases)) # WBI - 线性调频 lfm np.exp(1j * 2*np.pi * (50*t 0.5*400*t**2)) # WBI - 正弦调频 sfm np.exp(1j * (2*np.pi*300*t 10*np.sin(2*np.pi*5*t))) return {NBI: nbi, LFM: lfm, SFM: sfm} def analyze_signals(signals): 分析并绘制信号特性 for name, sig in signals.items(): # 时域波形(实部) plt.figure(figsize(10, 7)) plt.subplot(3, 1, 1) plt.plot(t, np.real(sig)) plt.title(f{name} Signal - Time Domain) plt.xlabel(Time (s)) plt.ylabel(Amplitude) # 频谱 plt.subplot(3, 1, 2) freq np.fft.fftshift(np.fft.fftfreq(len(t), d1/sample_rate)) spectrum np.abs(fftshift(fft(sig))) plt.plot(freq, spectrum) plt.title(Frequency Spectrum) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.xlim(-500, 500) # 时频图 plt.subplot(3, 1, 3) f, t_stft, Zxx signal.stft(sig, fssample_rate, nperseg64) plt.pcolormesh(t_stft, f, np.abs(Zxx), shadinggouraud) plt.title(STFT (Time-Frequency Analysis)) plt.xlabel(Time (s)) plt.ylabel(Frequency (Hz)) plt.ylim(0, 500) plt.colorbar() plt.tight_layout() plt.show() if __name__ __main__: signals generate_signals() analyze_signals(signals)运行这段代码将生成三种干扰信号的时域波形、频谱和时频图方便对比分析它们的特性差异。
用Python和Matplotlib动手画:窄带与宽带干扰的频谱/时频图(附完整代码)
用Python和Matplotlib动手画窄带与宽带干扰的频谱/时频图附完整代码在信号处理领域理解干扰信号的特性至关重要。无论是通信系统设计还是雷达信号分析工程师们都需要直观掌握窄带干扰(NBI)和宽带干扰(WBI)的时频特性。本文将带你用Python从零开始构建这些干扰模型并通过Matplotlib可视化它们的频谱和时频特征。1. 环境准备与基础设置在开始之前确保你的Python环境已安装以下库import numpy as np import matplotlib.pyplot as plt from scipy import signal from scipy.fft import fft, fftshift我们将使用这些工具来生成信号并进行分析。首先设置一些基础参数# 通用参数设置 sample_rate 1000 # 采样率(Hz) duration 1.0 # 信号持续时间(s) n_samples int(sample_rate * duration) # 总采样点数 t np.linspace(0, duration, n_samples, endpointFalse) # 时间轴注意采样率应至少是信号最高频率的两倍奈奎斯特准则这里设为1000Hz足够我们演示大多数情况。2. 窄带干扰(NBI)的建模与可视化窄带干扰通常表现为频谱中的尖峰我们可以用单频信号或多个单频信号的叠加来模拟。2.1 单频干扰生成def generate_nbi(freq, amplitude1.0, phase0.0): 生成单频窄带干扰信号 return amplitude * np.exp(1j * (2 * np.pi * freq * t phase)) # 示例生成1个250Hz的单频干扰 nbi_single generate_nbi(freq250)2.2 多频窄带干扰现实中的NBI往往是多个单频信号的组合# 生成包含3个频率分量的NBI frequencies [100, 250, 400] # 各分量频率(Hz) amplitudes [0.8, 1.0, 0.6] # 各分量幅度 phases [0, np.pi/4, np.pi/2] # 各分量相位 nbi_multi sum(generate_nbi(f, a, p) for f, a, p in zip(frequencies, amplitudes, phases))2.3 频谱分析使用FFT分析NBI的频谱特性def plot_spectrum(signal, title): n len(signal) freq np.fft.fftfreq(n, d1/sample_rate) spectrum np.abs(fftshift(fft(signal))) plt.figure(figsize(10, 4)) plt.plot(freq, spectrum) plt.title(title) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.grid(True) plt.xlim(-500, 500) plt.show() plot_spectrum(nbi_multi, Multi-component NBI Spectrum)关键观察点频谱在100Hz、250Hz和400Hz处出现明显尖峰各峰高度对应信号分量幅度其他频率成分几乎为零体现窄带特性2.4 时频分析(STFT)短时傅里叶变换能展示频率随时间的变化def plot_stft(signal, title): f, t, Zxx signal.stft(signal, fssample_rate, nperseg64) plt.figure(figsize(10, 4)) plt.pcolormesh(t, f, np.abs(Zxx), shadinggouraud) plt.title(title) plt.ylabel(Frequency [Hz]) plt.xlabel(Time [sec]) plt.colorbar(labelMagnitude) plt.ylim(0, 500) plt.show() plot_stft(nbi_multi, NBI STFT (Time-Frequency Representation))NBI的时频图显示三条平行于时间轴的直线印证了其频率不随时间变化的特性。3. 宽带干扰(WBI)的建模与可视化宽带干扰占据较宽的频带我们重点分析两种常见类型线性调频和正弦调频干扰。3.1 线性调频宽带干扰线性调频信号的频率随时间线性变化def generate_lfm(start_freq, end_freq, amplitude1.0): 生成线性调频宽带干扰 rate (end_freq - start_freq) / duration # 调频斜率 phase 2 * np.pi * (start_freq * t 0.5 * rate * t**2) return amplitude * np.exp(1j * phase) # 示例生成从50Hz到450Hz的线性调频信号 wbi_lfm generate_lfm(50, 450)3.1.1 频谱特性分析plot_spectrum(wbi_lfm, LFM WBI Spectrum)LFM信号的频谱呈现宽带的平台特征能量分布在50-450Hz范围内。3.1.2 时频特性分析plot_stft(wbi_lfm, LFM WBI STFT)时频图上清晰显示频率随时间线性增长的特征斜率为(450-50)/1400 Hz/s。3.2 正弦调频宽带干扰正弦调频信号的频率按正弦规律变化def generate_sfm(carrier_freq, mod_freq, mod_index, amplitude1.0): 生成正弦调频宽带干扰 phase 2 * np.pi * carrier_freq * t mod_index * np.sin(2 * np.pi * mod_freq * t) return amplitude * np.exp(1j * phase) # 示例载波300Hz调制频率5Hz调制指数10 wbi_sfm generate_sfm(300, 5, 10)3.2.1 频谱特性分析plot_spectrum(wbi_sfm, SFM WBI Spectrum)SFM信号的频谱呈现对称的宽带特征中心在载波频率300Hz处带宽与调制指数相关。3.2.2 时频特性分析plot_stft(wbi_sfm, SFM WBI STFT)时频图显示频率按正弦规律波动直观反映了调制过程。4. 干扰信号参数影响分析不同参数会显著影响干扰信号的时频特性理解这些影响对干扰识别和抑制至关重要。4.1 调制系数对SFM信号的影响我们固定载波频率(300Hz)和调制频率(5Hz)改变调制指数调制指数频谱宽度时频图波动幅度2窄小5中等中等10宽大# 对比不同调制指数的SFM信号 for beta in [2, 5, 10]: sfm_signal generate_sfm(300, 5, beta) plot_stft(sfm_signal, fSFM with β{beta})4.2 调频率对LFM信号的影响保持起止频率不变(50-450Hz)改变信号持续时间持续时间(s)调频率(Hz/s)频谱形状0.5800更陡峭1.0400中等2.0200更平缓# 不同持续时间的LFM信号比较 for dur in [0.5, 1.0, 2.0]: t_local np.linspace(0, dur, int(sample_rate * dur), endpointFalse) rate 400 / dur # 保持总频率跨度相同 phase 2 * np.pi * (50 * t_local 0.5 * rate * t_local**2) lfm_signal np.exp(1j * phase) plot_stft(lfm_signal, fLFM with duration{dur}s)5. 实际应用中的注意事项在真实场景中分析干扰信号时有几个关键点需要考虑频谱泄露问题使用合适的窗函数(如Hamming窗)确保采样时长包含完整周期# 加窗减少频谱泄露 window np.hamming(len(nbi_multi)) plot_spectrum(nbi_multi * window, Windowed NBI Spectrum)时频分辨率权衡STFT中时间分辨率与频率分辨率相互制约需要根据具体需求选择适当的窗口长度# 不同窗口长度的STFT比较 for nperseg in [32, 64, 128]: f, t, Zxx signal.stft(nbi_multi, fssample_rate, npersegnperseg) plt.pcolormesh(t, f, np.abs(Zxx), shadinggouraud) plt.title(fSTFT with {nperseg} points window) plt.show()计算效率优化对于长信号考虑使用重叠分段处理合理选择FFT点数(通常是2的幂次)# 优化FFT计算 def efficient_spectrum(signal, n_fft1024): spectrum np.abs(fftshift(fft(signal, nn_ff4))) freq np.fft.fftshift(np.fft.fftfreq(n_fft, d1/sample_rate)) return freq, spectrum6. 完整代码示例以下是生成和可视化NBI与WBI的完整Python脚本import numpy as np import matplotlib.pyplot as plt from scipy import signal from scipy.fft import fft, fftshift # 参设置 sample_rate 1000 # Hz duration 1.0 # seconds t np.linspace(0, duration, int(sample_rate * duration), endpointFalse) def generate_signals(): 生成各种干扰信号 # NBI - 多频 nbi_freqs [100, 250, 400] nbi_amps [0.8, 1.0, 0.6] nbi_phases [0, np.pi/4, np.pi/2] nbi sum(a * np.exp(1j*(2*np.pi*f*t p)) for f, a, p in zip(nbi_freqs, nbi_amps, nbi_phases)) # WBI - 线性调频 lfm np.exp(1j * 2*np.pi * (50*t 0.5*400*t**2)) # WBI - 正弦调频 sfm np.exp(1j * (2*np.pi*300*t 10*np.sin(2*np.pi*5*t))) return {NBI: nbi, LFM: lfm, SFM: sfm} def analyze_signals(signals): 分析并绘制信号特性 for name, sig in signals.items(): # 时域波形(实部) plt.figure(figsize(10, 7)) plt.subplot(3, 1, 1) plt.plot(t, np.real(sig)) plt.title(f{name} Signal - Time Domain) plt.xlabel(Time (s)) plt.ylabel(Amplitude) # 频谱 plt.subplot(3, 1, 2) freq np.fft.fftshift(np.fft.fftfreq(len(t), d1/sample_rate)) spectrum np.abs(fftshift(fft(sig))) plt.plot(freq, spectrum) plt.title(Frequency Spectrum) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.xlim(-500, 500) # 时频图 plt.subplot(3, 1, 3) f, t_stft, Zxx signal.stft(sig, fssample_rate, nperseg64) plt.pcolormesh(t_stft, f, np.abs(Zxx), shadinggouraud) plt.title(STFT (Time-Frequency Analysis)) plt.xlabel(Time (s)) plt.ylabel(Frequency (Hz)) plt.ylim(0, 500) plt.colorbar() plt.tight_layout() plt.show() if __name__ __main__: signals generate_signals() analyze_signals(signals)运行这段代码将生成三种干扰信号的时域波形、频谱和时频图方便对比分析它们的特性差异。