www.astrostatistics.org/projects/meteors

An implementation in R is given of the meteor detection method developed by Tom Roelandts. More information on the method can be found in Roelandts’ IMC 2014 paper and on his website which includes a python implementation of the method.

While generally applicable, the parametrisation and file naming conventions of the Belgian RAdio Meteor Stations (BRAMS) system are assumed in the present implementation of the detection method. BRAMS is a network of receiving radio stations and two emitting beacons. Forward scattering techniques are used to study the observed meteor population. The BRAMS system is operated by the Belgian Institute for Space Aeronomy (BISA). More information can be found on the BRAMS and BISA websites.

Data recorded at the receiving stations can be downloaded freely from the BRAMS website, after simple registration as a user. The data available for download are either .wav files containing the signal received, or .png image files of spectrograms. No meteor detections or counts are available.

The R script presented here demonstrates how meteors can be detected from the BRAMS .wav files using Roelandts’ method.

The detection function

The following code chunck shows the implementation of the entire detection method into a single function returning the meteor detections and a spectrogram plot.

Scroll down to after the function definition to see it in action.

detectMeteors = function(fileName,          # wav file following BRAMS naming conventions
                         Fmin = 0.8,        # min range to seek carrier in kHz
                         Fmax = 1.2,        # max range to seek carrier in kHz
                         passBw2 = 30,      # half bandwidth of bandpass filter in Hz
                         raShort = 101,     # length running average short filter
                         raLong = 30001,    # length running average long filter
                         thresh = 0.025,    # detection threshold for indicator signal
                         minBetween = 5512)  # minimum number of samples between meteors
{
   # Load required libraries
   require(tuneR)
   require(seewave)
   require(caTools)
   require(ggplot2)
   require(reshape)
   require(RColorBrewer)
   
   # Read wav file, normalize, and get properties
   W = readWave(fileName)
   WN = normalize(W)
   Fs = attr(W,"samp.rate") # sampling frequency
   n = length(W)
   duration = n / Fs
   fileComp = unlist(strsplit(fileName, "_", fixed="TRUE"))
   bedouri = which(fileComp=="BEDOUR")
   thisDate = as.integer(fileComp[bedouri+1])
   thisTime = as.integer(fileComp[bedouri+2])
   thisHour = thisTime %/% 100
   thisMin = thisTime %% 100
   thisLoc = fileComp[bedouri+3]
   
   # Find carrier frequency Fc
   Wx = W[1:2 ** floor(log2(n))]  # discard end of W so to obtain length power of 2
   nn = (2 ** floor(log2(Fs * 100)))
   Wx = W[1:nn] 
   WxN = normalize(Wx)
   Wspec = as.data.frame(spec(WxN, wl = nn, plot = FALSE))
   WspecBand = subset(Wspec,x>Fmin & x <Fmax)
   rowMax = which.max(WspecBand$y)
   Fc = 1000 * WspecBand$x[rowMax]  # carrier freq in Hz
   rm(Wx, nn, WxN, Wspec, WspecBand, rowMax)
   gc()
   
   # Band pass filter
   sigbp = as.numeric(ffilter(WN, from = Fc - passBw2, to = Fc + passBw2, 
                              wl = 4000, ovlp = 90, wn = "blackman"))
   sigbp = c(sigbp, rep(0,length(W) - length(sigbp))) # pad with zeros to full original length
   
   # Define indicator signal
   pow = as.numeric(sigbp)^2 # power signal
   powS = runmean(pow, k = raShort, alg = "fast", endrule = "constant", align = "center")
   powL = runmean(pow, k = raLong, alg = "fast", endrule = "constant", align = "center")
   rmsamp = (raLong + 1) / 2  # samples to remove at ends due to running average windowing
   pow = pow[(rmsamp+1):(n-rmsamp)]
   powS = powS[(rmsamp+1):(n-rmsamp)]
   powL = powL[(rmsamp+1):(n-rmsamp)]
   powL[powL < 1e-5] = 1e-5 
   indSig = (powS * raShort) / (powL * raLong)
   rm(powS, sigbp, WN)
   gc()
   
   # Detect meteors through thresholding of the indicator signal
   met = (indSig > thresh)
   metPlus = c(FALSE,met[-length(met)]) # shifted right one position
   metStart = met & (!metPlus) # start indices of TRUE sequence (=meteor)
   metStop = metPlus & (!met)  # end indices of TRUE sequence
   atLeast1mtr = sum(metStart) > 0
   if (atLeast1mtr) {
      i = 1:length(met)
      starti = i[metStart]
      stopi = i[metStop]
      if (length(starti) != length(stopi)) {
         stop("Meteor detection failed.")
         # @@@ try to fix the problem!
      }
      rm(i)
   }
   rm(met, metPlus, metStart, metStop)
   gc()
   
   # Collapse intervals that are not separated enough to define meteors
   metDf = data.frame(starti=NA,stopi=NA)
   if (atLeast1mtr) {
      valid = stopi >= starti # only keep valid intervals
      starti = starti[valid]
      stopi = stopi[valid]
      merge = FALSE
      counter = 1
      for (j in 1:length(starti)) {
         if (!merge) {
            metDf[counter,"starti"] = starti[j]
            metDf[counter,"stopi"] = stopi[j]
         } else {
            metDf[counter,"stopi"] = stopi[j]
         }
         if ( (j < length(starti)) & (starti[j+1] < stopi[j] + minBetween) ) {
            merge = TRUE
         } else {
            merge = FALSE
            counter = counter + 1
         }
      }
      metDf$mid = metDf$starti + floor((metDf$stopi - metDf$starti) / 2)
      metDf$date = thisDate
      metDf$starttime = thisTime
      metDf$hour = thisHour
      metDf$station = thisLoc
      
      # Define additional properties of meteors 
      for (j in 1:nrow(metDf)) {
         metDf$centre[j] = metDf$starti[j] + which.max(indSig[metDf$starti[j]:metDf$stopi[j]]) - 1
         thisRange = (metDf$starti[j]):(metDf$stopi[j])
         metDf$duration[j] = (metDf$stopi[j] - metDf$starti[j]) / Fs
      }
      # offset all indices with removed number of samples, rmsamp, and find corresponding times
      metDf$midadj = metDf$mid + rmsamp
      metDf$centreadj = metDf$centre + rmsamp
      metDf$time = metDf$centreadj / Fs
   } 
   rm(indSig, pow, powL)
   gc()
   
   # Create spectrogram plot, similar to BRAMS plots
   specWinLen = 16384
   Wspec = spectro(W, plot=FALSE, wl = specWinLen, ovlp = 80, complex = FALSE, norm = TRUE, dB = NULL)
   freqMin = 0.001 * (Fc - 100) # select band of 200 Hz around carrier
   freqMax = 0.001 * (Fc + 100)
   freqN = length(Wspec[["freq"]])
   freqAll = Wspec[["freq"]]
   myFreqIndex = (1:freqN)[freqAll > freqMin & freqAll < freqMax]
   myFreq = freqAll[myFreqIndex]
   myMat = Wspec[["amp"]][myFreqIndex,]
   myTime = Wspec[["time"]]
   dbVec = 10 * log10(as.vector(t(myMat)))
   X = data.frame(time=rep(myTime,length(myFreq)),
                  frequency=rep(myFreq,each=length(myTime)),
                  power=dbVec)
   
   # Adjust meteor times to time range in frequency domain
   if (atLeast1mtr) {
      metDf$timeAmp = metDf$time
      durationFreq = duration - specWinLen / Fs
      timeShift = (specWinLen / 2) / Fs
      timeScale = duration / durationFreq
      metDf$time = timeScale * (metDf$timeAmp - timeShift)
   }
   
   myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")), space="Lab")
   myColors = myPalette(100)
   plotRange = c(-20, 0)
   restrictRange = function(v, range) {
      v[v < range[1]] = range[1]
      v[v > range[2]] = range[2]
      return(v)
   }
   X$power = restrictRange(X$power, plotRange)
   yBreaks = seq(freqMin, freqMax,by = 0.02)
   yLabels = as.character(round(1000 * yBreaks,0))
   fillBreaks = seq(plotRange[1],plotRange[2],length.out = 6)
   fillLabels = as.character(round(fillBreaks,1))
   p = ggplot(data=X, aes(x = time, y = frequency)) + 
      geom_tile(aes(fill = power), limits = plotRange) +
      scale_fill_gradientn(colours = myColors, 
                           limits = plotRange,
                           breaks = fillBreaks,
                           labels = fillLabels) + 
      scale_x_continuous(limits=c(myTime[1], myTime[length(myTime)]), 
                         expand = c(0,0), 
                         breaks = seq(0,300,by = 30)) + 
      scale_y_continuous(limits=c(freqMin, freqMax), 
                         expand = c(0,0), 
                         breaks = yBreaks,
                         labels = yLabels) +
      xlab("Time (s)") + 
      ylab("Frequency (Hz)") +
      guides(fill = guide_colorbar(barwidth = 0.5, barheight = 12, title = "Power (dB)")) +
      ggtitle(paste(thisLoc, ", ", 
                    thisDate, ", ",
                    sprintf("%02d", thisHour), ":",
                    sprintf("%02d", thisMin),
                    sep=""))
   if (atLeast1mtr) {
      p = p + geom_vline(data=metDf, aes(xintercept=time),colour="red") 
   }
   return(list(meteors = metDf, meteorPlot = p))
}

Applying the detection method

The function can be used on .wav files obtained from the BRAMS website. Here we use the same example file as that used by Tom Roelandts on his website.

Provide the file path,

filePath = "C:/data/BRAMS/example"  # change as appropriate

and call the detection routine:

theFile = paste(filePath, "RAD_BEDOUR_20111007_0420_BEUCCL_SYS001.wav", sep="/" )
M = detectMeteors(theFile)

Inspecting results

The result is the object M, a list of length 2.

class(M)
## [1] "list"
length(M)
## [1] 2
names(M)
## [1] "meteors"    "meteorPlot"

First look at the list of meteors and some of their properties. Some are not very useful but are stored in the function as intermediate results and possibly future use.

M$meteors[, c("station", "date", "starttime", "time", "duration")]
##    station     date starttime      time    duration
## 1   BEUCCL 20111007       420  19.79573 0.063486305
## 2   BEUCCL 20111007       420  33.99398 0.133321241
## 3   BEUCCL 20111007       420  41.04923 0.076909124
## 4   BEUCCL 20111007       420  53.86740 0.024668964
## 5   BEUCCL 20111007       420  99.39930 0.024487575
## 6   BEUCCL 20111007       420 111.48061 0.009795030
## 7   BEUCCL 20111007       420 138.29660 0.077090513
## 8   BEUCCL 20111007       420 149.61468 0.069653546
## 9   BEUCCL 20111007       420 222.01625 0.006167241
## 10  BEUCCL 20111007       420 288.46071 0.022855070

A total of 10 meteors are detected. The time column gives the time in seconds since 4:20 on 2011-10-07 for the Uccle station (code BEUCCL).

Now look at the spectrogram. The red vertical lines indicate the meteor detections.

M$meteorPlot

When applying the detectMeteors() function to data from other stations, some of the parameters may need adjusting. At present the parameter values are set on an ad-hoc basis.

Check here for more on meteor detections: www.astrostatistics.org/projects/meteors