Calculating and Plotting Average Spectra
Basic Usage
PAMpal makes it easy to calculate the average spectrum for a set of clicks, since this is often a useful image to create for assisting in species classification. You simply need to provide PAMpal an AcousticStudy object, and tell it which event (or events) to create the image for. It will then take all the clicks from all channels in that event, calculate their spectra, and average them (converting to/from log space where appropriate). This function creates two plots - first a concatenated spectrogram of all the clicks, and second the average spectrum. Here we calculate and plot the average spectrum of the first event in our study.
avSpec <- calculateAverageSpectra(myStudy, evNum = 1)
Other ways to specify events
The evNum argument can be a vector to calculate and plot the average of multiple events. In this case the average will be calculated over all events, and the concatenated spectrogram will have black bars marking the start and end of each event. Here we see that our second event is much shorter and our third event is much longer.
calculateAverageSpectra(myStudy, evNum=1:3)
You can also specify the evNum argument as the name of the event you want to plot rather than the index number of that event. This is usually easier to do when you are trying to plot a specific event. Here we’ll also turn off the plotting option since the plot we are creating is identical to the first event, we are just specifying it by name this time. Vectors of names are also accepted.
avSpecName <- calculateAverageSpectra(myStudy, evNum="HB1603_MF_Master_Leg3_Skala.OE115", plot=FALSE)
# Double check that this is the same whether evNum = 1 or the name of event #1
identical(avSpec, avSpecName)## [1] TRUEDetails on the stored output, noise level, and calculation method
Note that calculateAverageSpectra also creates an output in addition to the plot (which is why we were able to compare those two examples to make sure they were the same). It stores all the data used to create the plots in the output.
str(avSpec)## List of 6
## $ freq : num [1:256] 375 750 1125 1500 1875 ...
## $ UID : chr [1:144] "1186" "1195" "1221" "1235" ...
## $ avgSpec : num [1:256] -43.1 -42.4 -42.1 -42.1 -41.9 ...
## $ allSpec : num [1:256, 1:144] -43.3 -46.6 -53.4 -48 -53.9 ...
## $ avgNoise: num [1:256] -43.1 -42.5 -42.1 -42.1 -41.9 ...
## $ allNoise: num [1:256, 1:144] -43.4 -46.6 -53.3 -48.1 -54.1 ...freq is the frequency value (Hz) on the axes of the plots (here the default window length is 512, so we have 256 frequency values), UID is the UIDs of the detections included in the plot, avgSpec is the computed average spectrum of those detections, and allSpec is a matrix of each individual spectrum used to calculate the average (each column corresponds to each UID). You’ll also note that there are similar outputs avgNoise and allNoise. calculateAverageSpectra internally attempts to compute a noise level for each clip, lets turn that option on and see what the plot shows. We’ll also make use of the ability to only plot one of the two default plots at a time since this will only show up on the average spectrum.
avNoise <- calculateAverageSpectra(myStudy, evNum=1, noise=TRUE, plot=c(FALSE, TRUE))
The estimated noise level shows up as a dotted line, but there’s a problem! The dotted line is exactly the same as the average spectrum. To figure out why, we need to talk about how PAMpal is calculating each individual spectra. For each detection, PAMpal only has access to the short click clip stored in the PAMGuard binary files. For the click spectra, it will take a window (default length is 512, can be specified with parameter wl) centered on the maximum value of the waveform. For the estimated noise spectra, it will take a window of the same length starting at the beginning of the clip. If the window length specified is longer than the duration of the clip stored in the binaries, the clip will be zero-padded until it is of appropriate length. Let’s take a look at one of the click clips to see what is going on here.
# Grab the binary data for the first click in this study
binData <- getBinaryData(myStudy, UID = getClickData(myStudy)$UID[1])
# plot the waveform
plot(binData[[1]]$wave[, 1], type='l')
The problem here is that a window of 512 centered around the peak of the waveform (~200) ends up being the exact same as the noise clip that starts at the beginning of the clip. This will also happen if your window length is too large and it results in zero-padding. The solution here is to shorten our wl argument so that our noise window will be independent from our click window. Looking at the plot, it seems like a wl of 128 should ensure this.
avNoise <- calculateAverageSpectra(myStudy, evNum=1, noise=TRUE, wl=128, plot=c(FALSE, TRUE))
Much better! That’s what we would expect the noise to be like. You might be wondering why we can’t just use the back half of that click clip to get a better estimate of noise with longer windows. The problem with this is that in some circumstances you might have an echo or reflection after the initial click, so taking the noise before the click is much more stable. The number of samples before the click can be changed within PAMGuard’s Click Detector settings.
SNR thresholding
There is an option to try and filter out clicks with a low signal-to-noise ratio (SNR). For each click, we estimate the SNR as the difference between the click spectrum and the noise spectrum at the frequency with the highest amplitude in the click spectrum. Since this depends on a relatively accurate noise spectrum, you must be sure that you have a reasonable noise line before trying to filter by SNR (ie. it looks like the one just above, and not the first plot we made). Once you are satisfied with that, then you can set the snr parameter, and all clicks with an SNR below that value will be removed from the results. Let’s see if we can remove some of the quieter clicks from one of our events.
par(mfrow=c(2,2))
calculateAverageSpectra(myStudy, evNum=3, wl=128, noise=TRUE, snr=0)
calculateAverageSpectra(myStudy, evNum=3, wl=128, noise=TRUE, snr=20)
Success! We removed about half of the clicks, and the noise level has dropped from ~-25dB to ~-30dB. Note that it is possible to set the snr parameter to too high of a value which might remove all clicks from an event, in which case you’ll get a warning and nothing will happen.
Sorting by peak frequency
Sometimes it can be easier to see patterns in the concatenated spectrogram if the clicks are sorted by peak frequency rather than the order in which they occur. Setting sort=TRUE changes the first plot to display in this manner, but changes none of the calculations. Note that if multiple events are being run at once, the black bars that would normally separate events are no longer present - all clicks are sorted as one group, rather than sorted within events.
calculateAverageSpectra(myStudy, evNum=1, sort=TRUE, plot=c(TRUE, FALSE))
calculateAverageSpectra(myStudy, evNum=1:3, sort=TRUE, plot=c(TRUE, FALSE))
Display options - Contrast, brightness, “q”, and axis limits
In some circumstances it can be tricky to create a concatenated spectrogram image with a satisfactory appearance. Since it is displaying a large number of clicks, one particularly loud or quiet click can throw off the scaling of the whole image. To help deal with this issue, there are two visual parameters contrast and brightness that adjust this image. They both range from -255 to 255, with positive values increasing contrast/brightness and negative values decreasing. There is usually a combination of these parameters that can overcome any image scaling issues with some trial and error. Note that these do not affect any of the numerical results, only the plots. The following are meant only as illustrative examples, the original plot is pretty decent.
calculateAverageSpectra(myStudy, evNum=1, plot=c(TRUE, FALSE), contrast=100)
calculateAverageSpectra(myStudy, evNum=1, plot=c(TRUE, FALSE), brightness=-40)
calculateAverageSpectra(myStudy, evNum=1, plot=c(TRUE, FALSE), contrast=100, brightness=50)
Apart from the brightness and contrast settings, there is also another graphical parameter you can tweak for the concatenated spectrogram, but its use is a little less obvious. This parameter is q, and it sets a quantile threshold for the plotted colors. The default value of q=.01 means that the bottom 1% and top 1% of values will be displayed with the same color value. Without this setting, a few extremely loud or extremely quiet clicks could throw off the scaling of the entire image, resulting in an image with apparently lower contrast. To illustrate this, here are plots with q values of 0%, 1%, and 5%.
par(mfrow=c(1, 3))
calculateAverageSpectra(myStudy, evNum=1, plot=c(TRUE, FALSE), q=0)
calculateAverageSpectra(myStudy, evNum=1, plot=c(TRUE, FALSE), q=.01) #default
calculateAverageSpectra(myStudy, evNum=1, plot=c(TRUE, FALSE), q=.05)
You can see that this has the effect of increasing the contrast in the image, but with higher q values you lose detail in the brightest and darkest points of the image.
There are also two options to adjust the displayed axis limits, flim and ylim. flim will adjust the frequency limits of both plots (units of Hz). ylim will adjust y-axis of the average spectra plot only.
calculateAverageSpectra(myStudy, evNum=1, plot=TRUE, flim=c(15e3, 80e3))
calculateAverageSpectra(myStudy, evNum=1, plot=c(FALSE, TRUE), ylim=c(-20, 0))
Decimation and Filters
There are a couple of signal processing options built in to the function. First, you can specify decimate to reduce the sample rate by integer factors (only integers are allowed). This uses the decimate function from the signal package and applies a lowpass filter to avoid aliasing.
calculateAverageSpectra(myStudy, evNum=1, decimate=2)
You can also specify high- or band-pass filters by settings filterfrom_khz and filterto_khz parameters (low-pass currently not supported). It is usually necessary to combine these with the ylim parameter in order to create a useful average spectrum plot since they will greatly reduce the magnitude of the signal.
calculateAverageSpectra(myStudy, evNum=1, filterfrom_khz = 10, plot=c(TRUE, FALSE))
# average spectrum scale will be way off
calculateAverageSpectra(myStudy, evNum=1, filterfrom_khz = 10, filterto_khz = 60)
calculateAverageSpectra(myStudy, evNum=1, filterfrom_khz = 10, filterto_khz = 60, ylim=c(-30, 0), plot=c(FALSE, TRUE))
Cepstrum
One more obscure feature of the calculateAverageSpectra function is that it can calculate the cepstrum of each click instead of the spectrum. In most cases (this example included) this is not particularly useful, but it can be used to find things like echoes of clicks or estimate the IPI in sperm whale clicks. Due to the way that the cepstrum works when finding these features, it is often necessary to apply a decimation factor, increase the window length, or both. This example is purely illustrative of how you would do that, but results are not likely to be meaningful.
calculateAverageSpectra(myStudy, evNum=1, mode='ceps', decimate=2, wl=1024)
Miscellaneous other parameters
channel can limit the calculations to only considering the first or second channel by setting channel=1 or channel=2, the default is to use both channels for each detection and average them.
sr can be used to manually set the samplerate for detections. This is sometimes necessary if data was decimated in Pamguard, but is not needed if you had supplied an XML settings file with addSettings before processing your data in PAMpal
calibration can be used to provide a frequency response curve calibration after calcuating the spectrum. This does not need to be specified if one was added before processing your data with addCalibration (recommended to do this if a calibration is needed/available), in this case it will detect that calibration and use it.
norm is a logical flag (default TRUE) to normalize the magnitude of the average spectra to 0. This is only shifting the displayed curve, not affecting any of the calculations.