[Dplr-commits] r671 - in branches/redfit: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Aug 29 19:07:57 CEST 2013
Author: mvkorpel
Date: 2013-08-29 19:07:57 +0200 (Thu, 29 Aug 2013)
New Revision: 671
Modified:
branches/redfit/R/redfit.R
branches/redfit/src/redfit.c
Log:
* Variable 'wz' (difference between successive angular frequencies) is
not used anymore. Instead, angular frequencies are computed as '2 *
pi * freq', where 'freq' is a frequency from the frequency vector
obtained with 'seq(from = 0, to = fnyq, length.out = nfreq)': It's
best to let R create a vector of evenly spaced frequencies.
* The formula for 'gredth' was also simplified in the same spirit.
* Instead of using 'freq[2]' with the comment 'NB: freq[2] = df', just
use 'df'.
Change visible to the user: the return value of redfit() no longer
contains $params$wz.
Modified: branches/redfit/R/redfit.R
===================================================================
--- branches/redfit/R/redfit.R 2013-08-29 16:31:57 UTC (rev 670)
+++ branches/redfit/R/redfit.R 2013-08-29 17:07:57 UTC (rev 671)
@@ -211,20 +211,19 @@
fnyq <- params[["fnyq"]]
nfreq <- params[["nfreq"]]
df <- params[["df"]]
- wz <- params[["wz"]]
ofac <- params[["ofac"]]
segskip <- params[["segskip"]]
- ia <- redfitInitArrays(t, x, params)
+ freq <- seq(from = 0, to = fnyq, length.out = nfreq)
+ ia <- redfitInitArrays(t, x, freq, params)
## determine autospectrum of input data
dn50 <- as.numeric(n50)
cbindfun <- match.fun("cbind")
lmfitfun <- tryCatch(match.fun(".lm.fit"),
error = function(...) match.fun("lm.fit"))
gxx <- .Call(dplR.spectr, t, x, np, ia[[1]], ia[[2]], ia[[3]], ia[[4]],
- nseg, nfreq, avgdt, wz, dn50, segskip, cbindfun, lmfitfun)
- freq <- seq(from = 0, by = 1, length.out = nfreq) * df
+ nseg, nfreq, avgdt, freq, dn50, segskip, cbindfun, lmfitfun)
## estimate data variance from autospectrum
- varx <- freq[2] * sum(gxx) # NB: freq[2] = df
+ varx <- df * sum(gxx)
## dplR: estimate lag-1 autocorrelation coefficient unless prescribed
if (is.null(rhopre) || rhopre < 0) {
rho <- redfitGetrho(t, x, np, n50, nseg, avgdt, segskip)
@@ -246,9 +245,9 @@
grr[, i] <-
.Call(dplR.spectr, t, .Call(dplR.makear1, difft, np, tau), np,
ia[[1]], ia[[2]], ia[[3]], ia[[4]], nseg, nfreq, avgdt,
- wz, dn50, segskip, cbindfun, lmfitfun)
+ freq, dn50, segskip, cbindfun, lmfitfun)
## scale and sum red-noise spectra
- varr1 <- freq[2] * sum(grr[, i]) # NB: freq[2] = df
+ varr1 <- df * sum(grr[, i])
grr[, i] <- varx / varr1 * grr[, i]
}
grrsum <- rowSums(grr)
@@ -261,9 +260,9 @@
## setup AR(1) time series and estimate its spectrum
grr <- .Call(dplR.spectr, t, .Call(dplR.makear1, difft, np, tau),
np, ia[[1]], ia[[2]], ia[[3]], ia[[4]], nseg, nfreq,
- avgdt, wz, dn50, segskip, cbindfun, lmfitfun)
+ avgdt, freq, dn50, segskip, cbindfun, lmfitfun)
## scale and sum red-noise spectra
- varr1 <- freq[2] * sum(grr) # NB: freq[2] = df
+ varr1 <- df * sum(grr)
grr <- varx / varr1 * grr
grrsum <- grrsum + grr
}
@@ -272,13 +271,14 @@
## determine average red-noise spectrum; scale average again to
## make sure that roundoff errors do not affect the scaling
grravg <- grrsum / nsim
- varr2 <- freq[2] * sum(grravg)
+ varr2 <- df * sum(grravg)
grravg <- varx / varr2 * grravg
rhosq <- rho * rho
## set theoretical spectrum (e.g., Mann and Lees, 1996, Eq. 4)
## make area equal to that of the input time series
- gredth <- (1 - rhosq) / (1 + rhosq - 2 * rho * cos(pi / fnyq * freq))
- varr3 <- freq[2] * sum(gredth)
+ gredth <- (1 - rhosq) /
+ (1 + rhosq - 2 * rho * cos(seq(from = 0, to = pi, length.out = nfreq)))
+ varr3 <- df * sum(gredth)
gredth <- varx / varr3 * gredth
## determine correction factor
corr <- grravg / gredth
@@ -565,21 +565,20 @@
invisible(x)
}
-redfitInitArrays <- function(t, x, params) {
+redfitInitArrays <- function(t, x, freq, params) {
np <- params[["np"]]
nseg <- params[["nseg"]]
- nfreq <- params[["nfreq"]]
+ nfreqM1 <- length(freq) - 1
n50 <- params[["n50"]]
iwin <- params[["iwin"]]
- wz <- params[["wz"]]
segskip <- params[["segskip"]]
ww <- matrix(NA_real_, nseg, n50)
- tsin <- array(NA_real_, c(nseg, nfreq - 1, n50))
- tcos <- array(NA_real_, c(nseg, nfreq - 1, n50))
- wtau <- matrix(NA_real_, nfreq - 1, n50)
+ tsin <- array(NA_real_, c(nseg, nfreqM1, n50))
+ tcos <- array(NA_real_, c(nseg, nfreqM1, n50))
+ wtau <- matrix(NA_real_, nfreqM1, n50)
for (i in as.numeric(seq_len(n50))) {
twk <- t[.Call(dplR.seg50, i, nseg, segskip, np)]
- tr <- redfitTrig(twk, nseg, wz, nfreq)
+ tr <- redfitTrig(twk, freq)
ww[, i] <- redfitWinwgt(twk, iwin)
wtau[, i] <- tr[[3]]
tsin[, , i] <- tr[[1]]
@@ -616,7 +615,6 @@
fnyq <- hifac / (2 * avgdt) # average Nyquist freq.
nfreq <- floor(hifac * ofac * nseg / 2 + 1) # f[1] == f0; f[nfreq] == fNyq
df <- fnyq / (nfreq - 1) # freq. spacing
- wz <- 2 * pi * fnyq / (nfreq - 1) # omega == 2*pi*f
if (verbose) {
cat(" N = ", np, "\n", sep="")
cat(" t[1] = ", t[1], "\n", sep="")
@@ -627,7 +625,7 @@
}
## dplR: ditched nout (nout == nfreq)
res <- list(np = np, nseg = nseg, nfreq = nfreq, avgdt = avgdt, df = df,
- wz = wz, fnyq = fnyq, n50 = n50, ofac = ofac, hifac = hifac,
+ fnyq = fnyq, n50 = n50, ofac = ofac, hifac = hifac,
segskip = segskip)
args <- list(...)
argnames <- names(args)
@@ -644,19 +642,21 @@
res
}
-redfitTrig <- function(tsamp, nn, wz, nfreq) {
+redfitTrig <- function(tsamp, freq) {
tol1 <- 1.0e-4
- nfreqM1 <- nfreq - 1
+ nfreqM1 <- length(freq) - 1
+ nn <- length(tsamp)
tcos <- matrix(NA_real_, nn, nfreqM1)
tsin <- matrix(NA_real_, nn, nfreqM1)
wtau <- numeric(nfreqM1)
+ wfac <- 2 * pi # omega == 2*pi*f
## start frequency loop
## dplR: In the original Fortran code, the variables ww (not used
## in this function), wtau, tsin and tcos have unused elements
## (one extra frequency). The unused elements have now been
## dropped.
for (k in seq_len(nfreqM1)) {
- wrun <- k * wz
+ wrun <- wfac * freq[k + 1]
## calc. tau
arg2 <- wrun * tsamp
arg1 <- arg2 + arg2
Modified: branches/redfit/src/redfit.c
===================================================================
--- branches/redfit/src/redfit.c 2013-08-29 16:31:57 UTC (rev 670)
+++ branches/redfit/src/redfit.c 2013-08-29 17:07:57 UTC (rev 671)
@@ -26,10 +26,10 @@
SEXP seg50(SEXP k, SEXP nseg, SEXP segskip, SEXP np);
void rmtrend(SEXP x, SEXP y, SEXP lmfit);
SEXP spectr(SEXP t, SEXP x, SEXP np, SEXP ww, SEXP tsin, SEXP tcos, SEXP wtau,
- SEXP nseg, SEXP nfreq, SEXP avgdt, SEXP wz, SEXP n50,
+ SEXP nseg, SEXP nfreq, SEXP avgdt, SEXP freq, SEXP n50,
SEXP segskip, SEXP cbind, SEXP lmfit);
void ftfix(const double *xx, const double *tsamp, const size_t nxx,
- const double wz, const size_t nfreq, const long double si,
+ const double *freq, const size_t nfreq, const long double si,
const size_t lfreq, const double tzero, const double *tcos,
const double *tsin, const double *wtau, const long double sumbysqrt,
double *ftrx, double *ftix);
@@ -142,16 +142,16 @@
/* dplR: Returns the spectrum of x(t), a vector of length nfreq.
*/
SEXP spectr(SEXP t, SEXP x, SEXP np, SEXP ww, SEXP tsin, SEXP tcos, SEXP wtau,
- SEXP nseg, SEXP nfreq, SEXP avgdt, SEXP wz, SEXP n50,
+ SEXP nseg, SEXP nfreq, SEXP avgdt, SEXP freq, SEXP n50,
SEXP segskip, SEXP cbind, SEXP lmfit) {
SEXP gxx, twk, xwk, ftrx, ftix, tmp, cbindcall;
- double wz_val, dnseg, sqrt_nseg, segskip_val, scal, np_val;
+ double dnseg, sqrt_nseg, segskip_val, scal, np_val;
long double sumx;
size_t i, j, nseg_val, nfreq_val, n50_val, segstart;
size_t sincos_skip, wtau_skip;
size_t wwidx = 0;
double *t_data, *x_data, *ww_data, *tsin_data, *tcos_data, *wtau_data;
- double *gxx_data, *twk_data, *xwk_data, *ftrx_data, *ftix_data;
+ double *gxx_data, *twk_data, *xwk_data, *ftrx_data, *ftix_data, *freq_data;
const long double si = 1.0;
const double tzero = 0.0;
const size_t lfreq = 0;
@@ -161,7 +161,6 @@
nseg_val = (size_t) dnseg;
nfreq_val = (size_t) *REAL(nfreq);
np_val = *REAL(np);
- wz_val = *REAL(wz);
n50_val = (size_t) *REAL(n50);
segskip_val = *REAL(segskip);
t_data = REAL(t);
@@ -170,6 +169,7 @@
tsin_data = REAL(tsin);
tcos_data = REAL(tcos);
wtau_data = REAL(wtau);
+ freq_data = REAL(freq);
PROTECT(gxx = allocVector(REALSXP, nfreq_val));
PROTECT_WITH_INDEX(twk = allocVector(REALSXP, nseg_val), &pidx);
@@ -220,7 +220,7 @@
sumx += xwk_data[j];
}
/* Lomb–Scargle Fourier transform */
- ftfix(xwk_data, twk_data, nseg_val, wz_val, nfreq_val, si,
+ ftfix(xwk_data, twk_data, nseg_val, freq_data, nfreq_val, si,
lfreq, tzero, tcos_data, tsin_data, wtau_data,
sumx / (long double) sqrt_nseg, ftrx_data, ftix_data);
/* dplR: adjust tsin, tcos, wtau for next segment */
@@ -263,7 +263,7 @@
* index corresponding to nxx runs faster.
*/
void ftfix(const double *xx, const double *tsamp, const size_t nxx,
- const double wz, const size_t nfreq, const long double si,
+ const double *freq, const size_t nfreq, const long double si,
const size_t lfreq, const double tzero, const double *tcos,
const double *tsin, const double *wtau, const long double sumbysqrt,
double *ftrx, double *ftix) {
@@ -284,7 +284,7 @@
ftix[0] = 0.0;
/* start frequency loop */
for (ii = 1; ii < nfreq; ii++) {
- wrun = (double) ii * wz;
+ wrun = M_2PI * freq[ii]; /* omega = 2 * pi * freq */
wtnew = wtau[ii - 1];
/* summations over the sample */
cross = 0.0;
More information about the Dplr-commits
mailing list