Reading Sintela DAS V3 Raw Data Files

P. McGill
10Mar25

Write a function to read  Sintela DAS V3 format files according to the documentation at https://sintela.atlassian.net/wiki/spaces/DEMO/pages/336917/Phase+Data+Sensor+Data+Sintela+Binary+Format.

sintelaReadRawV3[fileName_String]:= Module[ (* version 1.0.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,header,rdBuf,data={}},

(* 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[
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 *)

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 *)
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["Start Time: ",DateString[FromUnixTime[sampTime/10^9,TimeZone->"GMT"]],".",Mod[sampTime,10^9]," GMT"];

packetSizeBytes=(hdrSize+(numChans*numSamps*4));
totalPackets=Round[fileSizeBytes/packetSizeBytes];
totalSamps=totalPackets*numSamps;
Print["Total Samples Per Channel: ",totalSamps];
Print["Total Packets: ",totalPackets];
Print["Bytes per Packet: ",packetSizeBytes];

(* Read and extract data samples (skipping over header blocks) *)
SetStreamPosition[inStream,0];(* Rewind file *)
Monitor[ (* show a progress indicator as we read the packets *)
For[packet=1,packet<=totalPackets,packet++,
(* skip over header and read samples *)
SetStreamPosition[inStream,StreamPosition[inStream]+hdrSize];
rdBuf=BinaryReadList[inStream,"Real32",numChans*numSamps];
AppendTo[data,rdBuf];
] (* end For *),
ProgressIndicator[packet,{1,totalPackets}]
] (* end Monitor *);

(* data are now in sublists of packetSizeBytes length, and we must rearrange to sublists of numChans length *);
data=Partition[Flatten[data],numChans];

Close[inStream];
Return[data]
] (* end module *)

Read a random file from the Jan 2025 deployment.

(* filename="/Volumes/ProjectLibrary/902401_GeoSense/Project.Data/Deployments/DAS_2025_01_28/Geo-Sense\ data/decimator_2025-02-03_08.38.39_UTC_000721.raw"; *)

filename="/Volumes/Geo-Sense/2025\ Jan\ Data/decimator_2025-01-29_05.31.40_UTC_000146.raw";

dasDat=sintelaReadRawV3[filename];

Dimensions@dasDat

Check the first second of data from the first channel to see if it looks reasonable.

ListPlot[dasDat[[1;;1000,1]]]

Make a data browser to view the first 10 seconds of data for any channel.

MinMax@dasDat

Manipulate[ListPlot[dasDat[[1;;10000,chan]],PlotRange->{-50,250},ImageSize->Large],{chan,1,162,1}]

Plot a phasogram for the first 2 minutes of data by finding the RMS value for each 100 ms segment of time.

ArrayPlot[Transpose[BlockMap[RootMeanSquare[#]&,#,100]&/@Transpose[dasDat[[1;;120000,All]]]],Frame->True,FrameLabel->{"Time (s)","DAS Channel"},FrameTicks->Automatic,PlotRange->{-10,20},ColorFunction->"BlueGreenYellow",ImageSize->Full,AspectRatio->0.7]

The high DC offset of many of the channels increases the dynamic range of the signal and makes it difficult to see small variations in the signals. Improve this by subtracting the mean for each 100 ms segment before taking the RMS value.

ArrayPlot[Transpose[BlockMap[RootMeanSquare[#-Mean[#]]&,#,100]&/@Transpose[dasDat[[1;;120000,All]]]],Frame->True,FrameLabel->{"Time (s)","DAS Channel"},FrameTicks->Automatic,PlotRange->{0,0.5},ColorFunction->"BlueGreenYellow",ImageSize->Full,AspectRatio->0.7]

Instead of removing the mean of each time segment, apply a 2 Hz highpass filter to each channel.

dasDatHp=Transpose[HighpassFilter[#,(Quantity[2, "Hz", hertz, "Hertz"]),SampleRate->1000]&/@Transpose[dasDat]];

ByteCount@dasDat

ByteCount@dasDatHp

Dimensions@dasDat

Dimensions@dasDatHp

ArrayPlot[Transpose[BlockMap[RootMeanSquare[#]&,#,100]&/@Transpose[dasDatHp[[1;;120000,All]]]],Frame->True,FrameLabel->{"Time (s)","DAS Channel"},FrameTicks->Automatic,PlotRange->{-10,20},ColorFunction->"BlueGreenYellow",ImageSize->Full,AspectRatio->0.7]

ListPlot[{dasDat[[1;;10000,126]],HighpassFilter[dasDat[[1;;10000,126]],Quantity[2,"Hertz"],SampleRate->1000]},PlotRange->All]

Periodogram[{dasDat[[1;;10000,126]],HighpassFilter[dasDat[[1;;10000,126]],Quantity[2,"Hertz"],1001,SampleRate->1000]},SampleRate->1000,Frame->True,GridLines->Automatic,PlotRange->{{0,10},All}]

Periodogram[{dasDat[[1;;10000,126]],dasDatHp[[1;;10000,126]]},SampleRate->1000,Frame->True,GridLines->Automatic,PlotRange->{{0,10},All},ImageSize->Large]

Manipulate[ListPlot[dasDatHp[[1;;10000,chan]],Frame->True,PlotRange->{-10,25},ImageSize->Large],{chan,1,162,1}]