MARS DAS Channel Browser

P. McGill
9May25

We need a way to look at the data from the MARS Sintela Onyx DAS so we can compare events to the Geo-Sense Sintela Nano DAS.

Import the function to read Sintela  DAS files from the notebook “Sintela DAS Raw Data v3 File Reader v1.2.0.nb”.

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

Import MARS DAS data, stored locally on an SD card.

filename="/Volumes/Geo-Sense/MARS\ Data/2024_12_14/onyxAcquisition_2024-12-14_16.43.51_UTC_004401.raw";

(*filename="/Volumes/Geo-Sense/test_data/decimator_2025-01-28_23.31.57_UTC_000118.raw";*)

{time,dasDat}=sintelaReadRawV3[filename,lastChan->150,verbose->True];

Make a data browser that shows the first 10 seconds of data for any channel.

Manipulate[ListPlot[dasDat[[1;;10000,chan]],PlotRange->{-500,500},Frame->True,ImageSize->Large],{chan,1,150,1}]

The signals are much larger than the Geo-Sense DAS data, probably due to wave energy close to shore. Any signal from a gravity-flow event will be a tiny AC variation riding on a huge low-frequency offset.

Plot a spectrogram for all of the channel 50 data.

Spectrogram[dasDat[[All,50]] 100,SampleRate->1000,FrameLabel->{"Time (s)","Frequency (Hz)"},ImageSize->Full]

Focus on the frequency range from 0–50 Hz.

Spectrogram[dasDat[[All,50]] 100,SampleRate->1000,PlotRange->{0,50},FrameLabel->{"Time (s)","Frequency (Hz)"},ImageSize->Full]

There seems to be some artifacts from processing. Try processing 1024-sample partitions with a 512-sample offset and a Blackman smoothing window.

Spectrogram[dasDat[[All,50]] 20,1024,512,BlackmanWindow,SampleRate->1000,PlotRange->{0,50},FrameLabel->{"Time (s)","Frequency (Hz)"},ImageSize->Full]

Try eliminating the large, low-frequency signals by applying a 2 Hz highpass filter.

Spectrogram[HighpassFilter[dasDat[[All,50]],(Quantity[2, "Hz", hertz, "Hertz"]),1023,SampleRate->1000] 5,1024,512,BlackmanWindow,SampleRate->1000,PlotRange->{0,50},FrameLabel->{"Time (s)","Frequency (Hz)"},ImageSize->Full]