I have a signal looking approximatively as the one in the first subplot below (*), and I would like to:

  • extract the periods of the main components of the signal;
  • associate an amplitude to these components.

(*) Here, the original signal is constructed from two main curves: a sine with a period ~170 and an amplitude ~7 and ~0.5 for the first and second half of the X space respectively; and a sine with a period ~20 and amplitude ~1.1 over the whole X space. Complete code reproduced below.

enter image description here

Because I suspect components at larger (>50) and smaller (<50) periods which I all want to catch, I first decompose my signal into a smoothed one (myySmooth for larger periods, in red) and its residuals (myySmallForms for smaller periods, in blue).

A wavelet analysis is conducted over these two signals (myySmooth, upper panels, and myySmallForms, lower panels; the right panels show the average power over the X axis as a function of the period):

enter image description here

The smaller period component is very well caught (19.7). The larger period component is reasonably well caught (main at period 184.82), and I can imagine that additional components are detected (here, around 131 and 321) given how the signal looks like (see code below: the X domain was split into 10 parts to generate this weird looking signal).

My question is:

How can I estimate the average/global amplitude of each of these components?

In other words: how can I find back that a component with a period around ~184.82 (for my initial ~170) should have an amplitude around 7, and that a component with a period ~19.7 (for my initial ~20) should have an amplitude around 1.1?

If I use the Fast Fourier Transform (with which I am more comfortable to compute the amplitude - see here for example), several components are found and I am afraid the actual amplitude is somehow "diluted" over them (?) --> here, a few components barely reach 2 for the larger periods (how to get back to ~7?) or 0.3 for the lower period (how to get back to ~1.1?)...

enter image description here


Reproducible example:

N = 2000
totTime = 2000
tstep = totTime/N
myx = seq(0, totTime-tstep, by=tstep)

n=10 # the X axis is split into 10 segments, the global signal is built  
# with different components within each of these segments
periofactor = 1/cbind(rnorm(n, 170, 15), rnorm(n, 20, 1)) 
amplitudes = rbind(cbind(rnorm(n/2, 7, 0.9), rnorm(n/2, 1.1, 0.03)),
                   cbind(rnorm(n/2, 0.5, 0.1), rnorm(n/2, 1.1, 0.1)))
periofactor = split(periofactor, 1:n)
amplitudes = split(amplitudes, 1:n)
myxsplit = split(myx, rep(1:n, each=N/n))

myy = c(mapply(function(amps, pers, x) 
    a*sin(p*x*2*pi), amps, pers)),
  amplitudes, periofactor, myxsplit)) + 0.5 * rnorm(N)

# Get a smoother to decompose the signal into a smooth one (larger periods) 
# and the residuals
myM = gam(Y ~ s(X, fx=F, bs='cr', k=100), data=data.frame(X=myx, Y=myy))
myySmooth = predict(myM, se=T, type='response')$fit
myySmallForms = myy - myySmooth

par(mfrow=c(2,1), mar=c(2, 3.5, 1.5, 0.5), ps = 8, mgp=c(0.8,0.1,0))
plot(myx, myy, type='l', main='Original signal')
plot(myx, myySmooth, type='l', col='red', lwd=2,
     main='Decomposed into myySmooth (red) and myySmallForms (blue)')
lines(myx, myySmallForms, col='blue')

Wavelet analysis:

getWavelets = function(mysignal) {
  mydata = data.frame(x=mysignal)
  myw = analyze.wavelet(mydata, "x", loess.span = 0, dt=1, dj=1/100,
                        lowerPeriod = 2, upperPeriod = 1024, n.sim = 10)
  wt.image(myw, n.levels=250,
           legend.params = 
             list(lab="wavelet power levels", mar=4.7,label.digits=3),
wavSmooth = getWavelets(myySmooth)
plot(wavSmooth$Period, wavSmooth$Power.avg)
wavSmallF = getWavelets(myySmallForms)
plot(wavSmallF$Period, wavSmallF$Power.avg)
mtext("myySmooth", side = 3, line = -1, outer = TRUE, col='red')
mtext("myySmallForms", side = 3, line = -13.5, outer = TRUE, col='blue')

FFT analysis:

par(mfrow=c(2,2), mar=c(2, 3.5, 1.5, 0.5), ps = 8, mgp=c(1.3,0.5,0))
myfft <- fft(myySmooth)
myfft <- myfft[2:((N / 2) + 1)]
freq <- (1:(N / 2)) / totTime
periods = 1/freq
ampl <- Mod(myfft) / length(myfft)
plot(periods, ampl, main='myySmooth', col.main='red')
myfft <- fft(myySmallForms)
myfft <- myfft[2:((N / 2) + 1)]
freq <- (1:(N / 2)) / totTime
periods = 1/freq
ampl <- Mod(myfft) / length(myfft)
plot(periods, ampl, main='myySmallForms', col.main='blue')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.