Newer
Older
Q <- function(wT) {sin(wT)/(wT)*(pi^2)/(pi^2-(wT)^2)}
Qprime <- function(y) {ifelse(y==0, 0, pi^2/(pi^2-y^2)/y*(cos(y)+(sin(y)/y)*(3*y^2-pi^2)/(pi^2-y^2)))}
Qsecond0 <- 2/pi^2 - 1./3.
mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, maxfreq=NULL){
# nu can be provided, in which case the frequencies
# are considered to be given
if (!is.null(nu))
{
# we will need two versions of nfreq, one for the cosine and one for the cosine
# if if there is a zero frequency we only need it once
nutemp <- unlist ( lapply( nu, function(a) if (a==0) return (a) else return ( c(a,-a))) )
phase <- unlist ( lapply( nu, function(a) if (a==0) return (0) else return ( c(0, pi/2))) )
nu <- nutemp
if (length(nu) < 2*nfreq) {
# this will typically occur if one zero frequency was provided
nu[2*nfreq] <- NA
phase[2*nfreq] <- NA
}
} else
{
nu <- rep(NA,2*nfreq)
phase <- rep(NA,2*nfreq)
}
N <- length(xdata)
T <- N / 2
hann <- function(N) {return (1-cos(2*pi*seq(0,(N-1))/(N)))};
hN <- hann(N)
t <- seq(N)-1
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
hprod <- function(x,y) sum(x * hN * y)/N
# B_i = sum[m<= i](A_im) f_m the orthogonal base signals
# x_m what remains of the signal at time step m
S <- rep(0,2*nfreq)
Prod <- rep(0,2*nfreq)
amp <- rep(NA,2*nfreq)
A <- matrix(0,2*nfreq,2*nfreq)
Q <- matrix(0,2*nfreq,2*nfreq)
f <- list()
x <- list()
B <- list()
# S[m] = hprod (f[m], B[m])
freqs = 2.*pi*seq(0,(N-1))/N
x[[1]] = xdata
for (m in seq(2*nfreq)){
hx <- hN*x[[m]]
# there is a little tweak here:
# if we have reached 2*nfreq and not nu[m] is na, then
# it means that we have encountered a constant somewhere, and that last step should be skipped.
if ( ! ( m == 2*nfreq && is.na(nu[m]))) {
# ok we are on business
if (is.na(nu[m])){
# is the "m" frequency already provided ?
# either there were provided from the start and developped in the first lines of this routine
# or either they were not provided, but we are in this case where it was set in "m-1" because
# we identified a non-null frequency
# in both configurations, we look at a frequency with the fourier transform
fbase <- freqs[which.max(power(hx))]
brackets <- c(fbase-pi/N, fbase+pi/N);
# and then further corrects the brackets if minfreq and maxfreq where provided
brackets[1] <- max(minfreq, brackets[1])
brackets[2] <- max(maxfreq, brackets[1])
thx <- t(hx)
# after profiling, the fastest seems the first option below
# this is the weak link
# coding the function in c might be an option
fmax = cmna::goldsectmax(function(f)
{
ft <- f*t
(thx %*% cos(f*t))^2 + (thx %*% sin(f*t))^2},
brackets[1], brackets[2], tol=1.e-10, m=9999)
# fmax = cmna::goldsectmax(function(f) {
# ft <- f*t
# (sum(hx * cos(ft)))^2 + (sum(hx %*% sin(ft)))^2},
# brackets[1], brackets[2], tol=1.e-10, m=9999)
# fmax = cmna::goldsectmax(function(f) {
# ft <- f*t
# Mod(sum(hx * exp(1i*ft)))},
# brackets[1], brackets[2], tol=1.e-10, m=9999)
if (fmax > freqs[2]/2){
# if we really identified a frequency
# we are in this case where the frequency was not a priori provided
# we se set this phase to zero, and the next one to pi/2, and also set the frequencies accordinly
phase[m] <- 0.
nu[m] <- fmax
phase[m+1] <- pi/2.
nu[m+1] <- -fmax
} else {
# again we are still in thes case where a freqency was not a priori provided
# but the more particular case where the frequency is (with tolerance) considered as zero
# f[[m]] will effectively be constant.
phase[m] <- 0.
nu[m] <- 0.
}}
f[[m]] <- cos(nu[m]*t + phase[m])
if (fast) # we use analytical forms of the products
{
for (i in seq(m))
{
num <- (nu[m] - nu[i])*T
nup <- (nu[m] + nu[i])*T
phim <- (phase[m] - phase[i])
phip <- (phase[m] + phase[i])
Qm <- ifelse(num == 0, 1 , Q(num))
Qp <- ifelse(nup == 0, 1 , Q(nup))
cosm <- cos(phim) * cos(num) * Qm
sinm <- sin(phim) * sin(num) * Qm
cosp <- cos(phip) * cos(nup) * Qp
sinp <- sin(phip) * sin(nup) * Qp
Q[m,i] = 0.5 * (cosm + cosp - sinm - sinp)
}
} else {
for (i in seq(m))
{Q[m,i] = hprod(f[[m]],f[[i]]) # so remember the convetion
# these are symmetric matrices
# and we only populate the j <= m
}
}
#
# before normalisation
# B[[m]] = sum_j=1,m A[m,j]*f[j] et
# B[[m]] = ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i]
# B[[m]] = ( f[m] - sum_(1,m-1)
# (f[m] * sum(j=1,i A_j,i f[i]) ) *
# ( sum (A_j=1, i) A_j,i f[i] )
#/ ||B[[m]]||
# one can then verify that:
# B[[m]] %*% B[j]] =
# ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i] ) %*% B[j] =
# ( f[m] - (f[m] %*% B[j]) * B[j] ) %*% B[j] =
# ( f[m] %*% B[j] - f[m] %*% B[j]) = 0
# before normalisation
A[m,] = 0
A[m, m] = 1.
if (m>1){
# fmbi = f[m] %*% B[i]
fmbi <- rep(0,(m-1))
for (j in seq(m-1)) for (i in seq(j)) fmbi[j] = fmbi[j] + A[j,i]*Q[m,i]
for (j in seq(m-1)) for (i in seq(j,(m-1))) A[m,j] = A[m,j] - fmbi[i]*A[i,j]
}
# so, at this stage, we have B[[m]] = sum [A[m,j]] * f[j]
# B[[2]] = (A[2,1]*f[1] + A[2,2]*f[2])
# B[[3]] = (A[3,1]*f[1] + A[3,2]*f[2] + A[3,3] * f[3] )
norm = 0
if (m > 1){
for (i in seq(m)) {
norm = norm + (A[m,i]*A[m,i]*Q[i,i])
if (i>1) for (j in seq(i-1)) { norm = norm + 2*A[m,i]*A[m,j]*Q[i,j]
}}
} else {
norm <- A[m,m]*A[m,m]*Q[m,m]
}
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
Prod[m] = hprod(x[[1]],f[[m]])
# les prods sont corrects
S[m] = 0.
for (j in seq(m)) S[m] = S[m] + Prod[j]*A[m,j]
# for (j in seq(m)) S[m] = S[m] + A[m,j] * hprod(x[[1]], f[[j]])
# les Sm sont les projections dans la nouvelle base
# not necessary, for verification only
# computes the B and verify they are orthonormal
B[[m]]=0
for (j in seq(m)) B[[m]] = B[[m]] + A[m,j] * f[[j]]
#print ('S[m] computed two different ways')
#print (m)
#print (S[m])
#print (hprod(x[[1]], B[[m]]))
#print ('end test')
# for (j in seq(m)) print (sprintf("B%i%i : %3g", i, j, hprod(B[[j]], B[[m]])))
# if you are curious, the amplitude of the signal will be
# amp[j] = contribution of the sum( S[m]*B[m]) to the f[[m]]
# amp[j] = sum(m in seq(j)) S(m) * A[m,j]
# for (i in seq(m)) amp[j] = amp[j]+S[m]*A[m,j]
# and the residual signal
# x[[m+1]] = x[[m]] - S[m] * B[[m]]
x[[m+1]] = x[[m]]
# for (j in seq(m)) x[[m+1]] = x[[m+1]] - S[m] * A[m,j]*f[[j]]
x[[m+1]] = x[[m+1]] - S[m] * B[[m]]
# print('residu final')
# print(m)
# print (x[[m+1]][seq(20)])
}
#xtest <- 0
#xtest2 <- 0
#coeftests <- rep(0,m)
#phi1 <- 0.4 ; phi2 <- 0.1234
# coefreals <- c(cos(phi1), -sin(phi1), cos(phi2), -sin(phi2))
# for (j in seq(m)) xtest <- xtest + coefreals[j] * f[[j]]
# for (j in seq(m)) for (i in seq(j)) xtest2 <- xtest2 + S[j] * A[j,i] * f[[i]]
# for (j in seq(m)) for (i in seq(j)) coeftests[i] <- coeftests[i] + S[j] * A[j,i]
# print ((xtest - x[[1]])[seq(20)])
# print ((xtest2 - x[[1]])[seq(20)])
# print ("cofreals")
# print (coefreals)
# print ("coftests")
# print (coeftests)
# donc: j'ai demontre que les deux xtests sont les memes
# coefreals est la bonne solution
# et je sais ussi que B = sum A_mj f_j
# donc je devroais simplement avour les Sj A_mj
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
# at htis point I would like to make a littel check
# make A full
# is A Q t(A) = BtB orthonormal ?
# for (m in seq(2,2*nfreq)) for (j in seq(m-1)) A[j,m]=A[m,j]
# for (m in seq(2,2*nfreq)) for (j in seq(m-1)) Q[j,m]=Q[m,j]
# print (A*t(Q)*t(A))
# print ('END TEST')
#print (A)
#print (A %*% S)
# amplitudes
mmax <- 2*nfreq
if (is.na(nu[mmax])) mmax <- mmax - 1
amp[1:mmax] = 0;
for (m in seq(mmax)) for (j in seq(m)) amp[j] <- amp[j] + S[m] * A[m,j]
amp2m <- amp;
for (m in seq(mmax)){
# if the perivous "m" was already of the same frequency, we can merge amplitudes and set phase
if ((m > 1) && (nu[m-1] == -nu[m])){
phase[m] <- Arg (-1i * amp[m] + amp[m-1])
amp[m] <- sqrt(amp[m-1]^2 + amp[m]^2)
amp[m-1] <- NA
phase[m-1] <- NA
}
}
# at this point we should have exactly nfreq non na values
# we print a message if this is not right, and in that case I suspect some debugging will be needed.
valid_frequencies <- which(!is.na(amp))
nu <- -nu[valid_frequencies]
amp <- amp[valid_frequencies]
phase <- phase[valid_frequencies]
if (length(valid_frequencies) != nfreq) message (sprintf("something goes wrong : %i valid frequencies, and nfreq = %i", valid_frequencies, nfreq))
OUT = data.frame(nu=nu, amp=amp, phase=phase)
return(OUT)
}
#' Frequency Modified Fourier transform for real series
#'
#' R-coded version of the Modified Fourier Transform
#' with frequency correction, adapted to R.
#' much slower than mfft (for complex numbers) as the latter is
#' mainly written in C, but is physically
#' more interpretable if signal is real, because
#' it is designed to have no imaginary part in the residual
#' A C-version should be supplied one day.
#'
#' @importFrom cmna goldsectmax
#' @param xdata The data provided either as a time series (advised), or as a vector.
#' @param minfreq,maxfreq If provided, bracket the frequencies to be probed. Note this are
#' more precisely angular velocities (2\pi / period), expressed in time-inverse units
#' with the time resolution encoded in `xdata` if the latter is a time series.
#' @param correction: 0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data (both are documented in the Sidlichovsky and Nesvorny paper.
#' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata.
#' @param fast (default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines.
#' note: this is not really faster because the bottleneck is actually the goden section search. But more elegant.
#' @return a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases".
#' @author Michel Crucifix
#' @references
#' \insertRef{sidlichovsky97aa}{gtseries}
#' @examples
#'
#' data(harmonic_sample)
#' spectrum <- mfft_real(harmonic_sample$data)
#' print(spectrum)
mfft_real <- function(xdata, nfreq=5, minfreq=NULL, maxfreq=NULL, correction = 1 , fast=TRUE){
N <- length(xdata)
N2 <- N/2.
my_minfreq <- ifelse (is.null(minfreq), 0, minfreq*dt)
my_maxfreq <- ifelse (is.null(maxfreq), pi, maxfreq*dt)
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
startx = stats::start(xdata)[1]
N <- length(xdata)
OUT <- mfft_analyse(xdata, nfreq, fast, NULL, my_minfreq, my_maxfreq)
# correction (methode 2)
if (correction == 2){
xdata_synthetic <- rep(0,N)
t <- seq(N)-1
for (i in seq(nfreq)) xdata_synthetic = xdata_synthetic + OUT$amp[i]*cos(OUT$nu[i]*t + OUT$phase[i])
OUT2 <- mfft_analyse(xdata_synthetic, nfreq, fast, NULL, my_minfreq, my_maxfreq)
OUT$nu = OUT$nu + (OUT$nu - OUT2$nu)
OUT$amp = OUT$amp + (OUT$amp - OUT2$amp)
OUT$phase = OUT$phase + (OUT$phase - OUT2$phase)
} else if (correction == 1){
for (j in seq(nfreq)){
epsilon = OUT$amp[j] * Qprime(-2 * OUT$nu[j] * N2)*cos(2 * OUT$nu[j] * N2 + 2 * OUT$phase[j])
if ((j+1) <= nfreq) {
for (s in seq(j+1, nfreq)) {
epsilon = epsilon + OUT$amp[s] *
(
Qprime( (OUT$nu[s] - OUT$nu[j])*N2)*cos((OUT$nu[j] - OUT$nu[s])*N2 + OUT$phase[j] - OUT$phase[s] ) -
Qprime(( OUT$nu[s] + OUT$nu[j])*N2)*cos((OUT$nu[j] + OUT$nu[s])*N2 + OUT$phase[j] + OUT$phase[s] )
)
}
}
epsilon = epsilon / Qsecond0 / N2 / OUT$amp[j]
OUT$nu[j] = OUT$nu[j] - epsilon
}
OUT <- mfft_analyse(xdata, nfreq, fast, nu = OUT$nu, my_minfreq, my_maxfreq)
}
# account for tseries scaling (i.e. uses the vaue of deltat and start encoded
# in the time series object
OUT$nu <- OUT$nu / dt
OUT$phase <- OUT$phase - startx*OUT$nu
# rearrange terms to avoid negative amplitudes and have phases in the right quandrant
# these are cosines so even functions
to_be_corrected <- which (OUT$amp < 0)
if (length(to_be_corrected)){
OUT$amp[to_be_corrected] <- - OUT$amp[to_be_corrected]
OUT$phase[to_be_corrected] <- OUT$phase[to_be_corrected] + pi
}
to_be_corrected <- which (OUT$nu < 0)
if (length(to_be_corrected)){
OUT$nu[to_be_corrected] <- - OUT$nu[to_be_corrected]
OUT$phase[to_be_corrected] <- - OUT$phase[to_be_corrected]
}
# finally, order termes by decreasing amplitude
o <- order(OUT$amp, decreasing = TRUE)
OUT$amp <- OUT$amp[o]
OUT$nu <- OUT$nu[o]
OUT$phase <- OUT$phase[o]
# rename for class compatibility
names(OUT) <- c("Freq","Amp","Phases")
attr(OUT, "class") <- "mfft_deco"
attr(OUT, "data") <- xdata
attr(OUT, "nfreq") <- nfreq
return(OUT)