本帖最后由 天路 于 2018-2-25 21:16 编辑

本人正在学习心率变异性方面的内容，但是按照文献上的方法做出来的结果并不是很理想，文献上说的是心率变异性的频率的范围是0.4以内，但是我做的功率谱上显示频率分布在整个频域内，试了很多种方法也没有效果。不知道怎么回事，下面是我的代码，请各位帮个忙吧，这个问题困扰很长时间了.

clc

clear

load data

data1=data(260000:360000);

%% 去除肌电信号干扰

Fs=1000;                        %采样频率

fp=10;fs=12;                    %通带截止频率，阻带截止频率

rp=1.4;rs=5;                    %通带、阻带衰减

wp=2*pi*fp;ws=2*pi*fs;

[n,wn]=buttord(wp,ws,rp,rs,'s');     %'s'是确定巴特沃斯模拟滤波器阶次和3dB截止模拟频率

Fs=500;             %采样频率

fp=20;fs=25; wp=2*pi*fp;ws=2*pi*fs; rp=1.4;rs=1.6; %通带截止频率，阻带截止频率 通带、阻带衰减

[n,wn]=buttord(wp,ws,rp,rs,'s');     %'s'是确定巴特沃斯模拟滤波器阶次和3dB截止模拟频率

[z,P,k]=buttap(n);   %设计归一化巴特沃斯模拟低通滤波器，z为极点，p为零点和k为增益

[bp,ap]=zp2tf(z,P,k);  %转换为Ha(p),bp为分子系数，ap为分母系数

[bs,as]=lp2lp(bp,ap,wp) ;%Ha(p)转换为低通Ha(s)并去归一化，bs为分子系数，as为分母系数

[bz,az]=bilinear(bs,as,Fs);     %对模拟滤波器双线性变换

data_lvbo=filter(bz,az,data1);

%% RR间期序列提取

data_qjx=qjxpj1(data_lvbo);  %去除基线漂移

[data_Rf,t]=RFtiqu(data_qjx);  %提取R峰

s=1:length(data1);

data_rr=diff(t)/Fs; %生成RR间期序列

figure

plot (s,data_qjx,t,data_Rf,'r+') %画R峰图

data_rr(find(data_rr>1.3*mean(data_rr)))=[];

data_rr(find(data_rr<0.7*mean(data_rr)))=[];

t(find(data_rr>1.3*mean(data_rr)))=[];

t(find(data_rr<0.7*mean(data_rr)))=[];

%%

%均匀重采样

fs_resample=10;%重采样频率

f_RR_resample=t(1):Fs/fs_resample:t(end-1);%确定插值频率

t(end)=[];

RR_interp=interp1(t,data_rr,f_RR_resample,'spline');

% RR_interp=interp1(t,real_RR,f_RR_resample,'spline');%插值，重采样

figure;

plot(f_RR_resample/fs,RR_interp,'.');

% RR_interp=RR_interp-mean(RR_interp);   %该语句的实质作用？去除直流分量

RR_interp=QZL(RR_interp);

%%窗函数法计算功率谱并画图

nfft=2^(nextpow2(length(RR_interp)-1));

window=hamming(length(RR_interp));

[Pxx,f]=pwelch(RR_interp,window,0,nfft,fs_resample,'onesided');

Pxx=Pxx*2;

figure;

plot(f,Pxx);%实信号，频谱有正负 ，乘以2是真实功率谱

xlabel('f/Hz');ylabel('功率');title('瞬时心率的功率谱')

[px3,f3]=pyulear(data_rr,16,nfft,10);

figure

plot(f3,px3);title('AR MODEL')

下面附计算结果的图片

AR模型功率谱法.jpg

(34.52 KB, 下载次数: 3)

2018-2-25 16:56 上传

AR模型求功率谱

R峰监测.jpg

(51.4 KB, 下载次数: 0)

2018-2-25 16:56 上传

R峰监测

窗函数法功率谱.jpg

(45.86 KB, 下载次数: 0)

2018-2-25 16:59 上传

汉明窗求功率谱