How to calculate 10 dB duration in R

Question

I'm working with dolphins echolocation clicks and I'm trying to get the duration of each click. One criteria used for the duration is to get the points where the signal declines 10dB to both sides from the peak frequency. I'm having trouble finding how to do this with R, do you have any clues about it?

Answer 1

If you are using PAMGuard click detectors, consider checking out the R-package called PAMpal, developed by Taiki Sakai with NOAA Southwest Fisheries.

PAMpal automatically calculates a number of useful features from click (and whistle) detections, including 10 dB bandwidth. Check out the user-guide here

Answer 2

This can be done fairly easy in R. Have a look at the function Q() in the "seewave" package, and the excellent article by Bennet-Clark:

Bennet-Clark, H. C. (1999). Which Qs to choose: questions of quality in bioacoustics?. Bioacoustics, 9(4), 351-359.

I show an example below with a synthetic rain forest rocket frog note. This is how the oscillogram and spectrogram of the note looks like:

If I understood correctly, you want to obtain the bandwidth for the complete signal. Therefore, the first step is to compute the mean power spectrum for the complete sound using the meanspec() function (the spec() function would also work if you are interested in the peak at a given time point of your signal):

spectroRf <- meanspec(Rf_synth, #replace by your sound here.
                  wl = 512,     #define your own windows length depending on your signal.
                  ovlp = 90,    #define your own overlap depending on your signal.
                  from = 0, to = 0.2,
                  flim = c(5,8),
                  dB = "max0") # scale the spectrum to max. 0 dB. This is necessary for using it later in the Q() function.

Then use this mean power spectrum as the first argument of the function Q(). Define the dB level below the peak that you are interested in (level = -10 in your case):

Q(spectroRf, 
   level = -10,
   xlim = c(5, 8),
   las = 1)

This will return both an image of the mean power spectrum with the computed Q, and a list containing the values of the minimum, maximum and peak frequencies together with the bandwidth:

    $Q
       x 
9.822804 

$dfreq
       x 
6.459972 

$fmin
       x 
6.118846 

$fmax
       x 
6.776497 

$bdw
        x 
0.6576505 

Hope this helps!

EDIT: I subscribe to the previous answers saying that defining the duration of a sound signal out of spectral analyses is not the best way to go. But I hope this post will help you find the frequency values 10dB below the peak of your clicks.

Answer 3

I have a different answer, which addresses finding signal duration at 10dB below the peak amplitude value of the oscillogram. Note the peak amplitude value of the oscillogram is NOT the same as the peak frequency of the power spectrum I described in my previous answer. This is the source of the confusion. The first answer is for those trying to find bandwidth values at 10dB below the peak of a power spectrum. This answer is for those trying to detect signal durations 10dB below the peak amplitude of the signal envelope, which seems to me was the original question.

Have a look at the function timer() in the package "seewave". This functions computes the envelope of the waveform of a signal (check the function env()) and uses a user defined threshold to find the begin and end times of the sound.

You need to define a windows length, and overlap percent, the type of windowing you will be using, and the amplitude threshold.

For the synthetic rocket frog signal, the duration of the signal at 10dB below the peak amplitude of the envelope can be obtained with:

detectRf <- timer(Rf_synth,
                  tlim = c(0, 0.2),
                  wl = 512,
                  ovlp = 90,
                  envt = "hil",
                  threshold = 31.62)
axis(2, las = 1)

This will return you a figure showing the envelope, and the begin and end limits of the detected signal in red, and a list with the exact values:

:

> detectRf
$s
[1] 0.01777979

$p
[1] 0.07960086 0.10261934

$r
[1] 0.09757312

$s.start
[1] 0.07960086

$s.end
[1] 0.09738066

$first
[1] "pause"

Note that the timer() function uses an user-defined amplitude threshold to automatically detect the begin and end time of the envelope of a signal (not on the raw waveform, therefore how you construct the envelope will influence the detection limits). The threshold is expressed in percent relative to the peak, which is normalised to a maximum value of 1 and expressed on linear scale. Therefore, for -10dB below the peak amplitude, the percent value can be obtained with the following arithmetic:

$$-10 dB = 20 \log_{10} \frac{X}{X_{ref}}$$

Where $X$ is the amplitude value you are looking for, and $X_{ref}$ is the reference value, which in this case is the peak amplitude and equals to 1.

This becomes

$10^{-0.5} = X$

$0.3162278 = X$

This is you threshold value at -10 dB expressed as a proportion, and therefore your threshold for the timer() function is 31.62%. Try it with 6dB and you will find your threshold to be 50%.

Hope this helps!

@Editors: feel free to join this answer with the previous one if you think its better. Although I consider they address different questions, and therefore should be two separate answers.

Answer 4

Maybe the following example may help

#Create a sequence of 100 equally spaced numbers between -4 and 4
x <- seq(-4, 4, length=100)

# simulate a signal
y <- dnorm(x)

# length based on 10 dB down from peak level
iu <- which(y>=max(y)/sqrt(10))

dev.new()
plot(x,y, type = "l")
lines(x[iu],y[iu],type="l",lwd=3,col="red")

print(max(x[iu])-min(x[iu]))

# length based on 10%, 90% of energy
z <- cumsum(y^2)
vu1<-approx(z,x,max(z)*0.1)
vu2<-approx(z,x,max(z)*0.9)

dev.new()
plot(x,z, type="l",col="black")
lines(vu1$y,vu1$x, type = "o", col = "blue")
lines(vu2$y,vu2$x, type = "o", col = "red")

print(vu2$y-vu1$y)

For this I simulated a signal and estimated the length by 10dB down (sqrt(10) down for amplitude) and the more common estimate 10 to 90 % of signal energy.

Edit: as the context of the question is still confusing (see comment of @Rasmus) I would like to point out that the one should not estimate the duration of a broadband signal in the spectral domain (or after narrow-band filtering) as the estimated duration includes the filter window length. For a good spectral estimation one needs a long window (1 Hz resolution at 48 kHz means 1 second window). Only if your signal is much longer than the FFT/DFT window, narrowband filtered duration estimates could make sense (+/- filterlength)

Caveat: I'm not an R programmer, so script can be improved.

Edit: For those of us that simulate nice signals, the Madsen RMS paper may be of interest on why -3dB or -10dB down is not always a good idea for real signals.