Last active
July 5, 2016 12:58
-
-
Save cosacog/20ccaf6a1f080b6105f3 to your computer and use it in GitHub Desktop.
HRV解析用Rスクリプト。xlsxファイルから読み込み。要RHRV, XLConnect
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# ## データ読み込み | |
# setwd('C:/directory/to/xlsx/file') | |
# rrdata_pre<-read.table('before.txt',sep='/t',header=F) | |
# rrdata_pst<-read.table('after.txt',sep='/t',header=F) | |
# time_onset <- sum(rrdata_pre)/1000 # ここで前後を区切る | |
# doPlot<-TRUE | |
# names(rrdata_pre) <- names(rrdata_pst) <- 'Time' | |
# rrdata <- rbind(rrdata_pre,rrdata_pst) | |
## define function | |
r_figure <- function(){ | |
# osに依存せずウィンドウを出す | |
os_name <- tolower(.Platform$OS.type) | |
switch(os_name, | |
windows = {windows()}, | |
unix = {X11()} | |
) | |
} | |
calculate_hrv_strRR <- function(strRR,doPlot){ | |
## calculate_hrv修正 | |
## strRRからデータを取り出す | |
rrdata <- strRR$rrdata | |
time_onset <- strRR$time_onset | |
name_str_rr <- names(strRR) | |
is_t_ofset_exist <- !is.na(any(charmatch(name_str_rr,'time_ofset'))) | |
if (is_t_ofset_exist){ | |
time_ofset <- strRR$time_ofset | |
} | |
library(RHRV) | |
hrv_data <- CreateHRVData() | |
hrv_data$Verbose <- T | |
hrv_data$Beat <- rbind(0,cumsum(rrdata))/1000 | |
# heart rate data | |
hrv_dataRR <- BuildNIHR(hrv_data) | |
# filter | |
hrvRrFilt <- FilterNIHR(hrv_dataRR) | |
# interpolation | |
freq_hr =4 | |
hrvRrFI <- InterpolateNIHR(hrvRrFilt,freqhr = freq_hr) | |
# heart rate at onset time | |
hr_intrp <- hrvRrFI$HR | |
time_vector <- seq(0,by=1/freq_hr, length=length(hr_intrp)) | |
idx_onset <- min(which(time_vector>time_onset)) | |
hr_onset <- hr_intrp[idx_onset] | |
# separate hr data | |
# pre: 発作??秒前~発作開始時刻 | |
hrvRrFIPre <- hrvRrFI | |
hrvRrFIPre$Beat <- hrvRrFI$Beat[hrvRrFI$Beat$Time<=time_onset,] | |
hrvRrFIPre$HR <- hrvRrFI$HR[1:idx_onset] | |
# post: 発作開始時刻~発作終了**後 | |
hrvRrFIPst <- hrvRrFI | |
hrvRrFIPst$Beat <- hrvRrFI$Beat[hrvRrFI$Beat$Time>time_onset,] | |
hrvRrFIPst$HR <- hrvRrFI$HR[idx_onset:length(time_vector)] | |
# post offset: 発作終了時刻~発作終了**後 | |
if (is_t_ofset_exist){ | |
hrvRrFIPstOfset <- hrvRrFI | |
hrvRrFIPstOfset$Beat <- hrvRrFI$Beat[hrvRrFI$Beat$Time>time_ofset,] | |
idx_ofset <- min(which(time_vector>time_ofset)) | |
hrvRrFIPstOfset$HR <- hrvRrFI$HR[idx_ofset:length(time_vector)] | |
} | |
## plot | |
# r_figure() | |
# par(las=1) | |
# PlotNIHR(hrv_dataRR) | |
if (doPlot){ | |
r_figure() | |
par(las=1) | |
PlotHR(hrvRrFI) | |
# arrows(time_onset,hr_onset + 10,time_onset,hr_onset +5, | |
# code=2,angle=45, length=0.1,col='red') | |
abline(v=time_onset,col='red') | |
if (is_t_ofset_exist){ | |
abline(v= time_ofset,col='blue') | |
} | |
} | |
# r_figure() | |
# PlotNIHR(hrvRrFI) | |
################### analysis ##################################### | |
# time analysis | |
hrvRrFITanaPre <- CreateTimeAnalysis(hrvRrFIPre, size=time_onset) | |
t_size_pst <- diff(range(time_vector[(idx_onset+1):length(time_vector)])) | |
hrvRrFITanaPst <- CreateTimeAnalysis(hrvRrFIPst, size=t_size_pst) | |
list_hrvRrFITana <- list(hrvRrFITanaPre,hrvRrFITanaPst) | |
t_size <- c(time_onset,t_size_pst) | |
if (is_t_ofset_exist){ | |
t_size_pst_ofset <- diff(range(time_vector[(idx_ofset+1):length(time_vector)])) | |
hrvRrFITanaPstOfset <- CreateTimeAnalysis(hrvRrFIPstOfset, size=t_size_pst_ofset) | |
list_hrvRrFITana <-c(list_hrvRrFITana, list(hrvRrFITanaPstOfset)) | |
t_size <- c(t_size, t_size_pst_ofset) | |
} | |
# frequency analysis:settings | |
ULFmin=0; ULFmax= 0.03 | |
VLFmin=ULFmax; VLFmax= 0.05 | |
LFmin =VLFmax; LFmax = 0.15 | |
HFmin =LFmax; HFmax = 0.4 | |
list_hrvRrFITFana <- as.list(NULL) # frequency analysis | |
list_hrvRrFITFanaNl<-as.list(NULL) # poincare plot, approximate entropy | |
for (ii in c(1:length(list_hrvRrFITana))){ | |
# frequency analysis | |
list_hrvRrFITFana[[ii]] <- CreateFreqAnalysis(list_hrvRrFITana[[ii]]) | |
list_hrvRrFITFana[[ii]] <- CalculatePowerBand(list_hrvRrFITFana[[ii]], | |
indexFreqAnalysis =1, | |
size = t_size[ii], | |
shift = 1, | |
type='fourier', | |
ULFmin=ULFmin, ULFmax=ULFmax, | |
VLFmin=VLFmin, VLFmax=VLFmax, | |
LFmin=LFmin, LFmax=LFmax, | |
HFmin=HFmin, HFmax=HFmax) | |
# 結果参照 | |
# list_hrvRrFITFana[[1]]$FreqAnalysis # 1~3まで | |
# poincare plot | |
list_hrvRrFITFanaNl[[ii]] <- CreateNonLinearAnalysis(list_hrvRrFITFana[[ii]]) | |
# bugfix(niHRを使ってSD1,SD2を計算している) | |
list_hrvRrFITFanaNl[[ii]]$Beat$HR <-list_hrvRrFITFanaNl[[ii]]$Beat$niHR | |
# list_hrvRrFITFanaNl[[ii]]$Beat$niHR <-list_hrvRrFITFanaNl[[ii]]$Beat$RR | |
if (doPlot){r_figure()} | |
list_hrvRrFITFanaNl[[ii]] <- PoincarePlot(list_hrvRrFITFanaNl[[ii]], | |
indexNonLinearAnalysis = 1, | |
timeLag=1, doPlot=doPlot) | |
# approximate entropy: warning- deprecated function. CalculateSampleEntropy recommended instead | |
list_hrvRrFITFanaNl[[ii]] <- CalculateApEn(list_hrvRrFITFanaNl[[ii]], | |
indexNonLinearAnalysis = length(list_hrvRrFITFanaNl[[ii]]$NonLinearAnalysis), | |
m = 2, tau = 1, | |
r = 0.2, N = 1000, verbose=NULL) | |
} | |
# データまとめ:pre,post, post_ofset | |
df_result <- NULL | |
for (ii in c(1:length(list_hrvRrFITana))){ | |
# TimeAnalysis | |
hrvTrslt <- unlist(list_hrvRrFITFanaNl[[ii]]$TimeAnalysis[[1]]) | |
hrvTrslt3 <- hrvTrslt[c(2,7,5)] | |
# FreqAnalysis | |
hrvFrslt <- unlist(list_hrvRrFITFanaNl[[ii]]$FreqAnalysis[[1]]) | |
hrvFrslt3 <- hrvFrslt[c(4:6)] | |
# poincare plot | |
hrvNlrslt <- unlist(list_hrvRrFITFanaNl[[ii]]$NonLinearAnalysis[[1]]$PoincarePlot) | |
# SD1SD2 | |
sd1sd2 <- hrvNlrslt[1]/hrvNlrslt[2] | |
names(sd1sd2)<- 'SD1SD2' | |
# approximate entropy | |
ApEn <- unlist(list_hrvRrFITFanaNl[[ii]]$NonLinearAnalysis[[1]]$ApEn) | |
names(ApEn) <- 'ApEn' | |
# まとめ | |
hrvRslt <- c(hrvTrslt3, hrvFrslt3, sd1sd2, hrvNlrslt,ApEn) | |
print(hrvRslt) | |
df_result <- rbind(df_result, hrvRslt) | |
} | |
# row name変更 | |
row_names <- c('pre','pst','pst_offset') | |
row.names(df_result)<-row_names[c(1:length(list_hrvRrFITana))] | |
df_result <- as.data.frame(df_result) | |
return(df_result) | |
} | |
extract_rrdata_from_xlsx <- function(path_xlsx,pre_sec,pst_sec){ | |
# 前提条件 | |
# 1列目(6行目以降)が時刻 | |
# 2列目(6行目以降)がrr間隔(ms) | |
# 3行3列に発作開始時刻 | |
library(XLConnect) | |
rrdata_xlsx<-readWorksheet(loadWorkbook(path_xlsx), | |
sheet=1,header=F) | |
vec_rr <- as.numeric(rrdata_xlsx[6:nrow(rrdata_xlsx),2]) | |
# 時刻データ | |
xlsx_onset_str <- strsplit(rrdata_xlsx[3,3],' ')[[1]][2] | |
xlsx_onset<-as.numeric(strptime(xlsx_onset_str,'%H:%M:%S')) | |
xlsx_time <- rrdata_xlsx[6:nrow(rrdata_xlsx),1] | |
xlsx_time_split <- strsplit(xlsx_time,' ') | |
vec_time <- NULL # 時刻ベクトル | |
for (ii in xlsx_time_split){vec_time <- c(vec_time, | |
as.numeric(strptime(ii[2],'%H:%M:%S')))} | |
# 区切りのindex | |
idx_onset <- min(which(vec_time==xlsx_onset)) | |
#vec_rr_pre_cumsum <- rev(cumsum(rev(vec_rr[1:idx_onset]))-vec_rr[idx_onset]) | |
vec_rr_pre_cumsum <- rev(cumsum(rev(vec_rr[1:idx_onset]))) | |
idx_start <- max(which(vec_rr_pre_cumsum > pre_sec*1000)) | |
vec_rr_pst_cumsum <- cumsum(vec_rr[idx_onset:length(vec_rr)])-vec_rr[idx_onset] | |
idx_end <- min(which(vec_rr_pst_cumsum> pst_sec*1000)) | |
# rrデータ取り出し | |
rrdata <- as.data.frame(matrix(vec_rr[idx_start:(idx_onset+idx_end)],ncol=1)) | |
colnames(rrdata)<-'Time' | |
time_onset <- sum(vec_rr[(idx_start):idx_onset])/1000 | |
str<-NULL | |
str$rrdata <- rrdata | |
str$time_onset <- time_onset | |
return(str) | |
} | |
extract_rrdata_from_xlsx_end_plus1min <- function(path_xlsx,pre_sec,pst_sec){ | |
# 前提条件 | |
# 1列目(6行目以降)が時刻 | |
# 2列目(6行目以降)がrr間隔(ms) | |
# 3行3列に発作開始時刻 | |
# 3行6列に発作終了時刻 | |
# 発作終了+pst_sec(通常60sec)まで取り出し | |
library(XLConnect) | |
rrdata_xlsx<-readWorksheet(loadWorkbook(path_xlsx), | |
sheet=1,header=F) | |
vec_rr <- as.numeric(rrdata_xlsx[6:nrow(rrdata_xlsx),2]) | |
# 時刻データ | |
# onset/offset | |
xlsx_onset_str <- strsplit(as.character(rrdata_xlsx[3,3]),' ')[[1]][2] | |
xlsx_ofset_str <- strsplit(as.character(rrdata_xlsx[3,6]),' ')[[1]][2] | |
xlsx_onset<-as.numeric(strptime(xlsx_onset_str,'%H:%M:%S'))# 処理した日の日付が含まれる | |
xlsx_ofset<-as.numeric(strptime(xlsx_ofset_str,'%H:%M:%S')) | |
xlsx_time <- rrdata_xlsx[6:nrow(rrdata_xlsx),1] | |
xlsx_time_split <- strsplit(xlsx_time,' ') | |
vec_time <- NULL # 時刻ベクトル | |
for (ii in xlsx_time_split){vec_time <- c(vec_time, | |
as.numeric(strptime(ii[2],'%H:%M:%S')))} | |
# 区切りのindex | |
idx_onset <- min(which(vec_time==xlsx_onset)) | |
idx_ofset <- min(which(vec_time==xlsx_ofset)) | |
#vec_rr_pre_cumsum <- rev(cumsum(rev(vec_rr[1:idx_onset]))-vec_rr[idx_onset]) | |
vec_rr_pre_cumsum <- rev(cumsum(rev(vec_rr[1:idx_onset]))) | |
idx_start <- max(which(vec_rr_pre_cumsum > pre_sec*1000)) | |
#vec_rr_pst_cumsum <- cumsum(vec_rr[idx_onset:length(vec_rr)])-vec_rr[idx_onset] | |
vec_rr_pst_cumsum <- cumsum(vec_rr[idx_ofset:length(vec_rr)])-vec_rr[idx_ofset] | |
idx_end <- min(which(vec_rr_pst_cumsum> pst_sec*1000)) | |
# rrデータ取り出し | |
#rrdata <- as.data.frame(matrix(vec_rr[idx_start:(idx_onset+idx_end)],ncol=1)) | |
rrdata <- as.data.frame(matrix(vec_rr[idx_start:(idx_ofset+idx_end)],ncol=1)) | |
colnames(rrdata)<-'Time' | |
time_onset <- sum(vec_rr[(idx_start):idx_onset])/1000 | |
time_ofset <- sum(vec_rr[(idx_start):idx_ofset])/1000 | |
str<-NULL | |
str$rrdata <- rrdata | |
str$time_onset <- time_onset | |
str$time_ofset <- time_ofset | |
return(str) | |
} | |
extract_rrdata_from_xlsx_end_plus1min160621 <- function(path_xlsx,pre_sec,pst_sec){ | |
# 前提条件 | |
# 1列目(6行目以降)が時刻 => 5行目以降 | |
# 2列目(6行目以降)がrr間隔(ms) => 5行目以降 | |
# 3行3列に発作開始時刻 => 1行2列目 | |
# 3行6列に発作終了時刻 => 2行2列目 | |
# 発作終了+pst_sec(通常60sec)まで取り出し | |
library(XLConnect) | |
rrdata_xlsx<-readWorksheet(loadWorkbook(path_xlsx), | |
sheet=1,header=F) | |
# vec_rr <- as.numeric(rrdata_xlsx[6:nrow(rrdata_xlsx),2]) | |
vec_rr <- as.numeric(rrdata_xlsx[5:nrow(rrdata_xlsx),2]) | |
# 時刻データ | |
# onset/offset | |
# xlsx_onset_str <- strsplit(as.character(rrdata_xlsx[3,3]),' ')[[1]][2] | |
# xlsx_ofset_str <- strsplit(as.character(rrdata_xlsx[3,6]),' ')[[1]][2] | |
xlsx_onset_str <- strsplit(as.character(rrdata_xlsx[1,2]),' ')[[1]][2] | |
xlsx_ofset_str <- strsplit(as.character(rrdata_xlsx[2,2]),' ')[[1]][2] | |
xlsx_onset<-as.numeric(strptime(xlsx_onset_str,'%H:%M:%S'))# 処理した日の日付が含まれる | |
xlsx_ofset<-as.numeric(strptime(xlsx_ofset_str,'%H:%M:%S')) | |
# xlsx_time <- rrdata_xlsx[6:nrow(rrdata_xlsx),1] | |
xlsx_time <- rrdata_xlsx[5:nrow(rrdata_xlsx),1] | |
xlsx_time_split <- strsplit(xlsx_time,' ') | |
vec_time <- NULL # 時刻ベクトル | |
for (ii in xlsx_time_split){vec_time <- c(vec_time, | |
as.numeric(strptime(ii[2],'%H:%M:%S')))} | |
# 区切りのindex | |
idx_onset <- min(which(vec_time==xlsx_onset)) | |
idx_ofset <- min(which(vec_time==xlsx_ofset)) | |
#vec_rr_pre_cumsum <- rev(cumsum(rev(vec_rr[1:idx_onset]))-vec_rr[idx_onset]) | |
vec_rr_pre_cumsum <- rev(cumsum(rev(vec_rr[1:idx_onset]))) | |
idx_start <- max(which(vec_rr_pre_cumsum > pre_sec*1000)) | |
#vec_rr_pst_cumsum <- cumsum(vec_rr[idx_onset:length(vec_rr)])-vec_rr[idx_onset] | |
vec_rr_pst_cumsum <- cumsum(vec_rr[idx_ofset:length(vec_rr)])-vec_rr[idx_ofset] | |
idx_end <- min(which(vec_rr_pst_cumsum> pst_sec*1000)) | |
# rrデータ取り出し | |
#rrdata <- as.data.frame(matrix(vec_rr[idx_start:(idx_onset+idx_end)],ncol=1)) | |
rrdata <- as.data.frame(matrix(vec_rr[idx_start:(idx_ofset+idx_end)],ncol=1)) | |
colnames(rrdata)<-'Time' | |
time_onset <- sum(vec_rr[(idx_start):idx_onset])/1000 | |
time_ofset <- sum(vec_rr[(idx_start):idx_ofset])/1000 | |
str<-NULL | |
str$rrdata <- rrdata | |
str$time_onset <- time_onset | |
str$time_ofset <- time_ofset | |
return(str) | |
} | |
### hf, lf/hf 時系列解析 ################# | |
analyzeTimeSeriesLfhf <- function(strRR, t_duration, t_step, doPlot){ | |
# HF, LF/HFを時系列で計算. 要RHRV | |
# @param ... strRR, t_duration, t_step, doPlot | |
# @param strRR: stract_rrdata_from_xlsx_end_plus1minによるデータ | |
# @param t_duration: 解析時間幅(sec) | |
# @param t_step: 解析ステップ (sec) | |
# @param doPlot: boolean, プロットするかどうか | |
# | |
# @return strRtrn: LF/HF, HFとかのデータと時刻、その周波数の設定 | |
# データ取り込み | |
rrdata_orig <- strRR$rrdata | |
time_onset <- strRR$time_onset | |
name_str_rr <- names(strRR) | |
is_t_ofset_exist <- !is.na(any(charmatch(name_str_rr,'time_ofset'))) | |
if (is_t_ofset_exist){ | |
time_ofset <- strRR$time_ofset | |
} | |
library(RHRV) | |
rr_data <- CreateHRVData() | |
rr_data$Verbose <- T | |
rr_data$Beat <- rbind(0,cumsum(rrdata_orig))/1000 | |
# heart rate data | |
hr_rr_data <- BuildNIHR(rr_data) | |
# filter | |
hr_rr_filt <- FilterNIHR(hr_rr_data) | |
# interpolation | |
freq_hr =4 # max 60*freq_hr beat/minということ | |
hr_rr_filt_intrp <- InterpolateNIHR(hr_rr_filt,freqhr = freq_hr) | |
# heart rate at onset time | |
hr_intrp <- hr_rr_filt_intrp$HR | |
time_vector <- seq(0,by=1/freq_hr, length=length(hr_intrp)) | |
# frequency analysis:settings | |
ULFmin=0; ULFmax= 0.03 | |
VLFmin=ULFmax; VLFmax= 0.05 | |
LFmin =VLFmax; LFmax = 0.15 | |
HFmin =LFmax; HFmax = 0.4 | |
hrv_Freq <- CreateFreqAnalysis(hr_rr_filt_intrp) | |
hrv_Freq <- CalculatePowerBand(hrv_Freq, | |
indexFreqAnalysis =1, | |
size = t_duration, | |
shift = t_step, | |
type='fourier', | |
ULFmin=ULFmin, ULFmax=ULFmax, | |
VLFmin=VLFmin, VLFmax=VLFmax, | |
LFmin=LFmin, LFmax=LFmax, | |
HFmin=HFmin, HFmax=HFmax) | |
# データをまとめる | |
# print(names(hrv_Freq$FreqAnalysis[[1]])) # 'Time'が生成されないことがある | |
# print(!any(names(hrv_Freq$FreqAnalysis[[1]])=='Time')) | |
strRtrn <- NULL | |
if (!any(names(hrv_Freq$FreqAnalysis[[1]])=='Time')){ | |
t_xaxis <- seq(from = t_duration/2, by = t_step, | |
length=length(hrv_Freq$FreqAnalysis[[1]]$HF)) | |
strRtrn$Time <- t_xaxis | |
} else { | |
strRtrn$Time <- hrv_Freq$FreqAnalysis[[1]]$Time | |
} | |
strRtrn$LFHF <- hrv_Freq$FreqAnalysis[[1]]$LFHF | |
strRtrn$HF <- hrv_Freq$FreqAnalysis[[1]]$HF | |
strRtrn$LF <- hrv_Freq$FreqAnalysis[[1]]$LF | |
strRtrn$VLF <- hrv_Freq$FreqAnalysis[[1]]$VLF | |
strRtrn$ULF <- hrv_Freq$FreqAnalysis[[1]]$ULF | |
strRtrn$t_size_analysis <- hrv_Freq$FreqAnalysis[[1]]$size | |
strRtrn$t_shift <- hrv_Freq$FreqAnalysis[[1]]$shift | |
strRtrn$FreqHF <- c(hrv_Freq$FreqAnalysis[[1]]$HFmin, hrv_Freq$FreqAnalysis[[1]]$HFmax) | |
strRtrn$FreqLF <- c(hrv_Freq$FreqAnalysis[[1]]$LFmin, hrv_Freq$FreqAnalysis[[1]]$LFmax) | |
strRtrn$FreqVLF <- c(hrv_Freq$FreqAnalysis[[1]]$VLFmin,hrv_Freq$FreqAnalysis[[1]]$VLFmax) | |
strRtrn$FreqULF <- c(hrv_Freq$FreqAnalysis[[1]]$ULFmin,hrv_Freq$FreqAnalysis[[1]]$ULFmax) | |
# plot ---------------------------------------------------- | |
xlim_time <- range(hr_rr_filt_intrp$Beat$Time) | |
if (doPlot){ | |
r_figure() | |
par(mfrow=c(3,1)) | |
par(plt=c(0.1,0.9,0.2,0.8)) | |
par(las=1) | |
PlotHR(hr_rr_filt_intrp) | |
# arrows(time_onset,hr_onset + 10,time_onset,hr_onset +5, | |
# code=2,angle=45, length=0.1,col='red') | |
abline(v=time_onset,col='red') | |
if (is_t_ofset_exist){ | |
abline(v= time_ofset,col='blue') | |
} | |
plot(strRtrn$Time, strRtrn$LFHF, type='o', xlim=xlim_time, main='LF/HF', | |
xlab='',ylab='LF/HF') | |
plot(strRtrn$Time, strRtrn$HF, type='o', xlim=xlim_time, main='HF', | |
xlab='time (sec)', ylab='HF') | |
} | |
return(strRtrn) | |
} | |
## periodgramを返す ## | |
GetPeriodogram <- function(strRR, durPre, durPst, doPlot){ | |
## periodgramを返す | |
# params: | |
# strRR:extract_rrdata_from_xlsx_end_plus1minでできるデータ | |
# durPre: 発作前の解析区間(s) | |
# durPst: 発作後の解析区間(s) | |
# return: | |
# spc_out: spectrumの解析結果 | |
# spc_out$pre: 発作前の解析結果 | |
# spc_out$pst: 発作後の解析結果 | |
## データ取り込み | |
rrdata_orig <- strRR$rrdata | |
time_onset <- strRR$time_onset | |
name_str_rr <- names(strRR) | |
is_t_ofset_exist <- !is.na(any(charmatch(name_str_rr,'time_ofset'))) | |
if (is_t_ofset_exist){ | |
time_ofset <- strRR$time_ofset | |
} | |
library(RHRV) | |
rr_data <- CreateHRVData() | |
rr_data$Verbose <- T | |
rr_data$Beat <- rbind(0,cumsum(rrdata_orig))/1000 | |
# heart rate data | |
hr_rr_data <- BuildNIHR(rr_data) | |
# filter | |
hr_rr_filt <- FilterNIHR(hr_rr_data) | |
# interpolation | |
freq_hr =4 # max 60*freq_hr beat/minということ | |
hr_rr_filt_intrp <- InterpolateNIHR(hr_rr_filt,freqhr = freq_hr) | |
hr_intrp <- hr_rr_filt_intrp$HR | |
rr_intrp <- 60/hr_intrp * 1000 | |
time_vector <- seq(0,by=1/freq_hr, length=length(hr_intrp)) | |
t_pre_strt <- max(c(0, strRR$time_onset - durPre)) | |
t_pre_end <- strRR$time_onset | |
t_pst_strt <- strRR$time_ofset | |
t_pst_end <- min(c(max(time_vector), t_pst_strt + durPst)) | |
idx_pre <- which((time_vector >= t_pre_strt) & (time_vector <= t_pre_end)) | |
idx_pst <- which(time_vector >= t_pst_strt & time_vector <= t_pst_end) | |
## periodgram | |
spc_out_pre <- spectrum(rr_intrp[idx_pre], method='pgram', plot=F) # arというのもあり | |
spc_out_pst <- spectrum(rr_intrp[idx_pst], method='pgram', plot=F) | |
spc_out <- NULL | |
spc_out$pre <- spc_out_pre | |
spc_out$pst <- spc_out_pst | |
# plot | |
if (doPlot){ | |
r_figure() # windows()と同じ. os非依存でwindowを立ち上げる | |
par(bty='l') | |
par(las=1) | |
par(mfrow=c(2,1)) | |
plot(time_vector, rr_intrp, type='l', xlab='time (s)', ylab='RR interval (ms)', main ='RR interval') | |
abline(v=c(t_pre_strt, t_pre_end), col='blue') | |
abline(v=c(t_pst_strt, t_pst_end), col='red') | |
legend('topright',legend=c('pre解析区間','post解析区間'), col=c('blue','red'),lty=1) | |
plot(spc_out_pre$freq, spc_out_pre$spec, log='y',type='l', | |
col='blue',xlab='frequency (Hz)', ylab='spectrum', main='Raw Periodogram') | |
lines(spc_out_pst$freq, spc_out_pst$spec, col='red') | |
legend('topright',legend=c('pre','post'), col=c('blue','red'),lty=1) | |
} | |
return(spc_out) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
サンプルスクリプト3続き (4):上の複数ファイルのデータを更にまとめて(中央値)プロット