I am trying to create a set of wavelets in frequency space--namely Cauchy
wavelets for an intensity analysis (von Tscharner, 2000). The wavelets are
defined by the following formula:
[(f/cf)^(cf*scale)]*[e^((-f/cf)+1)^(cf*scale)]
where *f *is frequency of length *n*, *cf* is center frequency (defined
below) and is an array of *j *columns and *n* row, and scale is a constant.
cf = (1/scale)*(j + q)^r
where *j *is the number of center frequencies (i.e., columns), and *q* and *
r* are constants.
When I get to Q, the output returns numbers only up to cf. Beyond cf, Q and
subsequently W return "NA." I understand this is because when fq >/= cf U
becomes zero^B[j] or a negative^B[j]. From here, Q should equal 1/exp(U) to
compute W. I know that I need to create a loop here, but as a newbie to R,
I don't know how to write the correct loop that will work for me. Any
suggestions??
Here is the code I have been using (less the numerous if...else attempts):
> J <- 12
> r_ <- 1.959
> q_ <- 1.45
> scale_ <- 0.3
> N <- 500
> fq <- seq(0, N, 1)
> center_frequencies <- function(J = 12, r_ = 1.959, q_ = 1.45, scale_ =
0.3){
+ j <- seq(0,J-1,1)
+ fc <- (q_ + j)^r_/scale_
+ }
> fc <- center_frequencies(12,r_,q_,scale_)
> cf <- t(fc)
> lambda <- function(cf, J = 12, scale_ = 0.3){
+ B <- cf*scale_
+ }
> B <- lambda(cf, 12, 0.3)
> dummy <- fq/cf[1]
> Z <- dummy**B[1] # Note here I tried using '^' to raise dummy to the power
B[j], but it didn't work. I tried '**' which serves the same purpose in
LabVIEW, and it worked. There is no references to this in any of the R
documentation that I have read. Should '^' work here???
> U <- (-dummy+1)**B[1]
> Q <- exp(U)
> Q <- ifelse(Q>30,30,Q) # Not sure why I use this.
> W <- Z*Q
--
*W. Jeffrey Armstrong, Ph.D.
*Assistant Professor
Exercise Science
*Managing Editor
Clinical Kinesiology*
Official Journal of the American Kinesiotherapy Association
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
