这是一个信号采样数学实验,你可以直观感受到冲击信号的时域和频域特征
原始信号是一个频率为180Hz附近的一个冲击性信号:
它的频谱,可能会超出你的想象,它的1x频率幅度可能并不最高的。频谱在高频展开
低通滤波后的波形几乎不变。但是此时的频谱:
注意1x, 2x, 3x谱线的峰值都更低,但是出现了直流分量以抵消掉高次谐波对1x的贡献:
相关的./gphelper/calc/gpFFT参见:GitCode - 全球开发者的开源社区,开源代码托管平台
#!/usr/bin/env python3 # -*- coding: utf-8 -*- # 获取当前脚本文件所在目录的父目录,并构建相对路径 import os import sys current_dir = os.path.dirname(os.path.abspath(__file__)) project_path = os.path.join(current_dir, '..') sys.path.append(project_path) sys.path.append(current_dir) sys.path.append('./gphelper/calc/') import numpy as np import matplotlib.pyplot as plt from matplotlib import rcParams from scipy.signal import get_window import gpFFT # 设置中文字体 rcParams['font.sans-serif'] = ['SimHei'] # 指定中文字体 rcParams['axes.unicode_minus'] = False # 解决负号问题 # 参数设置 fs = 20000 # 采样频率 (Hz) f_signal = 181.5 # 信号频率 (Hz) rms_signal = 5e-2 # 信号 RMS (m/s) duration = 1 # 信号持续时间 (秒) pulse_duration = 1 / (15 * f_signal) # 脉冲持续时间 (秒) # 生成时间序列 t = np.arange(0, duration, 1/fs) # 生成冲击性振动信号(脉冲信号) signal = np.zeros_like(t) firstEdge = 3 cntOfCycle = 0 while True: pulse_start = int(firstEdge+fs/f_signal*cntOfCycle) pulse_end = int(pulse_start+pulse_duration*fs) if(pulse_start>=fs): break if(pulse_end>fs): pulse_end = fs-1 signal[pulse_start:pulse_end] = rms_signal * np.sqrt(2) if(pulse_end>=fs): break; cntOfCycle += 1 # 生成高斯噪声 rms_noise = rms_signal * 0.01 noise = rms_noise * np.random.randn(len(t)) # 将信号和噪声相加 signal_with_noise = signal + noise # 绘图 plt.subplot(2, 2, 1) plt.plot(t, signal_with_noise, label='含噪声信号', color='red') plt.title('含噪声的冲击性振动信号') plt.xlabel('时间 (秒)') plt.ylabel('幅值') plt.legend() #频谱展示 (freq, fft_toshow_with_noise, fft_ac) = gpFFT.GetFFTOfSignal(t, signal_with_noise) plt.subplot(2, 2, 2) plt.plot(freq, fft_toshow_with_noise, label='冲击性振动信号_频谱') plt.title('冲击性振动信号') plt.xlabel('频率 (Hz)') plt.ylabel('幅值') plt.legend() (x, signal_filtered) = gpFFT.lowFilterSignal(t, signal_with_noise, 1*f_signal+3) plt.subplot(2, 2, 3) plt.plot(t, signal_with_noise, label='低通滤波', color='red') plt.title('含噪声的冲击性振动信号') plt.xlabel('时间 (秒)') plt.ylabel('幅值') plt.legend() (freq, fft_signal_filtered, fft_ac) = gpFFT.GetFFTOfSignal(t, signal_filtered) plt.subplot(2, 2, 4) plt.plot(freq, fft_signal_filtered, label='低通滤波频谱') plt.title('含噪声的冲击性振动信号') plt.xlabel('频率 (Hz)') plt.ylabel('幅值') plt.legend() plt.tight_layout() plt.show()