模拟PPG生物信号，利用短时傅里叶变换滤波，提取心率

1、利用短时傅里叶变换滤波提取心率def test06(): import numpy as np import matplotlib.pyplot as plt from scipy import signal from scipy.signal import find_peaks # ---------- 参数设置 ---------- fs 1000 # 采样率 (Hz) duration 10 # 信号时长 (秒) t np.linspace(0, duration, int(fs * duration), endpointFalse) # ---------- 1. 生成模拟PPG信号同前---------- f_ppg 1.0 ppg (0.5 * np.sin(2 * np.pi * f_ppg * t) 0.2 * np.sin(2 * np.pi * 2 * f_ppg * t) 0.1 * np.sin(2 * np.pi * 3 * f_ppg * t)) baseline 0.3 * np.sin(2 * np.pi * 0.1 * t) noise 0.2 * np.random.randn(len(t)) power_line 0.15 * np.sin(2 * np.pi * 50 * t) signal_raw ppg baseline noise power_line # ---------- 2. STFT滤波同前---------- window hann nperseg 256 noverlap nperseg // 2 nfft nperseg f, t_stft, Zxx signal.stft(signal_raw, fsfs, windowwindow, npersegnperseg, noverlapnoverlap, nfftnfft, boundaryNone, paddedFalse) # 设计掩码保留0.5-5 Hz f_low, f_high 0.5, 5.0 mask np.ones_like(Zxx, dtypefloat) freq_indices np.where((f f_low) | (f f_high))[0] mask[freq_indices, :] 0.0 Zxx_filtered Zxx * mask # ISTFT重建 t_reconstructed, signal_filtered signal.istft(Zxx_filtered, fsfs, windowwindow, npersegnperseg, noverlapnoverlap, nfftnfft, boundaryFalse) # 对齐长度 min_len min(len(signal_raw), len(signal_filtered)) signal_raw signal_raw[:min_len] signal_filtered signal_filtered[:min_len] t t[:min_len] # ---------- 3. 峰谷检测 ---------- # 检测波峰收缩期峰值 # 设置最小间隔对应最大心率200 bpm 间隔至少 60/200 0.3秒 300个采样点 min_distance_peaks int(0.3 * fs) # 300 samples # 可选高度阈值这里设为0.2因为滤波后信号幅值约为0.5 peaks, _ find_peaks(signal_filtered, distancemin_distance_peaks, height0.2) # 检测波谷舒张期最低点反转信号 signal_inverted -signal_filtered valleys, _ find_peaks(signal_inverted, distancemin_distance_peaks, height0.2) # 提取峰谷的时间点 peak_times t[peaks] valley_times t[valleys] # 计算心率相邻峰值间隔转换为bpm if len(peaks) 2: peak_intervals np.diff(peak_times) # 间隔秒 instant_hr 60.0 / peak_intervals # 瞬时心率 (bpm) mean_hr np.mean(instant_hr) else: instant_hr np.array([]) mean_hr np.nan print(f检测到波峰数量{len(peaks)}) print(f平均心率{mean_hr:.1f} bpm) # ---------- 4. 结果可视化 ---------- plt.figure(figsize(16, 12)) # 子图1原始信号 plt.subplot(3, 3, 1) plt.plot(t, signal_raw) plt.title(原始模拟PPG信号) plt.xlabel(时间 (s)) plt.ylabel(幅值) # 子图2滤波后信号 plt.subplot(3, 3, 2) plt.plot(t, signal_filtered) plt.title(STFT滤波后信号 (0.5-5 Hz)) plt.xlabel(时间 (s)) plt.ylabel(幅值) # 子图3原始时频谱 plt.subplot(3, 3, 3) plt.pcolormesh(t_stft, f, np.abs(Zxx), shadinggouraud) plt.title(原始信号 STFT 幅度谱) plt.xlabel(时间 (s)) plt.ylabel(频率 (Hz)) plt.ylim(0, 100) plt.colorbar(label幅度) # 子图4滤波后时频谱 plt.subplot(3, 3, 4) plt.pcolormesh(t_stft, f, np.abs(Zxx_filtered), shadinggouraud) plt.title(滤波后 STFT 幅度谱) plt.xlabel(时间 (s)) plt.ylabel(频率 (Hz)) plt.ylim(0, 100) plt.colorbar(label幅度) # 子图5滤波后信号 峰谷标记 plt.subplot(3, 3, 5) plt.plot(t, signal_filtered, label滤波后信号) plt.plot(peak_times, signal_filtered[peaks], ro, label波峰) plt.plot(valley_times, signal_filtered[valleys], go, label波谷) plt.title(滤波后信号及峰谷检测) plt.xlabel(时间 (s)) plt.ylabel(幅值) plt.legend() # 子图6前2秒放大显示峰谷 plt.subplot(3, 3, 6) plt.plot(t[:2000], signal_filtered[:2000], label滤波后信号) plt.plot(peak_times[peak_times 2], signal_filtered[peaks][peak_times 2], ro, label波峰) plt.plot(valley_times[valley_times 2], signal_filtered[valleys][valley_times 2], go, label波谷) plt.title(前2秒峰谷细节) plt.xlabel(时间 (s)) plt.ylabel(幅值) plt.legend() # 子图7瞬时心率变化如果至少有两个峰 plt.subplot(3, 3, 7) if len(instant_hr) 0: plt.plot(peak_times[1:], instant_hr, b-o, markersize4) plt.axhline(ymean_hr, colorr, linestyle--, labelf平均心率: {mean_hr:.1f} bpm) plt.title(瞬时心率) plt.xlabel(时间 (s)) plt.ylabel(心率 (bpm)) plt.legend() else: plt.text(0.5, 0.5, 无法计算心率, hacenter, vacenter, transformplt.gca().transAxes) plt.title(瞬时心率) # 子图8残差原始-滤波后 plt.subplot(3, 3, 8) plt.plot(t, signal_raw - signal_filtered) plt.title(残差 (原始 - 滤波后)) plt.xlabel(时间 (s)) plt.ylabel(幅值) # 子图9峰谷间隔分布直方图可选 plt.subplot(3, 3, 9) if len(peak_intervals) 0: plt.hist(peak_intervals, bins10, edgecolorblack) plt.title(峰峰间隔分布) plt.xlabel(间隔 (s)) plt.ylabel(频数) else: plt.text(0.5, 0.5, 无间隔数据, hacenter, vacenter, transformplt.gca().transAxes) plt.title(峰峰间隔分布) plt.tight_layout() plt.show()2、信号的低频和高频的提取import numpy as np # 原始信号 x np.array([1, 2, 3, 7, 4, 3, 9, 8, 6], dtypefloat) # 1. 计算 FFT X np.fft.fft(x) N len(x) # 2. 构造低频频谱掩码 Y np.zeros(N, dtypecomplex) # 初始化全零 # 保留直流分量 k0 Y[0] X[0] # 保留 k1 和 k2 及其对应的对称分量对于实信号N-k 位置 Y[1] X[1] Y[2] X[2] Y[N-1] X[N-1] # k8 对应 N-1 Y[N-2] X[N-2] # k7 对应 N-2 # 其余 k3,4,5,6 保持为零 # 3. 对低频频谱进行逆变换得到低频信号 x_low np.fft.ifft(Y).real # 取实部理论上应为实数 # 4. 高频信号 原始信号 - 低频信号 x_high x - x_low # 输出结果保留三位小数 print(原始信号 x , x) print(低频信号 x_low , np.round(x_low, 3)) print(高频信号 x_high , np.round(x_high, 3)) # 验证低频高频是否等于原信号 print(\n低频高频 , np.round(x_low x_high, 3)) print(低频信号之和 , np.sum(x_low)) print(高频信号之和 , np.sum(x_high))实序列信号是共轭对称只需要计算其中的一半即可假设N8k: 0 1 2 3 4 5 6 7 实 复 复 复 实 复* 复* 复*假设N9k: 0 1 2 3 4 5 6 7 8 实 复 复 复 复 复* 复* 复* 复*3、主要区别总结项目偶数 NN奇数 NN独立频率点数单边谱N/21(N1)/2奈奎斯特频率点存在对应 kN/2为实数不存在频率成对情况除 k0和 kN/2外其余成对除 k0外其余全部成对频谱绘制通常画到 fs/2 为止包含奈奎斯特点画到略低于 fs/2的最大频率点物理意义可准确表示采样频率一半的成分无法准确表示正好 fs/2的频率a、总结表格序列类型时域对称性频域对称性频域数值特点实偶序列x[n]x[N−n]X[k]为偶对称实部偶虚部为零实奇序列x[n]−x[N−n]X[k] 为奇对称实部为零虚部奇一般实序列可分解为偶部奇部实部偶虚部奇同时存在实部和虚部b、真实的应用场景偶对称实序列只需计算余弦变换如 DCT可省去一半计算量。奇对称实序列只需计算正弦变换如 DST。一般实序列利用共轭对称性只需计算 k0k0 到 N/2N/2 的频率点即“单边谱”。理解这些区别能帮助你在信号处理中快速判断频谱结构并选用合适的变换如 DCT/DST来压缩或分析信号。共轭对称性就是幅度相等相位相反。a3j --------a-3j(实部相等虚部相反)