前言

之前讲了MTM（多锥形窗谱估计）的相关原理，现在来分析一下它的R语言的实现，这个实现是提出人的学生写的，和matlab的实现进行对照分析，加深理解，提高大家对这门技术的掌握程度，解析的顺序依旧是从下至上，先从简单的子程序，最后到复杂的主程序。

想要复习的可以参考一下之前的文件：

mtm原理：现代谱估计：多窗口谱
想要复习一下如何实现的可以参考：
MTM：matlab实现1MTM：matlab实现1
MTM：matlab实现2参数解析MTM参数解析
MTM：matlab实现3谱功率计算MTM谱功率计算
MTM：matlab实现4主函数解析MTM 主函数
R学习_multitaper包解析1：子函数centre，dpss：R 子函数

子函数：spec.mtm.dpss

使用slepian窗序列计算多窗口谱估计
代码对应论文在下面
##############################################################
##
##  .spec.mtm.dpss
##
##  Computes multitaper spectrum using Slepian tapers
##  References:
##    Percival and Walden "Spectral Analysis
##    for Physical Applications" 1993 and associated LISP code
##
##    Thomson, D.J. Spectrum Estimation and Harmonic Analysis,
##    Proceedings of the IEEE, 1982 and associated Fortran code
##
##############################################################
输入.spec.mtm.dpss <- function(timeSeries,时间序列nw,窗口阶数k,使用多少个窗口序列nFFT,   使用多少点计算              dpssIN,输入离散扁球序列returnZeroFreq,零频率的幅值需不需要Ftest,F检验jackknife,是否用剪刀法计算数据jkCIProb,t分布的概率值adaptiveWeighting, 自适应的权重值maxAdaptiveIterations,最大自适应的迭代数returnInternals,是否返回间隔n,时间序列长度deltaT,等待间隔sigma2,时间序列的方差series,解析好的时间序列dtUnits,时间间隔单位...) {# Complex check case复值检测if(is.complex(timeSeries)) {如果序列时复值且设定返回零值谱，则调整为返回点1处的谱。并且警告if(!returnZeroFreq) {returnZeroFreq <- 1 warning("Cannot set returnZeroFreq to 0 for complex time series.")} }初始化dw，ev，receivedDW为真dw <- NULLev <- NULLreceivedDW <- TRUE
如果没有初始化dpss序列，则设置为假，使用dpss生成DPssin相关序列参数if(!.is.dpss(dpssIN)) {receivedDW <- FALSEdpssIN <- dpss(n, k, nw=nw, returnEigenvalues=TRUE)dw等于v值乘以v值*采样间隔的均方根 是一个特征值矩阵dw <- dpssIN$v*sqrt(deltaT)eigen是对应的k个特征向量中心能量ev <- dpssIN$eigen }如果初始化了dpss，则将对应值域赋值给相应对象。else {dw <- .dpssV(dpssIN)ev <- .dpssEigen(dpssIN)如果ev是个空值，则按照公式生成对应的ev值if(all(is.null(ev))) {ev <- dpssToEigenvalues(dw, nw) }dw <- dw*sqrt(deltaT) }频率点是nfft点数的一半，加上偏移的零值矩阵nFreqs <- nFFT %/% 2 + as.numeric(returnZeroFreq)偏移点值如果是1，则反馈0offSet <- if(returnZeroFreq) 0 else 1 注意频率轴是使用默认值设定的单位值。# Note that the frequency axis is set by default to unit-less默认是0.5 hz，其他情况则是按照既定参数设置的。# scaling as 0 through 0.5 cycles/period. The user parameter# dtUnits modifies this scaling in the plot.mtm function.尺度频率是1/采样持续时间scaleFreq <- 1 / as.double(nFFT * deltaT)初始化k个中心能量swz <- NULL ## Percival and Walden H0初始化特征向量能量ssqswz <- NULLswz <- apply(dw, 2, sum)if(k >= 2) {swz[seq(2,k,2)] <- 08041}ssqswz <- as.numeric(t(swz)%*%swz)加窗数据初始化taperedData <- dw * 1需要补领的点=nfft-n点nPadLen <- nFFT - n 如果时间序列非复if(!is.complex(timeSeries)) {补零加窗数据为一个nfft*k的矩阵paddedTaperedData <- rbind(taperedData, matrix(0, nPadLen, k))} else {或者复值矩阵补零paddedTaperedData <- rbind(taperedData, matrix(complex(0,0), nPadLen, k)) }cft等于对加窗数据，使用快速傅里叶变换cft <- mvfft(paddedTaperedData)如果不是复值数据if(!is.complex(timeSeries)) {cft等于cft（1：nfreqs）的数据cft <- cft[(1+offSet):(nFreqs+offSet),]} else {复值数据双边展示cft <- rbind(cft[(nFreqs+offSet+1):nFFT,],cft[(1+offSet):(nFreqs+offSet),])}谱能量等于cft的平方sa <- abs(cft)^2如果时间序列不是复值数值if(!is.complex(timeSeries)) {结果频率点构造resultFreqs <- ((0+offSet):(nFreqs+offSet-1))*scaleFreq } else {结果频率点构造resultFreqs <- (-(nFreqs-1):(nFreqs-2))*scaleFreq}初始化自适应参数adaptive <-  NULL初始化剪刀参数jk <- NULL初始化功率谱频率密度PWdofs <- NULL如果不使用剪刀法if(!jackknife) {如果 是实数序列if(!is.complex(timeSeries)) {使用mw2wta法计算自适应的结果adaptive <- .mw2wta(sa, nFreqs, k, sigma2, deltaT, ev)} else {如果是复数序列adaptive <- .mw2wta(sa, nFFT, k, sigma2, deltaT, ev) }其他} else {如果概率密度不符合要求，则停止程序stopifnot(jkCIProb < 1, jkCIProb > .5)如果是实数序列，同时采用自适应方法if(!is.complex(timeSeries) & adaptiveWeighting) {自适应计算adaptive <- .mw2jkw(sa, nFreqs, k, sigma2, deltaT, ev)} else {adaptive <- .mw2jkw(sa, nFFT, k, sigma2, deltaT, ev)}计算对应的比例尺scl <- exp(qt(jkCIProb,adaptive$dofs)*sqrt(adaptive$varjk))上届=自适应s*sclupperCI <- adaptive$s*scllowerCI <- adaptive$s/scl下界，最小值最小的哪一个minVal = min(lowerCI)上界，最大值最大的哪一个maxVal = max(upperCI)jk剪刀值，因为我计算的时候没用，具体参数也不是很明白jk <- list(varjk=adaptive$varjk,bcjk=adaptive$bcjk,sjk=adaptive$sjk,upperCI=upperCI,lowerCI=lowerCI,maxVal=maxVal,minVal=minVal)} 自适应程序有可能没有很好的停止## Short term solution to address bug noted by Lenin Castillo noting that adaptive weights are not properly turned off (Karim 2017).     初始化特征谱，自由度resSpec <- NULLdofVal <- NULL如果，不是自适应的方法if(!adaptiveWeighting) {获得总能量的平均值，自由度2kresSpec <- apply(sa, 1, mean)dofVal <- 2*k} else {谱来自于自适应resSpec <- adaptive$s自由度来自于具体参数dofVal <- adaptive$dofs}                    f检验，空值ftestRes <- NULL
如果要进行f检验if(Ftest) {如果特征频率为空if(is.null(swz)) {则swz可以应用swz <- apply(dw, 2, sum)}使用。hf4mpl计算f检验的结果ftestRes <- .HF4mp1(cft,swz,k,ssqswz)}初始化特征权重，加权因子eigenCoef1 <- NULLwtCoef1 <- NULL如果返回频率间隔if(returnInternals) {特征因子是傅里叶变换eigenCoef1 <- cft如果使用自适应的方法if(adaptiveWeighting) {则加权因子是自适应因子的平方根wtCoef1 <- sqrt(adaptive$wt)} }自适应参数列表auxiliary <- list(dpss=dpssIN,eigenCoefs=eigenCoef1,eigenCoefWt=wtCoef1,nfreqs=nFreqs,nFFT=nFFT,jk=jk,Ftest=ftestRes$Ftest,cmv=ftestRes$cmv,dofs=dofVal,nw=nw,k=k,deltaT=deltaT,dtUnits=dtUnits,taper="dpss")##   Thomson, D.J. Spectrum Estimation and Harmonic Analysis,##   Proceedings of the IEEE, 1982.## note that the weights are squared, they are |d_k(f)^2 from equation## (5.4)## These weights correspond to Thomoson's 1982 Fortran code.## dof fix for one taper, only value.如果k=1自由度为2if(k==1) {auxiliary$dofs <- 2}spec谱的结果spec.out <- list(origin.n=n,method="Multitaper Spectral Estimate",pad= nFFT - n,spec=resSpec,freq=resultFreqs,series=series,adaptive=adaptiveWeighting, mtm=auxiliary)class(spec.out) <- c("mtm", "spec")if(Ftest) {class(spec.out) <- c("mtm", "Ftest", "spec")}返回谱计算的结果return(spec.out)
}

子函数：dpssHelper

获得对应的参数

.dpssV <- function(obj) obj$v.dpssEigen <- function(obj) obj$eigen
使用相应的计算手段。
.is.dpss <- function(obj) {return( sum ( "dpss"==class(obj) ) >= 1 )
}

R学习_multitaper包解析2：子函数spec.mtm.dpss，dpssHelper相关推荐

R学习_multitaper包解析1：子函数centre，dpss
前言之前讲了MTM(多锥形窗谱估计)的相关原理,现在来分析一下它的R语言的实现,这个实现是提出人的学生写的,和matlab的实现进行对照分析,加深理解,提高大家对这门技术的掌握程度. 想要复习的可以 ...
R语言dplyr包排序及序号函数实战（row_number、ntile、min_rank、dense_rank、percent_rank、cume_dist）
R语言dplyr包排序及序号函数实战(row_number.ntile.min_rank.dense_rank.percent_rank.cume_dist) 目录 R语言dplyr包排序及序号函数实 ...
R语言mgcv包中的gam函数拟合广义加性模型（Generalized Additive Model)GAM（对非线性变量进行样条处理、计算RMSE、R方、调整R方、可视化模型预测值与真实值的曲线）
R语言mgcv包中的gam函数拟合广义加性模型(Generalized Additive Model)GAM(对非线性变量进行样条处理.计算RMSE.R方.调整R方.可视化模型预测值与真实值的曲线) ...
R语言 car包recode（）函数被dplyr包里的同名recode（）函数覆盖导致出错
R语言 Error相关处理 recode函数报错 recode函数报错用r语言做分位数回归使用recode()函数对数据进行重新编码,本来应该调用car包里面的recode函数, 即 recode( ...
R语言dplyr包：高效数据处理函数（filter、group_by、mutate、summarise）
R语言dplyr包的数据处理.分析函数在日常数据处理过程中难免会遇到些难处理的,选取更适合的函数分割.筛选.合并等实在是大快人心! 利用dplyr包中的函数更高效的数据清洗.数据分析,及为后续数据建 ...
R语言psych包的corr.test函数计算相关性并给出所有相关性的显著性（Correlation matrix and tests of significance via corr.test()）
R语言使用psych包的corr.test函数计算所有变量组合对的相关性并给出所有相关性的显著性(Correlation matrix and tests of significance via co ...
R语言leaps包中的regsubsets函数实现全集子集回归(all subsets regression)、使用调整R方和Mallows Cp统计量筛选最优模型、并可视化不同组合参数下的模型指标
R语言使用leaps包中的regsubsets函数实现全集子集回归(All Subsets Regression,ASR).使用调整R方和Mallows Cp统计量筛选最佳模型.并可视化不同组合参数下 ...
R语言dplyr包超完整版函数指南
R语言dplyr包的使用一.常用函数功能速查二.常用函数详解 iris数据集 1. 取子集 filter/distinct slice select 2. 窗口函数 3. 连接合并 left_jo ...
R语言dplyr包：高效数据处理函数arrange、sample_n、n_distinct、select、compute等
今天是个特别的日子,小编在这里祝大家情人节快乐!本篇文章继续之前文章提到的关于dplyr包数据处理的函数.错了,小编是准备那天发的,忘发了 R语言在数据整理.分析上面的方法是很多的,并且通俗易懂,相信 ...

R学习_multitaper包解析2：子函数spec.mtm.dpss，dpssHelper

前言

目录

子函数：spec.mtm.dpss

子函数：dpssHelper

R学习_multitaper包解析2：子函数spec.mtm.dpss，dpssHelper相关推荐

最新文章

热门文章