基于ECG信号的HRV时域与频域分析Matlab代码实现
一、核心流程
二、完整Matlab代码实现
1. 数据预处理（去噪与基线校正）
functionclean_ecg=preprocess_ecg(ecg,fs)% 带通滤波（0.5-40Hz）去除工频干扰和基线漂移[b,a]=butter(3,[0.540]/(fs/2));filtered=filtfilt(b,a,ecg);% 小波去噪（db4小波，5层分解）[c,l]=wavedec(filtered,5,'db4');denoised=waverec(wthresh(c,'s',3*std(c)/2),l,'db4');% 移动平均平滑（窗口宽度=50ms）window=round(0.05*fs);clean_ecg=movmean(denoised,window);end2. R波检测（改进Pan-Tompkins算法）
function[r_peaks,locs]=detect_r_peaks(ecg,fs)N=length(ecg);window=round(0.15*fs);% 积分窗口% 平方与积分squared=ecg.^2;integrated=filter(ones(1,window)/window,1,squared);% 动态阈值threshold=0.7*max(integrated(1:1000))+0.3*std(integrated(1:1000));% 峰值检测[peaks,locs]=findpeaks(ecg,'MinPeakHeight',threshold,...'MinPeakDistance',round(0.2*fs),'Prominence',0.5*threshold);end3. RR间期提取
functionrr_intervals=get_rr_intervals(locs,fs)rr_samples=diff(locs);rr_intervals=rr_samples/fs*1000;% 转换为毫秒% 异常值处理（3σ原则）mean_rr=mean(rr_intervals);std_rr=std(rr_intervals);rr_intervals(rr_intervals<mean_rr-3*std_rr|rr_intervals>mean_rr+3*std_rr)=[];end4. 时域分析
function[sdnn,rmssd,pnn50]=time_domain(rr_intervals)N=length(rr_intervals);% 基本指标sdnn=std(rr_intervals);% 差分序列diffs=diff(rr_intervals);% RMSSDrmssd=sqrt(mean(diffs.^2));% pNN50nn50=sum(diffs>50);pnn50=nn50/N;end5. 频域分析（Welch方法）
function[lf_power,hf_power,lf_hf_ratio]=freq_domain(rr_intervals,fs)N=length(rr_intervals);window=hamming(256);noverlap=128;% 功率谱估计[Pxx,f]=pwelch(rr_intervals,window,noverlap,[],fs);% 频段划分lf_band=[0.04,0.15];hf_band=[0.15,0.4];% 频段积分lf_idx=f>=lf_band(1)&f<=lf_band(2);hf_idx=f>hf_band(1)&f<=hf_band(2);lf_power=trapz(f(lf_idx),Pxx(lf_idx));hf_power=trapz(f(hf_idx),Pxx(hf_idx));% LF/HF比值lf_hf_ratio=lf_power/hf_power;end6. 可视化模块
functionplot_hrv_results(ecg,rr_intervals,locs,sdnn,rmssd,lf_hf_ratio)figure;% ECG波形与R峰标记subplot(3,1,1);plot(ecg);hold on;plot(locs,ecg(locs),'ro');title('ECG波形与R峰检测');xlabel('时间 (s)');ylabel('幅值 (mV)');% RR间期序列subplot(3,1,2);stem(rr_intervals,'b','LineWidth',1.5);title('RR间期序列');xlabel('心跳序号');ylabel('间隔 (ms)');% 时频域指标subplot(3,1,3);bar([sdnn,rmssd,lf_hf_ratio]);set(gca,'XTickLabel',{'SDNN (ms)','RMSSD (ms)','LF/HF'});ylabel('数值');title('HRV时频域指标');end三、完整调用示例
%% 参数设置fs=250;% 采样率 (Hz)ecg=load('ecg_signal.mat');% 加载ECG数据（列向量）%% 预处理clean_ecg=preprocess_ecg(ecg,fs);%% R波检测[r_peaks,locs]=detect_r_peaks(clean_ecg,fs);%% RR间期提取rr_intervals=get_rr_intervals(locs,fs);%% 时域分析[sdnn,rmssd,pnn50]=time_domain(rr_intervals);%% 频域分析[lf_power,hf_power,lf_hf_ratio]=freq_domain(rr_intervals,fs);%% 可视化plot_hrv_results(clean_ecg,rr_intervals,locs,sdnn,rmssd,lf_hf_ratio);四、结果示例
| 指标 | 正常范围 | 临床意义 |
|---|---|---|
| SDNN | 50-100 ms | 总体自主神经活性 |
| RMSSD | 20-80 ms | 副交感神经张力 |
| LF/HF | 0.5-2.0 | 交感-副交感平衡 |
参考代码 采用ECG信号进行HRV的时域分析、频域分析的Matlab代码www.youwenfan.com/contentcnq/55064.html
五、扩展功能
呼吸耦合分析
% 呼吸信号同步处理[resp_rate,phase]=respiratory_analysis(ecg,fs);非线性分析
% Poincaré图[sd1,sd2]=poincare_plot(rr_intervals);机器学习集成
% 使用HRV特征训练分类模型features=[sdnn,rmssd,lf_hf_ratio];model=train_classifier(features,labels);
六、注意事项
数据质量验证：需检查QRS波群形态（使用ECG诊断工具箱）
采样率要求：建议≥250Hz（MIT-BIH标准）
临床验证：需与专业医疗设备（如Holter）对比
参考文献：
ECG信号预处理与R波检测算法优化
HRV时域频域特征计算方法
基于Welch方法的功率谱估计
运动伪影消除技术