前言

之前講了MTM（多錐形窗譜估計）的相關原理，現在來分析一下它的R語言的實現，這個實現是提出人的學生寫的，和matlab的實現進行對照分析，加深理解，提高大家對這門技術的掌握程度。

想要復習的可以參考一下之前的文件：

mtm原理：現代譜估計：多窗口譜
想要復習一下如何實現的可以參考：
MTM：matlab實現1MTM：matlab實現1
MTM：matlab實現2參數解析MTM參數解析
MTM：matlab實現3譜功率計算MTM譜功率計算
MTM：matlab實現4主函數解析MTM 主函數

目錄

• 前言
• 目錄
• centre函數
• DPSS函數

centre函數

該函數負責對輸入的時間序列進行處理，使其滿足零均值條件

######################################################################### ## ## centre ## ## Takes a time series and converts to zero-mean using one of three ## methods: Slepian projection, arithmetic mean, or trimmed mean. ## ######################################################################### r里面的賦值采用<-格式 輸入參數 時間序列x nw生成幾階dpss序列 k使用的窗口數 deltat 采樣間隔 trim截尾平均時丟掉的點數 centre <- function(x, nw=NULL, k=NULL, deltaT=NULL, trim=0) { 判斷x是否有非數值na.fail(x)初始化結果res <- NULL如果沒有輸入dpss參數的話，則默認為截尾平均if(is.null(nw) && is.null(k) ) {res <- x - mean(x, trim=trim)如果有dpss參數，則使用slepian投影，截尾數舍棄} else {if(trim != 0) {warning(paste("Ignoring trim =", trim))}如果采用dpss方法，且相關參數不符合條件，則停止函數stopifnot(nw >= 0.5, k >= 1, nw <= 500, k <= 1.5+2*nw)如果dpss階數和長度之比不合適，停止運算if (nw/length(x) > 0.5) { stop("half-bandwidth parameter (w) is greater than 1/2")}如果，deltat沒有值且x是時間序列，則獲得它的采樣間隔，否則，默認為1if(is.null(deltaT)) {if(is.ts(x)) {deltaT <- deltat(ts)} else {warning("deltaT not specified; using deltaT=1.")deltaT <- 1.0}}獲得時間序列x的長度nn <- length(x)獲得生成的DPSS序列結果，輸入參數為n，k，nw，返回特征值dpssRes <- dpss(n, k=k, nw=nw,returnEigenvalues=TRUE)dw值等于返回結果中的特征值*采樣間隔的單位根dw <- dpssRes$v*sqrt(deltaT)ev值等于返回結果中的特征加窗矩陣，每一列對應第k個。ev <- dpssRes$eigenswz是特征向量的列值和swz <- apply(dw, 2, sum)## zero swz where theoretically zero; odd tapers如果k》2則swz偶數位的值是0if(k >=2) {swz[seq(2,k,2)] <- 0.0}ssqswz是幅頻率功率之和。ssqswz <- sum(swz^2)如果不是復值x，對x采取meave函數進行處理if(!is.complex(x)) {res <- .mweave(x, dw, swz,n, k, ssqswz, deltaT)res <- x - res$cntr} else {res.r <- .mweave(Re(x), dw, swz,n, k, ssqswz, deltaT)res.i <- .mweave(Im(x), dw, swz,n, k, ssqswz, deltaT)res <- x - complex(real=res.r$cntr, imaginary=res.i$cntr)}}return(res) }

DPSS函數

負責生成給定的k個正交的離散扁球序列，即用來加權的窗
使用對角法生成。
參照論文和譜分析教材那本書，會講如何得出給定的值

########################################################## ## ## dpss ## ## Generates k orthogonal discrete prolate spheroidal ## sequences (dpss) using the tridiagonal method. See ## Slepian (1978) page 1379 and Percival and Walden ## chapter 8.4 ## ########################################################## 函數dpss 輸入 n值序列個數 k值加窗個數 nw值窗口給定值 ReturnEigenvalues=true真 返回特征值dpss <- function(n, k, nw, returnEigenvalues=TRUE) { 如果不滿足條件停止計算stopifnot(n >= 1, nw/n >0, nw/n < 0.5, k >= 1)將k變成一個整數。## if k is passed in as floating point, the cast to ## as.integer() in the Fortran call does not quite work properlyif(!is.integer(k)) {k<-as.integer(floor(k));} 使用專門的包計算特征值##eigen is of length for use by lapack functoins.## this will use lapack functions in place of the## eispack functions referenced in Percival and Waldern輸入相應參數，得到相應的結果out <- .Fortran("fdpss", as.integer(n), as.integer(k),as.double(nw), v=double(n*k), eigen=double(k),PACKAGE='multitaper')輸出v值矩陣 out$v <- matrix(data=out$v, nrow=n, ncol=k, byrow=FALSE)如果返回特征值為真if(returnEigenvalues) {輸出的特征值利用特定子函數生成out$eigen <- dpssToEigenvalues(out$v, nw)} else {特征值由對稱公式返回## eigen values returned from the tridiagonal formulation## Slepian eqn #13 (1978)out$eigen <- out$eigen}結果是輸出的v和特征值res <- list(v=out$v,eigen=out$eigen)結果類型 dpss? class(res) <- "dpss"return(res) }########################################################## ## ## dpssToEigenvalues ## 給定一個離散扁球窗，找到對應生成dpss序列的特征值 ## Given a set of dpss tapers, find the eigenvalues corresponding ## to the generated dpss's ## ## See Percival and Walden (1993) exercise 8.4, and ## associated LISP code. ## ##########################################################dpss序列到特征值矩陣 輸入 v nwdpssToEigenvalues <- function(v, nw) { 將v矩陣化v <- as.matrix(v)獲得v的第一維序列長度n <- length(v[,1])獲得v的第二維，k值k <- length(v[1,])獲得窗口百分比w <- nw/n找到序列長度最接近的2的次方的npot數量npot <- 2**(ceiling(log2(2*n)))## pad生成v值矩陣scratch <- rbind(v, matrix(data=0, nrow=npot-n, ncol=k))## n * acvs矩陣的值是對應特偵序列的功率。scratch <- Re(mvfft(abs(mvfft(scratch))**2))/npotj <- 1:(n-1)生成比率向量vectorOfRatios <- c(2*w, sin(2*pi*w*j)/(pi*j))這兩個矩陣都隨著序號的上升，幅值下降## Note: both vector of ratios and scratch## roughy decrease in amplitude with increasing index,## so sum things up in reverse order特征值為空eigenvalues <- NULL如果k大于1if(k>1) {特征值按列相乘假和eigenvalues <- apply(scratch[n:2,] * vectorOfRatios[n:2], 2, sum)} else {其他情況下不用加eigenvalues <- sum(scratch[n:2,] * vectorOfRatios[n:2])}返回特征向量加上第一列的和。return(2*eigenvalues +vectorOfRatios[1]*scratch[1,1:k]) }

總結

以上是生活随笔為你收集整理的R学习_multitaper包解析1：子函数centre，dpss的全部內容，希望文章能夠幫你解決所遇到的問題。

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