Radio-Sky Spectrograph DPS Data File Manupulator
Dave Typinski, AJ4CO
May, 2013

DPS File Header & Raw Data

(* Open one file for input *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130513\\2215.sps";
infile=OpenRead[filename,BinaryFormat->True];

(* Define read string function *)
(* file = input file stream *)
(* len = number of characters to retrieve *)
ReadString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Read fixed length portion of file header *)
chartversion="";
chartversion=ReadString[infile,10];
chartstart="";
chartstart=BinaryRead[infile,"Real64"];
chartend="";
chartend=BinaryRead[infile,"Real64"];
chartlat="";
chartlat=BinaryRead[infile,"Real64"];
chartlong="";
chartlong=BinaryRead[infile,"Real64"];
chartmax="";
chartmax=BinaryRead[infile,"Real64"];
chartmin="";
chartmin=BinaryRead[infile,"Real64"];
chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
chartsource="";
chartsource=ReadString[infile,10];
chartauthor="";
chartauthor=ReadString[infile,20];
chartlocalname="";
chartlocalname=ReadString[infile,20];
chartlocation="";
chartlocation=ReadString[infile,40];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
chartnote="";
chartnote=ReadString[infile,chartnotelength];

chartduration=Round[(chartend-chartstart)*86400000]*0.001;

chartstartdate=DateList[(chartstart-2)*86400];
chartenddate=DateList[(chartend-2)*86400];

chartstartF=DateString[chartstartdate, "Day"]<>"-"<>DateString[chartstartdate, "MonthName"]<>"-"<>DateString[chartstartdate, "Year"]<>"  "<>DateString[chartstartdate, "Hour24"]<>":"<>DateString[chartstartdate, "Minute"]<>":"<>DateString[chartstartdate, "Second"]<>"."<>DateString[chartstartdate, "Millisecond"]<>" UTC";

chartendF=DateString[chartenddate, "Day"]<>"-"<>DateString[chartenddate, "MonthName"]<>"-"<>DateString[chartenddate, "Year"]<>"  "<>DateString[chartenddate, "Hour24"]<>":"<>DateString[chartenddate, "Minute"]<>":"<>DateString[chartenddate, "Second"]<>"."<>DateString[chartenddate, "Millisecond"]<>" UTC";

chartversionF=ToString[FromDigits[StringTake[chartversion,5]]]<>"."<>ToString[FromDigits[StringTake[chartversion,{6,7}]]]<>"."<>ToString[FromDigits[StringTake[chartversion,{8,10}]]];

sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* create temp table and read in spectro data *)
numrecs=sweeps;
indata=Table[0,{i,1,numrecs},{j,1,chartchannels*2+1}];
i=1;

Do[{
indata[[i,j]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,numrecs}
,{j,1,chartchannels*2+1}
];
Close[infile];

(* Find actual number of whole records present in data file *)
(* Check for EOR markers *)
(* The EOR marker is two consecutive 0xFE characters, ASCII "þ" *)
(* This is seens as 65278 decomal when read as a 2-byte integer  *)
chartrecs=Count[indata,65278,2];
chartsamplerate=chartrecs/chartduration;

chartmaxF=Max[indata[[1;;chartrecs,1;;chartchannels*2]]];
chartminF=Min[indata[[1;;chartrecs,1;;chartchannels*2]]];

freqlo=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"LOWF"][[1,2]]+1 ,StringPosition[chartnote,"LOWF"][[1,2]]+8}]]*10^-6;
freqhi=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"HIF"][[1,2]]+1 ,StringPosition[chartnote,"HIF"][[1,2]]+8}]]*10^-6;

channelsep=(freqhi-freqlo)/chartchannels*10^3;

channeldwelltime=Round[chartduration/(chartrecs*chartchannels)*10^6];

(* Display file header data to screen *)
Print["          File path: ", filename]
Print["        RSS version: ",chartversionF]
Print["        Chart start: ",chartstartF]
Print["          Chart end: ",chartendF]
Print[]
Print["         Local name: ",chartlocalname]
Print["     Chart location: ",chartlocation]
Print["           Latitude: ",ToString[PaddedForm[chartlat,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(S) -","(N) +"}]]<>" degrees"]
Print["          Longitude: ",ToString[PaddedForm[chartlong,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(W) -","(E) +"}]]<>" degrees"]
Print["          Time zone: UTC+",chartTZ]
Print["           Observer: ",chartauthor]
Print[]
Print["        Data Source: ",chartsource]
Print["      Data duration: ",ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]," s"]
Print["        Sweep Start: ",ToString[PaddedForm[freqhi,{10,6},NumberPadding->{"","0"}]]," MHz"]
Print["          Sweep End: ",ToString[PaddedForm[freqlo,{10,6},NumberPadding->{"","0"}]]," MHz"]
Print["             Sweeps: ",chartrecs]
Print["         Sweep rate: ",ToString[PaddedForm[chartsamplerate,{6,3},NumberPadding->{"","0"}]]," Hz"]
Print["           Channels: ",chartchannels]
Print[" Channel separation: ",ToString[PaddedForm[channelsep,{6,3},NumberPadding->{"","0"}]]," kHz"]
Print[" Channel dwell time: ",ToString[channeldwelltime]<>" μs"]
Print[]
Print["           Data Max: ",ToString[chartmaxF]<>" mV"]
Print["           Data Min: ",ToString[chartminF]<>" mV"]
Print[]
Print["  Chart note length: ",chartnotelength]
Print["       Chart note 1: ",chartnote]

(* create 12-bit data from input file array *)
indata12bit=indata[[1;;chartrecs,1;;chartchannels*2]];

(* verify record count *)
dataA=Table[0.0,{i,1,chartrecs},{j,1,chartchannels}];
dataB=Table[0.0,{i,1,chartrecs},{j,1,chartchannels}];

Do[
Do[
{dataA[[i,j]]=indata12bit[[i,j*2]],
dataB[[i,j]]=indata12bit[[i,j*2-1]]}
,{j,1,chartchannels}
],{i,1,chartrecs}
];

dataAM=N[Mean[dataA]];
dataBM=N[Mean[dataB]];

Print[]
Print["      Continuum display (mean values of ",ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]," seconds of data)"]

Needs["PlotLegends`"]
ListLinePlot[{dataAM,dataBM},PlotRange->{1600,2000},PlotStyle->{Darker[Green],Red},PlotLegend->{"RCP","LCP"},LegendPosition->{1.1,-0.4},LegendSize->0.3,LegendShadow->None]

dataAM//MatrixForm