HRV解析用Rスクリプト。xlsxファイルから読み込み。要RHRV, XLConnect

hrv_analysis_xlsx.R

# ## データ読み込み
# 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)
}

サンプルスクリプト2. analyzeTimeSeriesLfhf： LF/HFとかの時系列解析

要 RHRV, XLConnect

上のサンプルのstrRRを利用のこと

t_duration (例60, 180とか)とかt_step (10とか100とか)をいじってどう変わるか観察してみましょう。

# 設定
t_duration <- 120 # LF/HFの解析に使う区間の時間(sec)、60秒だと1/60 Hz (約0.0167 Hz)刻みということになる。
t_step <- 1 # 上の解析区間(窓とも言う)をどの何秒ずつずらして解析していくか(単位sec)。
doPlot <- T # True (真)の意味プロットが面倒な時はF(False: 偽)にしましょう

# 解析
hrv_freq <- analyzeTimeSeriesLfhf(strRR, t_duration, t_step, doPlot)

# csv書き込み：必要に応じてsetwd(dirname(path_xlsx))などで書き込みするディレクトリを指定のこと。
## getwd()で保存しようとしているディレクトリ(カレントディレクトリ)がわかります。
hrv_table <- cbind(hrv_freq$Time, hrv_freq$LFHF, hrv_freq$HF)
colnames(hrv_table) <- c('Time','LF/HF', 'HF')
write.table(file='hrv_lfhf_time_series.csv',hrv_table,sep=',',row.name=F, col.name=T)

サンプルスクリプト3. Periodgram

上と同様ですが、 extract_rrdata_from_xlsx_end_plus1min160621で出力されるstrRRを利用します (発作時刻のセルの位置が変わったため).

解析区間は発作開始時刻から以前(例えば60秒とか)と、発作終了後(例えば60秒とか)の2箇所解析してプロットするようにしてます.

## 設定
pre_sec <- 60 # sec ：発作時刻を基準に前の解析区間
pst_sec <- 60 # sec: 発作終了時刻を基準に後の解析区間
doPlot <- T # プロットをしたくない時はFにする

## データ読み込み：ファイル1個の時
path_xlsx <- 'D:/directory/to/data/rr_data.xlsx'
strRR <- extract_rrdata_from_xlsx_end_plus1min160621(path_xlsx, pre_sec, pst_sec)

## 解析: doPlot =Tの時(数行上参照) プロットが出力されます。なんとなく意味はわかると期待.
spc_out <- GetPeriodogram(strRR, pre_sec, pst_sec, doPlot)

# output:
#     spc_out$pre: 発作時刻を基準に前の解析結果
#     spc_out$pst: 発作終了時刻を基準に後の解析結果

サンプルスクリプト3の続き：(2) CSVへの出力

CSVに出力するときは下記を追加

## csvに出力する時はspc_outを利用します:
# 設定
fname_pre <- 'periodogram_pre.csv'
fname_pst <- 'periodogram_pst.csv'

# 以下基本変更不要
colname_periodo <- c('freq','spec')
data_pre <- cbind(spc_out$pre$freq, spc_out$pre$spec)
data_pre <- rbind(colname_periodo, data_pre)
data_pst <- cbind(spc_out$pst$freq, spc_out$pst$spec)
data_pst <- rbind(colname_periodo, data_pst)

write.table(file=fname_pre, data_pre, sep=',', row.names=F, col.names=F)
write.table(file=fname_pst, data_pst, sep=',', row.names=F, col.names=F)

サンプルスクリプト3 続き(3): 複数のxlsxファイルをまとめて処理する.

複数のxlsxがあるディレクトリを指定して、そこにあるファイルをまとめて処理します。中央値のプロットを出力します。
medianとあるところをmeanにすれば平均値になります。

"## 設定"としているところが2箇所あります。その中を適宜確認、変更すると作業できると思います。

最後のほうでCSVの出力をしています。

## 設定
dir_xlsx <- '/dir/to/xlsx/files' # .xlsxファイルをまとめているディレクトリ
pre_sec <- 60 # 発作開始時刻を基準に取り出す発作前解析区間(sec)
pst_sec <- 60 # 発作終了時刻を基準に取り出す発作後解析区間(sec)
doPlot <- T # TRUE(真)の意味。プロットを出したくないときはF(FALSE, 偽)。

## 解析
setwd(dir_xlsx)
files_xlsx <- dir(pattern='*.xlsx')
ii <- 0
for (fname in files_xlsx){
  strRR <- extract_rrdata_from_xlsx_end_plus1min160621(fname, pre_sec, pst_sec)
  spc_out <- GetPeriodogram(strRR, pre_sec, pst_sec, doPlot)
  if (ii == 0){
    mat_spc_pre_out <- spc_out$pre$spec
    mat_spc_pst_out <- spc_out$pst$spec
  } else {
    mat_spc_pre_out <- cbind(mat_spc_pre_out, spc_out$pre$spec)
    mat_spc_pst_out <- cbind(mat_spc_pst_out, spc_out$pst$spec)
  }
  ii = ii + 1
}

freq_pre <- spc_out$pre$freq
freq_pst <- spc_out$pst$freq

## csv出力
# 設定: ファイル名
fname_csv_pre <- 'periodogram_pre.csv' # ディレクトリを変えたいときはディレクトリを付加するか.
fname_csv_pst <- 'periodogram_pst.csv' #   この前にsetwd(hogehoge)でディレクトリを変える

# freqとspecまとめる
data_out_pre <- cbind(freq_pre, mat_spc_pre_out)
data_out_pst <- cbind(freq_pst, mat_spc_pst_out)

# 出力
write.table(file = fname_csv_pre, data_out_pre, sep=',', col.names=F, row.names=F)
write.table(file = fname_csv_pst, data_out_pst, sep=',', col.names=F, row.names=F)

サンプルスクリプト3続き (4)：上の複数ファイルのデータを更にまとめて(中央値)プロット

r_figure() # windows()と同じ. os非依存でwindowを立ち上げる
par(bty='l')
par(las=1)

# 下のmedianをmeanにすると平均値
plot(freq_pre, apply(mat_spc_pre_out, 1, median), log='y',type='l', 
     col='blue',xlab='frequency (Hz)', ylab='spectrum', main='Raw Periodogram')
lines(freq_pst, apply(mat_spc_pst_out, 1, median), col='red')
legend('topright',legend=c('pre','post'), col=c('blue','red'),lty=1)