前言

之前講了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
• 子函數：dpssHelper

子函數：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的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。