前言
之前讲了MTM(多锥形窗谱估计)的相关原理,现在来分析一下它的R语言的实现,这个实现是提出人的学生写的,和matlab的实现进行对照分析,加深理解,提高大家对这门技术的掌握程度。
想要复习的可以参考一下之前的文件:
mtm原理:现代谱估计:多窗口谱
想要复习一下如何实现的可以参考:
MTM:matlab实现1MTM:matlab实现1
MTM:matlab实现2参数解析MTM参数解析
MTM:matlab实现3谱功率计算MTM谱功率计算
MTM:matlab实现4主函数解析MTM 主函数
目录
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是时间序列,则获得它的采样间隔,否则,默认为1
if(is.null(deltaT)) {
if(is.ts(x)) {
deltaT <- deltat(ts)
} else {
warning("deltaT not specified; using deltaT=1.")
deltaT <- 1.0
}
}
获得时间序列x的长度n
n <- 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$eigen
swz是特征向量的列值和
swz <- apply(dw, 2, sum)
## zero swz where theoretically zero; odd tapers
如果k》2则swz偶数位的值是0
if(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 properly
if(!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
nw
dpssToEigenvalues <- 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))/npot
j <- 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大于1
if(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])
}