Sample Rate Reduction for Hydrophone Data

Paul McGill
July 29, 2019

Calculating decimation factors

Our goal is to decimate the 256 kHz hydrophone data down to 2 kHz, a factor of 128. This will require first low-pass filtering the signal to limit its bandwidth, followed by downsampling (i.e. throwing out 127 out of every 128 samples). Richard Lyons describes the process in [1], section 10.1.

Instead of decimating by a factor D of 128 in one stage, it’s more efficient to use two stages, D1 and D2, resulting in much shorter filters. The optimum factor for the first stage D1 is defined in Lyons as:

dopt[D_,F_]:=2 D ((1-(Sqrt[D F/(2-F)]))/(2-F(D+1)))

where D = D1 x D2 is the total decimation factor, and F = Δf/fstop is the ratio of the final filter transition region width over the stopband frequency. The stopband frequency fstop is 1 kHz, the Nyquist rate of the final 2 kHz sample rate. If we set Δf to 100 Hz, which is 10% of the final bandwidth, then F = 100 Hz/1 kHz = 0.1.

dopt[128,0.1]

We can now set D1 equal to the closest integer submultiple of D,

Divisors[128]

which is D1 = 32.

d1=32;d2=128/d1

Calculating filter parameters

Now we must find the FIR filter lengths and frequencies for each of the two stages. First define the original sample rate.

sr=256000

For the first stage, we want to start filtering at 90% of the final 1 kHz bandwidth.

f1cut=0.9 1000

We want the stopband frequency at the Nyquist rate after the first stage decimation.

f1stop=((sr/d1)/2)

Find the transition bandwidth of the first-stage filter, in radians. Note that this filter is applied to the 256 kHz data, so the Nyquist rate is 128 kHz (sr / 2).

tbw1=((f1stop-f1cut)/((sr/2))) Pi

Find the minimum length of the Blackman filter kernel for that transition bandwidth. Use the table on page 9 of [2].

m1=Ceiling[((12 Pi)/tbw1)]

For the second stage, we start filtering at the same frequency as the first stage.

f2cut=f1cut

We want the stopband frequency at the Nyquist rate after the second stage.

f2stop=((sr/(d1 d2))/2)

Find the transition bandwidth of the second-stage filter, in radians. Note that this filter is applied to the 8 kHz output of the first stage, so the Nyquist rate is 4 kHz (f1stop).

tbw2=((f2stop-f2cut)/f1stop) Pi

Find the minimum length of the Blackman filter kernel for that transition bandwidth.

m2=Ceiling[((12 Pi)/tbw2)]

In summary, here are the specifications for the two stages. The cutoff frequency and transition bandwidth is specified in radians over a range of 0 to Pi, with Pi being the Nyquist frequency.

Grid[{{"D1","Fc1","Tbw1","M1","D2","Fc2","Tbw2","M2"},{d1,f1cut  2 Pi/sr//N,tbw1,m1,d2,f2cut 2 Pi/(sr/d1)//N,tbw2,m2}},Frame->All,Dividers->{{True,True,True,True,Thickness[3],True,True,True},{True,Thickness[3],True}}]

Calculating filter kernels

Now design filter kernels for the two stages of decimation. We want to use windowed sinc low-pass filters because they’re very stable, linear-phase, and easy to calculate. Although Mathematica has many filter design functions, windowed sinc is not one of them and so we must write one. Steven Smith describes the filter design process in [3], chapter 16. We’ll use a Blackman window, which has linear phase, 0.02 % ripple in the passband, and provides 74 dB of attenuation in the stopband. The filter cutoff frequency ω is specified to be between 0 and Pi representing 0 Hz to the Nyquist frequency. The length of the kernel produced is m+1 elements.

windowedSincFilterKernel[ω_,m_]:=Module[{kern,win},
(* calculate sinc kernel *)
kern=Table[Sinc[ω x],{x,-m/2,m/2}];
(* calculate window *)
win=Array[BlackmanWindow,m+1,{-0.5,0.5}];
(* calculate windowed sinc kernel and normalize for 0 dB gain at DC *)
kern win/Total[kern win]
] (* end Module *)

Note that a Python implementation of windowed sinc filters can be found in [4].

Calculate the decimation filter kernels for the two stages of decimation and look at their responses.

{f1cut,f1stop,m1}

h1=windowedSincFilterKernel[2 Pi f1cut/sr,m1];

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi f/sr]]],{f,0,f1stop},PlotRange->{{0,f1stop},{-120,0}},Frame->True,GridLines->Automatic,FrameLabel->{{"Attenuation (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}}]

We’d like to have very little attenuation below, say, 850 Hz in order to preserve signal amplitudes in the final passband. We can zoom in to the upper left portion of the filter response to evaluate this better.

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi f/sr]]],{f,0,f1stop},PlotRange->{{0,1000},{-10,0}},Frame->True,GridLines->Automatic,FrameLabel->{{"Attenuation (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}}]

We’re about 5 dB down at 850 Hz, which is a bit too much. We can raise the cutoff frequency of the filter as long as we still have at least 74 dB of attenuation at 4 kHz.

f1cut=2000

h1=windowedSincFilterKernel[2 Pi f1cut/sr,m1];

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi f/sr]]],{f,0,f1stop},PlotRange->{{0,f1stop},{-120,0}},Frame->True,GridLines->Automatic,FrameLabel->{{"Attenuation (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}}]

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi f/sr]]],{f,0,f1stop},PlotRange->{{0,1000},{-10,0}},Frame->True,GridLines->Automatic,FrameLabel->{{"Attenuation (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}}]

Now we have very little attenuation anywhere in the final passband. This means we probably could’ve used a shorter FIR filter but we’ll leave it as is for now.

{f2cut,f2stop,m2}

h2=windowedSincFilterKernel[Pi f2cut/f1stop,m2];

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h2,Pi f/f1stop]]],{f,0,f1stop},PlotRange->{{0,f2stop},{-120,0}},Frame->True,GridLines->Automatic,FrameLabel->{{"Attenuation (dB)",None},{"Frequency (Hz)","Stage 2 Lowpass Filter Response"}}]

Once again we’ll zoom into the upper portion of the graph to verify that there’s not much attenuation in the passband.

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h2,Pi f/f1stop]]],{f,0,f1stop},PlotRange->{{0,f2stop},{-10,0}},Frame->True,GridLines->Automatic,FrameLabel->{{"Attenuation (dB)",None},{"Frequency (Hz)","Stage 2 Lowpass Filter Response"}}]

Here are the modified filter parameters.

Grid[{{"D1","Fc1","Tbw1","M1","D2","Tbw2","Fc2","M2"},{d1,(f1cut  2 Pi/sr)//N,tbw1,m1,d2,f2cut 2 Pi/(sr/d1),tbw2,m2}},Frame->All,Dividers->{{True,True,True,True,Thickness[3],True,True,True},{True,Thickness[3],True}}]

Filtering hydrophone data

Now we can apply this decimation to data from the hydrophone to produce 2 kHz data. Import the second minute of a ten-minute data file that contains whale calls.

whales=Take[Import["/Users/mcgill/Documents/My\ Projects/Passive\ Acoustic\ Monitoring/MARS_20171003_200022.wav","Data"][[1]],
{60 sr,120 sr}];

Look at both the spectrogram and periodogram of the data. Note that we must amplify the data to see it on the spectrogram. Also note the attenuation of frequencies above 100 kHz due to the antialiasing filter in the hydrophone.

Spectrogram[1 (10^4) whales,SampleRate->256000,PlotRange->All]

Periodogram[whales,1024,SampleRate->256000,PlotRange->{All,{-120,-20}},Frame->True,GridLines->Automatic]

Filter and downsample the data, with the first stage producing 8 kHz data and the second stage producing 2 kHz data.

whales8kHz=
Downsample[
ListConvolve[h1,whales],
d1];

whales2kHz=
Downsample[
ListConvolve[h2,whales8kHz],
d2];

Spectrogram[5 (10^4) whales8kHz,1024,SampleRate->8000,PlotRange->All,Frame->True,FrameLabel->{{"Time (s)","Stage 1 Output Spectrogram"},{"Frequency (Hz)",""}}]

Periodogram[whales8kHz,1024,SampleRate->8000,PlotRange->{All,{-120,-20}},Frame->True,GridLines->Automatic]

Spectrogram[5 (10^4) whales2kHz,1024,SampleRate->2000,PlotRange->All,Frame->True,FrameLabel->{{"Time (s)","Stage 2 Output Spectrogram"},{"Frequency (Hz)",""}}]

Periodogram[whales2kHz,1024,SampleRate->2000,PlotRange->{All,{-120,-20}},Frame->True,GridLines->Automatic]

Verifying the results

The downsampled data clearly show the whale calls and the attenuation above 900 Hz, but how do we know that the low-pass filters are properly preventing aliasing from strong signals just above the passband? To test this create a 1200 Hz tone at 256 ksps that lasts one minute.

tone=Table[Sin[2 Pi x 1200/sr],{x,1,Length[whales]}];

ListPlot[tone[[;;213]]]

Attenuate the tone and add it to the whale data.

whaletone=whales+(10^(-3)) tone;

Look at the spectrum of the combined signal.

Periodogram[whaletone[[;;256000]],SampleRate->256000,PlotRange->{{0,2000},{-120,0}},Frame->True,GridLines->Automatic]

The tone is very strong at about 70 db louder than the natural signals we see. Now do the same analysis on the data containing the tone and look at the resulting spectrum.

whaletone8kHz=
Downsample[
ListConvolve[h1,whaletone],
d1];

whaletone2kHz=
Downsample[
ListConvolve[h2,whaletone8kHz],
d2];

Spectrogram[5 (10^4) whaletone2kHz,1024,SampleRate->2000,PlotRange->All,Frame->True,FrameLabel->{{"Time (s)","Stage 2 Output Spectrogram\rwith 1200 Hz Tone"},{"Frequency (Hz)",""}}]

Periodogram[whaletone2kHz,1024,SampleRate->2000,PlotRange->{All,{-120,-20}},Frame->True,GridLines->Automatic]

The spectrogram shows no hint of a tone at 800 Hz, which is where the 1200 Hz tone would appear if it had not been sufficiently filtered.

Implementation notes

The analysis above doesn't address three issues that might be important in a practical implementation:

1) Although most implementations of  windowed sinc filters are linear phase, they're not zero phase, and therefor will introduce a delay in the output signal equal to one-half the filter length. If exact timing is important for subsequent analysis, then the output samples should be shifted earlier in time accordingly. See [3], chapter 19, page 328 for an explanation of filter phase response.

2) As with any digital filtering operation, the effects at the beginning and end of each file must be considered. In the analysis above, I didn't worry about this and there is no overhang between the filter kernel and either end of the data. To process long, continuous segments of data in multiple files, the overlap-add method should be used to ensure that no data are lost at the beginning and end of each file. A description of overlap-add can be found in [3], chapter 18.

3) The filtering operations in the analysis above are very inefficient. For instance in the case of decimation by a factor of 32, every output value of the filter is calculated in the filtering step. This is immediately followed by the downsampling step, where 31 out of 32 filter output values are discarded. It's much more efficient to calculate only every 32nd output value, which is possible because FIR filters don't depend on past filter outputs. This requires additional bookkeeping overhead that may not save time for small decimation factors.

References

[1] Lyons, R. G. (2004). Understanding digital signal processing (2nd ed.). Upper Saddle River, NJ: Prentice Hall.
[2] Notes on FIR filter design using window functions. (n.d.). Retrieved March 6, 2018, from  http://course.ece.cmu.edu/~ece396/lectures/L19/WindowFIRDesign.pdf
[3] Smith, S. W. (1997). The scientist and engineers guide to digital signal processing. San Diego (California): California Technical Publishing. (accessible for free at http://www.dspguide.com/pdfbook.htm)
[4] How to Create a Simple Low-Pass Filter. (n.d.). Retrieved March 06, 2018, from https://tomroelandts.com/articles/how-to-create-a-simple-low-pass-filter