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) | |
} |
サンプルスクリプト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)
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
サンプルスクリプト3. Periodgram
上と同様ですが、 extract_rrdata_from_xlsx_end_plus1min160621で出力されるstrRRを利用します (発作時刻のセルの位置が変わったため).
解析区間は発作開始時刻から以前(例えば60秒とか)と、発作終了後(例えば60秒とか)の2箇所解析してプロットするようにしてます.