Sample Rate Reduction for Geo-Sense DAS Data

Paul McGill
May 14, 2025

Calculating decimation factors

Our goal is to decimate the 1 kHz DAS data down to 1 Hz, a factor of 1000. The 1 Hz data can then be processed to look for long-period effects such as internal waves, infragravity waves, and tides. This will require first low-pass filtering the signal to limit its bandwidth, followed by downsampling (i.e. throwing out 999 out of every 1000 samples). Richard Lyons describes the process in [1], section 10.1.

Instead of decimating by a factor D of 1000 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 (Eq 10-2) 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 0.5 Hz, the Nyquist rate of the final 1 Hz sample rate. If we set Δf to 0.05 Hz, which is 10% of the final bandwidth, then F = 0.05 Hz/0.5 Hz = 0.1.

dopt[1000,0.1]

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

Divisors[1000]

which is D1 = 125.

d1=125;d2=1000/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 going into the first-stage filter, and the sample rate going into the second-stage filter.

sr1=1000; sr2=sr1/d1

We’ll need the stopband frequency of the second stage,  which is the Nyquist rate after the second stage.

f2stop=((sr2/d2)/2)//N

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

f1cut=0.9 f2stop

We need to set the stopband frequency of the first stage (f1stop) low enough to prevent aliasing in the second-stage output, but as high as possible to make the transition band wide, thus requiring a shorter filter. The stopband frequency of the first stage can be higher than the Nyquist rate after the first stage decimation, because we can allow aliasing in the band between f2stop (0.5 Hz) and 7.5 Hz, which is the sample rate after the first stage (8 Hz) minus the Nyquist after the second stage (0.5 Hz). For an explanation of this see Lyons Fig. 10-4.

f1stop=sr2-f2stop

Find the transition bandwidth of the first-stage filter, in radians. Note that this filter is applied to the 1 kHz data, so the transition bandwidth is calculated with respect to the Nyquist rate of 500 Hz (sr1/ 2).

tbw1=((f1stop-f1cut)/((sr1/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

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

sr2/2

tbw2=((f2stop-f2cut)/(sr2/2)) 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 are
 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/sr1,tbw1,m1,d2,f2cut 2 Pi/sr2,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/sr1,m1];

ListPlot[h1]

sr1

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi f/sr1]]],{f,0,8},PlotRange->{{0,8},{-150,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Gain (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}},ImageSize->Full]

We’d like to have very little attenuation below 0.45 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/sr1]]],{f,0,1},PlotRange->{{0,1},{-3,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Gain (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}},ImageSize->Full,Epilog->{Gray,Dashed,Line[{{0.45,-3},{0.45,0}}]}]

We’re about 0.63 dB down at 0.45 Hz, which is acceptable, but we should have a bit more attenuation at 7.5 Hz (at least 74 dB). Increase the length of the stage 1 filter a bit to achieve this.

f1cut=0.45

m1=1024

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

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi f/sr1]]],{f,0,8},PlotRange->{{0,sr2},{-120,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Attenuation (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}},ImageSize->Full,Epilog->{Gray,Dashed,Line[{{0,-74},{8,-74}}]}]

Find the attenuation at 7.5 Hz.

20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi 7.5/sr1]]]

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi f/sr1]]],{f,0,1},PlotRange->{{0,1},{-3,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Gain (dB)",None},{"Frequency (Hz)","Stage 1 Lowpass Filter Response"}},ImageSize->Full,Epilog->{Gray,Dashed,Line[{{0.45,-3},{0.45,0}}]}]

Find the attenuation at 0.45 Hz.

20 Log[10,Abs[ListFourierSequenceTransform[h1,2 Pi 0.45/sr1]]]

Now we have acceptable attenuation and aliasing in the final passband.

Design the second stage filter.

{f2cut,f2stop,m2}

h2=windowedSincFilterKernel[2 Pi f2cut/sr2,m2];

ListPlot[h2,Joined->True,PlotRange->All]

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h2,2 Pi f/sr2]]],{f,0,8},PlotRange->{{0,f2stop},{-120,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Gain (dB)",None},{"Frequency (Hz)","Stage 2 Lowpass Filter Response"}},Epilog->{Gray,Dashed,Line[{{0.45,-120},{0.45,0}}]}]

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,2 Pi f/sr2]]],{f,0,sr2/2},PlotRange->{{0.4,f2stop},{-5,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Gain (dB)",None},{"Frequency (Hz)","Stage 2 Lowpass Filter Response"}},Epilog->{Gray,Dashed,Line[{{0.45,-120},{0.45,0}}]}]

The attenuation in the passband is too much, so we can adjust the cutoff frequency to be slightly higher as long as we have sufficient attenuation at the stopband frequency f2stop.

f2cut=0.47

h2=windowedSincFilterKernel[2 Pi f2cut/sr2,m2];

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h2,2 Pi f/sr2]]],{f,0,1},PlotRange->{{0,1},{-120,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Gain (dB)",None},{"Frequency (Hz)","Stage 2 Lowpass Filter Response"}},Epilog->{Gray,Dashed,Line[{{0.45,-120},{0.45,0}}]}]

Find the attenuation at 0.5 Hz.

20 Log[10,Abs[ListFourierSequenceTransform[h2,2 Pi 0.5/sr2]]]

Plot[20 Log[10,Abs[ListFourierSequenceTransform[h2,2 Pi f/sr2]]],{f,0,sr2/2},PlotRange->{{0.4,f2stop},{-3,0}},Exclusions->None,Frame->True,GridLines->Automatic,FrameLabel->{{"Gain (dB)",None},{"Frequency (Hz)","Stage 2 Lowpass Filter Response"}},ImageSize->Full,Epilog->{Gray,Dashed,Line[{{0.45,-120},{0.45,0}}]}]

Find the attenuation at 0.45 Hz.

20 Log[10,Abs[ListFourierSequenceTransform[h2,2 Pi 0.45/sr2]]]

Here are the modified filter parameters.

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

Filtering Geo-Sense DAS data

Now we can apply this decimation to DAS data to produce 1 Hz data. We’ll do this in two stages, processing every file to produce an equal number of files containing data sampled at 8 Hz, and then processing those files to produce a single file sampled at 1 Hz.

First import the function to read Sintela raw data files.

sintelaReadRawV3[fileName_String,OptionsPattern[{firstChan->1,lastChan->-1,startTime->{1970},endTime->{2050},verbose->False}]]:= Module[ (* v1.2.0 of this function *)
{fileSizeBytes,inStream,syncWord,hdrSize,hdrVers,sampCnt,sampTime,numChans,numSamps,sampRate,chanSpacing,gaugeLen,trigger,startChan,chanStep,gpsPPSTime,gpsPPSSampOffs,gpsStat,trigFlags,statFlags,packetSizeBytes,totalPackets,totalSamps,packet,sampTimeOffset=20,startTimeSecs,endTimeSecs,rdBuf,data={},firstPacket=True,firstSampTime},

(* Decode header information *)
fileSizeBytes=UnitConvert[FileSize[fileName],"Bytes"]//QuantityMagnitude;
inStream=OpenRead[fileName,BinaryFormat->True];
syncWord=BinaryRead[inStream, "UnsignedInteger32"] (* offset 0 *);
hdrSize=BinaryRead[inStream, "UnsignedInteger32"] (* offset 4 *);
hdrVers=BinaryRead[inStream, "UnsignedInteger32"] (* offset 8 *);
If[OptionValue[verbose],
If[
syncWord==16^^11223344&&hdrSize==96,
Print["Sintela Data Type 3 - Raw 32 Bit."], 
(Print["Not a Sintela Type 3 File."];Close[inStream];Abort[])
] ;(* end if *)
]; (* end if *)
sampCnt=BinaryRead[inStream, "UnsignedInteger64"] (* offset 12 *);
sampTime=BinaryRead[inStream,"UnsignedInteger64"] (* offset 20 *);
numChans=BinaryRead[inStream,"UnsignedInteger32"] (* offset 28 *);
numSamps=BinaryRead[inStream,"UnsignedInteger32"] (* offset 32 *);
sampRate=BinaryRead[inStream,"Real32"] (* offset 36 *);
chanSpacing=BinaryRead[inStream,"Real32"] (* offset 40 *);
gaugeLen=BinaryRead[inStream,"Real32"] (* offset 44 *);
trigger=BinaryRead[inStream,"UnsignedInteger32"] (* offset 48 *);
startChan=BinaryRead[inStream,"UnsignedInteger32"] (* offset 52 *);
chanStep=BinaryRead[inStream,"UnsignedInteger32"] (* offset 56 *);
gpsPPSTime=BinaryRead[inStream,"UnsignedInteger32"] (* offset 60 *);
gpsPPSSampOffs=BinaryRead[inStream,"Integer32"] (* offset 64 *);
gpsStat=BinaryRead[inStream,"UnsignedInteger32"] (* offset 68 *);
trigFlags=BinaryRead[inStream,"UnsignedInteger8"] (* offset 72 *);
statFlags=BinaryRead[inStream,"UnsignedInteger8"] (* offset 73 *);

(* Display first header *)
If[OptionValue[verbose],
Print["Packet Channels x Samples : ",numChans," x ",numSamps];Print["Start Channel: ",startChan];
Print["Fs: ",NumberForm[sampRate,{8,3}]," Hz"] ;Print["Channel Spacing: ",NumberForm[chanSpacing,{4,1}]," m"];
Print["Gauge Length: ",NumberForm[gaugeLen,{4,1}]," m"];

Print["File Start Time: ",DateString[FromUnixTime[sampTime/10^9,TimeZone->"GMT"]],".",Mod[sampTime,10^9]," GMT"];
]; (* end If *)
packetSizeBytes=(hdrSize+(numChans*numSamps*4));
totalPackets=Round[fileSizeBytes/packetSizeBytes];
totalSamps=totalPackets*numSamps;
If[OptionValue[verbose],
Print["Total Samples Per Channel: ",totalSamps];
Print["Total Packets: ",totalPackets];
Print["Bytes per Packet: ",packetSizeBytes];
]; (* end If *)

(* Convert start and end times to UnixTime seconds *)
startTimeSecs=Block[{$TimeZone=0},UnixTime[OptionValue[startTime]]];
endTimeSecs=Block[{$TimeZone=0},UnixTime[OptionValue[endTime]]];

(* Read and extract data samples (skipping over header blocks) *)
SetStreamPosition[inStream,0];(* Rewind file *)
For[packet=1,packet<=totalPackets,packet++,
(* skip to sample time of current packet *)
SetStreamPosition[inStream,StreamPosition[inStream]+sampTimeOffset];
sampTime=BinaryRead[inStream,"UnsignedInteger64"]; (* read sample time *)
(* skip to end of header *)
SetStreamPosition[inStream,StreamPosition[inStream]+hdrSize-sampTimeOffset-8];
(* read samples if they're within requested time window, and keep only requested channels *)
If[startTimeSecs<(sampTime/10^9)<endTimeSecs,
rdBuf=Partition[BinaryReadList[inStream,"Real32",numChans*numSamps],numChans][[All,OptionValue[firstChan];;OptionValue[lastChan]]];
AppendTo[data,rdBuf];
(* if this is first saved packet, keep the time *)
If[firstPacket,firstSampTime=sampTime;firstPacket=False],
(* if packet not within requested time window, skip to start of next packet *)
SetStreamPosition[inStream,StreamPosition[inStream]+(numChans*numSamps*4)] 
]; (* end If *)
If[(sampTime/10^9)>endTimeSecs,Break[]]; (* stop if all requested samples have been read *)
]; (* end For *)
Close[inStream];

(* Data are now in sublists for each packet and we must rearrange to only sublists of requested channels. If no data at all were extracted from file, then just print message and quit *)
If[!firstPacket, (* if data in time range were found *)
If[OptionValue[firstChan]!=OptionValue[lastChan], (* if more than one channel selected *)
data=Flatten[data,1],
data=Flatten[data,2] (* else only one channel selected *)
]; (* end If *)
If[OptionValue[verbose],
If[OptionValue[firstChan]!=OptionValue[lastChan], (* if more than one channel selected *)
Print["Output Data Channels x Samples : ",Dimensions[data][[2]]," x ",Dimensions[data][[1]]],
Print["Output Data Channels x Samples : 1 x ",Length[data]] (* else only one channel selected *)
]; (* end If *)
Print["First Output Sample Time: ",DateString[FromUnixTime[firstSampTime/10^9,TimeZone->"GMT"]],".",Mod[firstSampTime,10^9]," GMT"]
]; (* end If *),
(* else if no data in time range were found *)
Print["### No data in this file for selected time range ###"];
Abort[];
] ;(* end If *)

Return[{firstSampTime,data}]
] (* end module *)

Now produce the 8 Hz files. This assumes that the files in the given directory are listed in time-increasing order with no time gaps between files. Because the input data are segmented into individual files, filtering each file alone would result in discontinuities between files. We must use the “overlap-add” method as documented in chapter 18 of [3].

inDir="/Volumes/1\ TB/Geo-Sense\ data";
outDir="/Volumes/1\ TB/decimated8Hz";

Module[{fnames,fnum,strtTime,dasDat,fltDat,overlap=Table[0,{162},{m1}],outName},
fnames=FileNames["*.raw",inDir]; (* get all the 1 kHz data file names *)
Monitor[ (* print input file names as we process them *)
For[fnum=1,fnum<=Length[fnames],fnum++, (* for each input file... *)
{strtTime,dasDat}=sintelaReadRawV3[fnames[[fnum]]]; (* read input data *)
(* apply first-stage filter to every channel of data *)
fltDat=ListConvolve[h1,#,{1,-1},0]&/@Transpose[dasDat];
fltDat[[All,1;;m1]]+=overlap; (* add overlap data from previous segment *)
(* save last m1 samples as overlap data for next segment *)
overlap=fltDat[[All,-m1;;-1]];
(* decimate filtered data after discarding the last m1 samples *)
outDat=Downsample[#,d1]&/@fltDat[[All,1;;-(m1+1)]];
outName=FileNameJoin[outDir,ToString[fnum]<>".dat"]; (* create output file name *)
Export[outName,outDat,"Real32"]; (* save data to output file *)
] (* end For *)
,fnames[[fnum]]
] ;(* end Monitor *)
] (* end Module *)

Now produce the 1 Hz file. This assumes that the 8 Hz files have been produced and are all in the directory “outDir”. Use the overlap-add method again to process the individual files as one continuous data stream.

Module[{fnames,fnum,fltDat,overlap=Table[0,{162},{m2}],inDat,out1HzSeg},
fnames=FileNames["*.dat",outDir]; (* get all the 8 Hz data file names *)
Monitor[ (* print input file names as we process them *)
For[fnum=1,fnum<=Length[fnames],fnum++, (* for each input file... *)
inDat=Import[FileNameJoin[outDir,ToString[fnum]<>".dat"],"Real32"]; (* read an 8 Hz data file *)
inDat=Partition[inDat,Length[inDat]/162]; (* partition data into 162 channels *)
(*Print["inDat "<>ToString[Dimensions@inDat]];*)
(* apply second-stage filter to every channel of data in this segment *)
fltDat=ListConvolve[h2,#,{1,-1},0]&/@inDat;
fltDat[[All,1;;m2]]+=overlap; (* add overlap data from previous segment *)
(* save last m2 samples as overlap data for next segment *)
overlap=fltDat[[All,-m2;;-1]];
(* decimate filtered data segment after discarding the last m2 samples *)
out1HzSeg=Downsample[#,d2]&/@fltDat[[All,1;;-(m2+1)]];
(*out1HzSeg=fltDat[[All,1;;-(m2+1)]];*)
If[fnum==1, (* if first segment of data *)
out1HzDat=out1HzSeg, (* make output array equal to first segment *)
(* otherwise, append the filtered segment to the output array *)
out1HzDat=MapIndexed[FlattenAt[Append[#1,out1HzSeg[[First[#2]]]],-1]&,out1HzDat];
]; (* end If *)
(*Print["out1HzDat "<>ToString[Dimensions@out1HzDat]];*)
] (* end For *)
,FileNameJoin[outDir,ToString[fnum]<>".dat"]
] ;(* end Monitor *)
(* save data to output file *)
Export[FileNameJoin[outDir,"decimated1Hz","out1HzDat.dat"],out1HzDat,"Real32"]; 
] (* end Module *)

filenames=FileNames["*.raw",inDir];

Length@filenames

filenames[[1000;;1029]]

Dimensions@out1HzDat

162 61670 4

ListPlot[out1HzDat[[6,1;;771]],PlotRange->{All,All},Joined->True]

ListPlot[out1HzDat[[6,All]],PlotRange->{All,All},Joined->True]

Manipulate[ListPlot[out1HzDat[[chan,1;;6999]],PlotRange->{All,All},Joined->True,ImageSize->Full],{chan,1,162,1}]

dat8HzSeg1=Import["/Volumes/1\ TB/testOut/1.dat","Real32"];
dat8HzSeg1=Partition[dat8HzSeg1,Length[dat8HzSeg1]/162];

ListPlot[dat8HzSeg1[[6,1;;8 770]],Joined->True]

dat8HzSeg2=Import["/Volumes/1\ TB/testOut/2.dat","Real32"];
dat8HzSeg2=Partition[dat8HzSeg2,Length[dat8HzSeg2]/162];

ListPlot[dat8HzSeg2[[6,1;;8 770]],Joined->True]

out8HzDat=MapIndexed[FlattenAt[Append[#1,dat8HzSeg2[[First[#2]]]],-1]&,dat8HzSeg1];

ListPlot[out8HzDat[[6,All]],Joined->True]

fltDat8Hz=ListConvolve[h2,#,{1,-1},0]&/@dat8Hz;

Manipulate[ListPlot[fltDat8Hz[[chan,1;;-(m2+1)]],PlotRange->{All,All},Joined->True,ImageSize->Full],{chan,1,162,1}]

fltDat1Hz=Downsample[#,d2]&/@fltDat8Hz[[All,1;;-(m2+1)]];

Dimensions@fltDat8Hz

Dimensions@fltDat1Hz

(7128-m2)/8//N

Manipulate[ListPlot[fltDat1Hz[[chan,All]],PlotRange->{All,All},Joined->True,ImageSize->Full],{chan,1,162,1}]

Length@h2

ListPlot[h1]

Manipulate[ListPlot[dat8Hz[[chan,All]],Joined->True,ImageSize->Full],{chan,1,162,1}]

ListPlot[outDat[[6,1;;100]],PlotRange->All,Joined->True]

{strtTime,ch6at1kHz}=sintelaReadRawV3["/Volumes/1\ TB/testIn/decimator_2025-01-28_19.40.42_UTC_000100.raw",firstChan->6,lastChan->6,verbose->True];

Dimensions@ch6at1kHz

ListPlot[ch6at1kHz[[All]],GridLines->Automatic]

ListPlot[ch6at1kHz[[200000;;208000]],GridLines->Automatic]

{strtTime,dat1kHz}=sintelaReadRawV3["/Volumes/1\ TB/testIn/decimator_2025-01-28_19.40.42_UTC_000100.raw",verbose->True];

ListPlot[dat1kHz[[All,6]],PlotRange->{All,All},Joined->True]

202000/125//N

Implementation notes

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

Filter delays

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.

The delay due to the first-stage filter is half the filter length times the time between samples at the original sample rate, as described in [1], section 5.8.

delay1=((Length[h1]-1)/2) (1/sr1)//N

The delay due to the second-stage filter is half the filter length times the time between samples after the first-stage decimation.

delay2=((Length[h2]-1)/2) (1/sr2)//N

The total delay is the sum of the first- and second-stage delays.

totalDelay=delay1+delay2

This delay represents a certain number of samples at the final decimated sample rate (2000 Hz in the examples above). The final output should be shifted forward in time this many samples to adjust for filter delays.

Round[totalDelay (sr2/d2)]

End effects during block processing

As with any digital filtering operation, the effects at the beginning and end of each file (i.e. block) 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.

Filter efficiency

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. (2011). Understanding digital signal processing (3rd 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