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 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))
}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 appropriateand call the detection routine:
theFile = paste(filePath, "RAD_BEDOUR_20111007_0420_BEUCCL_SYS001.wav", sep="/" )
M = detectMeteors(theFile)The result is the object M, a list of length 2.
class(M)## [1] "list"length(M)## [1] 2names(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.022855070A 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