Gage Digitizer Data File Utilities
Dave Typinski, AJ4CO Observatory
September, 2013

Plot one FFT from one SIG file - converts all the .sig files found in the specified directory to .wav files

(********************************************************************************************)
(*                                                                                          *)
(*  User configuration section                                                              *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Set directory for input                                      *)
filebase="R:\\Observation Records\\AJ4CO Observatory\\TWB\\Gage\\2013-09-23_04-07-10_CH01\\Folder.00001";

(********************************************************************************************)
(*                                                                                          *)
(*  End of user configuration section                                                       *)
(*                                                                                          *)
(********************************************************************************************)

(* Finda and count all Gage .sig files *)
fileext=".sig";
SetDirectory[filebase];
infiles=FileNames["*"<>fileext];
filecount=infiles//Length;

(* 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
];

(* Start of Do loop to cycle thru all files would go here *)
infile=OpenRead[infiles[[1]],BinaryFormat->True];

(* Read fixed length file header - 512 bytes total *)
chartversion="";
chartversion=ReadString[infile,14];

dummy=BinaryRead[infile,"Integer16"];
dummy=ReadString[infile,9];
dummy=BinaryRead[infile,"Integer16"];
dummy=ReadString[infile,256];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];

chartsamplerateindex="";
chartsamplerateindex=BinaryRead[infile,"Integer16"];
Print["Sample rate index = ",chartsamplerateindex," -- Note: index 22 = 10 MHz sample rate"]

chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];

dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];

chartsamples="";
chartsamples=BinaryRead[infile,"Integer32"];
Print["Chart sample count = ",chartsamples]

chartinputrangeindex="";
chartinputrangeindex=BinaryRead[infile,"Integer16",ByteOrdering->-1];

chartinputcoupling="";
chartinputcoupling=BinaryRead[infile,"Integer16"];

dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];

charttime="";
charttime=BinaryRead[infile,"UnsignedInteger16"];
Print["Chart time = ",charttime]

chartdate="";
chartdate=BinaryRead[infile,"UnsignedInteger16"];
Print["Chart date = ",chartdate]

dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];

chartboardtype="";
chartboardtype=BinaryRead[infile,"UnsignedInteger16"];

chartresolution="";
chartresolution=BinaryRead[infile,"Integer16"];

dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];

chartsampleoffset="";
chartsampleoffset=BinaryRead[infile,"Integer16"];

chartsampleresolution="";
chartsampleresolution=BinaryRead[infile,"Integer16"];

chartsamplebits="";
chartsamplebits=BinaryRead[infile,"Integer16"];
Print["Chart sample resolution = ",chartsamplebits," bits"]

chartextendedtriggertime="";
chartextendedtriggertime=BinaryRead[infile,"UnsignedInteger32"];

chartinputimedance="";
chartinputimedance=BinaryRead[infile,"Integer16"];

dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Real32"];
dummy=BinaryRead[infile,"Real32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"UnsignedInteger16"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=ReadString[infile,127];
(* End of header data *)



signal=Table[0.,{i,1,chartsamples}];

Do[
signal[[i]]=BinaryRead[infile,"Integer16"]
,{i,1,chartsamples}
];

Print["2048 samples, raw data from Gage"]
ListLinePlot[signal[[1;;2048]]]

fsignal=Abs[Fourier[signal[[1;;2048]]]]^2/50; (* power spectrum in uW into 50 ohms *)

fsignal = fsignal+fsignal[[1]]; (* add DC offset to spectrum *)

fsignal=fsignal/1000; (* convert uW to mW *)

fsignal=10*Log10[fsignal/1]; (* convert mW to dBm *)

Print[]
Print["FFT of 2 MHz slice of baseband output from 0 to 5 MHz in dBm, RBW = 4.88 kHz"]
ListLinePlot[fsignal[[1;;1025]],PlotRange->{-2,0}] (* plot 4.8 to 2.8 MHz slice *)

Dimensions[signal]

SIG to SPS converter - converts all the .sig files in one directory to a single long SPS file

(********************************************************************************************)
(*                                                                                          *)
(*  User configuration section                                                              *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

Print["Processing started "<>DateString[]]

(* Set UTC manual offset (in hours) to correct for file creation times on computers with different clock TZ settings *)
utcoffset=3.;

(* Set directory for input -- contains one or more .sig Gage signal files *)
filebase="R:\\Observation Records\\AJ4CO Observatory\\TWB\\Gage\\2013-11-16_22-31-14_CH01\\Folder.00001";

(* Set Gage data sample rate in MHz *)
samplerate=10;

(* Set measured cycle time in milliseconds between trigger events *)
cycletime=375;

(* Set number of FFT bins *)
fftbins=2048;

(* Set low freq and high freq limits of FFT display in MHz *)
(* This is the portion of the FFT that we want to see; i.e., the actual near-baseband output of the TWB *)
(* See below for the RF input frequency setting to be placed in the SPS file *)
displaylof=2.8;
displayhif=4.8;

(* Set the R8500 tuned frequency used during digitization, MHz *)
r8500=73.5;

(* SPS file detector sensitivity in ADC counts per dB *)
adcsens=50;

(* Set DC offset (in μW) to be added to FFT before calculatind dBm *)
dcoffset=100;

(* Set DC offset (in dBm) to be added to FFT before conversion to SPS data *)
postdcoffset=11;

(* Set DC offset (in ADC counts) to be added before export to SPS file *)
spsoffset=1000;

(* Set desired sweep rate in SPS data file in Hz *)
SPSsweeprate=10;

(********************************************************************************************)
(*                                                                                          *)
(*  End of user configuration section                                                       *)
(*                                                                                          *)
(********************************************************************************************)



(********************************************************************************************)
(*                                                                                          *)
(*  Function definitions section                                                            *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Define proper rounding function as MM only rounds 0.5 up when the result is an even number *)
SetAttributes[Round2,Listable];
Round2[a_]:=If[FractionalPart[a]≥1/2,Ceiling[a],Floor[a]];

(* 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
];

(********************************************************************************************)
(*                                                                                          *)
(*  End of function definitions section                                                     *)
(*                                                                                          *)
(********************************************************************************************)



(********************************************************************************************)
(*                                                                                          *)
(*  Main program code                                                                       *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Set the SPS file RF frequency range labels in MHz (derived from R8500 tuned frequency) *)
(* TWB 1st LO (in the upconverter) is 97 MHz, TWB output is 2 MHz wide *)
SPSlof=97-r8500-1.;
SPShif=97-r8500+1.;

(* calculate some FFT and SPS file parameters *)
(* find RBW *)
fftrbw=(samplerate*10^6)/fftbins;
(* fined low frequency element of FFT to export to SPS file *)
lofelement=Round2[fftbins*displaylof/samplerate]+1;
(* find high frequency element of FFT to export to SPS file *)
hifelement=Round2[fftbins*displayhif/samplerate]+1;
(* find total number of bins to export to SPS file (this is the SPS file channel count) *)
outputbincount=hifelement-lofelement+1;

(* Find and count all Gage .sig files in the specified directory *)
fileext=".sig";
SetDirectory[filebase];
infiles=FileNames["*"<>fileext];
filecount=infiles//Length;
Print["Processing ",ToString[filecount]," files..."]

(* Open Gage data file for read *)
(* Open SPS file for write *)
infile=OpenRead[infiles[[1]],BinaryFormat->True];

(* Read fixed length header portion of Gage sig file - 512 bytes total *)
chartversion="";
chartversion=ReadString[infile,14];
dummy=BinaryRead[infile,"Integer16"];
dummy=ReadString[infile,9];
dummy=BinaryRead[infile,"Integer16"];
dummy=ReadString[infile,256];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartsamplerateindex="";
chartsamplerateindex=BinaryRead[infile,"Integer16"];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartsamples="";
chartsamples=BinaryRead[infile,"Integer32"];
chartinputrangeindex="";
chartinputrangeindex=BinaryRead[infile,"Integer16",ByteOrdering->-1];
chartinputcoupling="";
chartinputcoupling=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
charttime="";
charttime=BinaryRead[infile,"UnsignedInteger16"];
chartdate="";
chartdate=BinaryRead[infile,"UnsignedInteger16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartboardtype="";
chartboardtype=BinaryRead[infile,"UnsignedInteger16"];
chartresolution="";
chartresolution=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartsampleoffset="";
chartsampleoffset=BinaryRead[infile,"Integer16"];
chartsampleresolution="";
chartsampleresolution=BinaryRead[infile,"Integer16"];
chartsamplebits="";
chartsamplebits=BinaryRead[infile,"Integer16"];
chartextendedtriggertime="";
chartextendedtriggertime=BinaryRead[infile,"UnsignedInteger32"];
chartinputimedance="";
chartinputimedance=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Real32"];
dummy=BinaryRead[infile,"Real32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"UnsignedInteger16"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=ReadString[infile,127];
(* End of header data *)

(* Calculate observation duration in seconds *)
observationduration=cycletime*filecount/1000.;

(* Calculate time length of digitized data in milliseconds *)
dburstduration=chartsamples/samplerate/1000.;

(* Calculate number of sweeps in the Gage data *)
fftsweeps=Floor[chartsamples/fftbins];

(* Calculate padding in milliseconds needed before start of next digitized burst of data *)
SPSpaddingtime=cycletime-dburstduration;
SPSpaddingsweeps=Round2[fftsweeps*cycletime/dburstduration-fftsweeps];
SPSsweepsperburst=Round2[fftsweeps*cycletime/dburstduration];
SPSsweepsperobservation=SPSsweepsperburst*filecount;

(* Calculate start and stop time to be placed in SPS header *)
SPSchartstart=DateDifference[{1900,1,1},FileDate[infiles[[1]],"Creation"]]+2; (* Time since Dec 30, 1899, 00:00:00.000 UTC *)
SPSchartstart=SPSchartstart+utcoffset/24.; (* apply manual UTC offset correction *)
SPSchartstartf=DateString[(SPSchartstart-2)*86400,{"Day"," ","MonthNameShort"," ","Year"," ","Hour24",":","Minute",":","Second"," UTC"}];
SPSchartend=SPSchartstart+SPSsweepsperobservation/SPSsweeprate/86400.;
SPSchartendf=DateString[(SPSchartend-2)*86400,{"Day"," ","MonthNameShort"," ","Year"," ","Hour24",":","Minute",":","Second"," UTC"}];
SPSfileduration=SPSsweepsperobservation/SPSsweeprate/86400.*24.;
SPStimeexpansionratio=SPSfileduration*3600./observationduration;
SPSminute=60./SPStimeexpansionratio;

(* Generate SPS & PDF filenames and open SPS output file for write *)
SPSoutfilename="AJ4CO-TWB-"<>ToString[DateList[(SPSchartstart-2)*86400][[1]]]<>ToString[DateList[(SPSchartstart-2)*86400][[2]]]<>ToString[DateList[(SPSchartstart-2)*86400][[3]]]<>ToString[DateList[(SPSchartstart-2)*86400][[4]]]<>ToString[DateList[(SPSchartstart-2)*86400][[5]]]<>ToString[Round2[DateList[(SPSchartstart-2)*86400][[6]]]];
PDFoutfilename=SPSoutfilename<>" Analysis";
outfile=OpenWrite[SPSoutfilename<>".sps",BinaryFormat->True];

(* Write SPS file header *)
(* Set RSS data file version number 2.2.19 *)
SPSchartversion="0000202019";
WriteString[outfile,SPSchartversion];
BinaryWrite[outfile,SPSchartstart,"Real64"];
BinaryWrite[outfile,SPSchartend,"Real64"];
SPSchartlat=0;
BinaryWrite[outfile,SPSchartlat,"Real64"];
SPSchartlong=0;
BinaryWrite[outfile,SPSchartlong,"Real64"];
SPSchartmax=4095;
BinaryWrite[outfile,SPSchartmax,"Real64"];
SPSchartmin=0;
BinaryWrite[outfile,SPSchartmin,"Real64"];
SPSchartTZ=0;
BinaryWrite[outfile,SPSchartTZ,"Integer16"];
SPSchartsource="TWB       ";
WriteString[outfile,SPSchartsource];
SPSchartauthor="Dave Typinski       ";
WriteString[outfile,SPSchartauthor];
SPSchartlocalname="AJ4CO_TWB           ";
WriteString[outfile,SPSchartlocalname];
SPSchartlocation="High Springs, Florida                   ";
WriteString[outfile,SPSchartlocation];
SPSchartchannels=outputbincount;
BinaryWrite[outfile,SPSchartchannels,"Integer16"];
SPSchartnote="*[[*SWEEPS"<>ToString[SPSsweepsperobservation]<>"ÿLOWF"<>ToString[Round[SPSlof*10^6]]<>"ÿHIF"<>ToString[Round[SPShif*10^6]]<>"ÿSTEPS411ÿDUALSPECFILEFalseÿRCVRÿCOLORRES1ÿ*]]*";
SPSchartnotelength=StringLength[SPSchartnote];
BinaryWrite[outfile,SPSchartnotelength,"Integer32"];
WriteString[outfile,SPSchartnote];

(* Generate text report entries *)
txtoutput=Table["",{i,1,60}];
txtoutput[[1]]=Style["TWB Gage Signal File to RSS SPS File Converter",Larger,Darker[Darker[Blue]],Bold];
txtoutput[[2]]="Analysis run date: "<>DateString[{"Day"," ","MonthNameShort"," ","Year"," ","Hour24",":","Minute",":","Second"," Local"}];

txtoutput[[10]]=Style["Data Conversion Analysis Report",Bold];

txtoutput[[12]]="Observation start time: "<>SPSchartstartf;
txtoutput[[13]]="Duration of observation: "<>ToString[observationduration]<>" real-time seconds";

txtoutput[[15]]="Data directory: "<>filebase;
txtoutput[[16]]="Number of digitized input files: "<> ToString[filecount];
txtoutput[[17]]="First input filename: "<>infiles[[1]];
txtoutput[[18]]="Last input filename: "<>infiles[[filecount]];

txtoutput[[20]]="Digitized burst file size: "<>ToString[chartsamples]<>" samples per file";
txtoutput[[21]]="Digitized burst file sample rate: "<>ToString[samplerate]<>" MHz";
txtoutput[[22]]="Digitized burst file duration: "<>ToString[dburstduration]<>" ms";
txtoutput[[23]]="Digitized burst cycle time: "<>ToString[cycletime]<>" ms";
txtoutput[[24]]="Dead time between data bursts: "<>ToString[SPSpaddingtime]<>" ms";
txtoutput[[25]]="Digitization coverage: "<>ToString[dburstduration/cycletime*100.]<>" percent";

txtoutput[[27]]="FFT bins: "<>ToString[fftbins];
txtoutput[[28]]="FFT sweep time: "<>ToString[fftbins/samplerate*1.]<>" μs";
txtoutput[[29]]="FFT sweeps per digitized data burst: "<>ToString[fftsweeps];
txtoutput[[30]]="Dead FFT sweeps between each digitized data burst: "<>ToString[SPSpaddingsweeps];
txtoutput[[31]]="FFT sweeps per digitized data burst including dead time padding: "<>ToString[SPSsweepsperburst];
txtoutput[[32]]="Total FFT sweeps for "<>ToString[filecount]<>" input files, including padding: "<>ToString[SPSsweepsperobservation];

txtoutput[[34]]="FFT BW: "<>ToString[fftrbw*(fftbins/2)/10^6]<>" MHz";
txtoutput[[35]]="FFT RBW: "<>ToString[fftrbw/1000.]<>" kHz";
txtoutput[[36]]="FFT Windowing: "<>"None (uniform window)";
txtoutput[[37]]="FFT display low frequency: "<>ToString[displaylof]<>" MHz ( FFT bin # "<>ToString[lofelement]<>" )";
txtoutput[[38]]="FFT display high frequency: "<>ToString[displayhif]<>" MHz ( FFT bin # "<>ToString[hifelement]<>" )";
txtoutput[[39]]="Total FFT bins exported to SPS file: "<>ToString[outputbincount];

(* txtoutput[[41]]="DC offset per FFT element zero: "<>ToString[fsignal[[1]]]<>" μW"; *)
txtoutput[[42]]="DC offset applied to FFT before calculating dBm: "<>ToString[dcoffset]<>" μW";
txtoutput[[43]]="DC offset applied to FFT after calculating dBm: "<>ToString[postdcoffset]<>" dBm";
txtoutput[[44]]="SPS file detector sensitivity: "<>ToString[adcsens]<>" ADC counts per dB";
txtoutput[[45]]="DC offset applied to SPS data before export to SPS file: "<>ToString[spsoffset]<>" ADC counts";

txtoutput[[47]]="SPS output file name: "<>SPSoutfilename<>".sps";
txtoutput[[48]]="SPS data file sweep rate: "<>ToString[SPSsweeprate]<>" Hz";
txtoutput[[49]]="SPS file start time: "<>SPSchartstartf;
txtoutput[[50]]="SPS file end time: "<>SPSchartendf;
txtoutput[[51]]="SPS file duration: "<>ToString[SPSfileduration]<>" hours";
txtoutput[[52]]="SPS file time expansion ratio: "<>ToString[SPStimeexpansionratio]<>" to 1";
txtoutput[[53]]="One minute of SPS file time represents: "<>ToString[SPSminute*1000.]<>" milliseconds of real time";

(* close first .sig file to reset processing -- main loop must start with the header read operations *)
Close[infile];

(***********************************************************)
(*                                                         *)
(* MAIN LOOP to cycle through all the available .sig files *)
(*                                                         *)
(***********************************************************)

Do[
{
infile=OpenRead[infiles[[k]],BinaryFormat->True];

(* Read fixed length header portion of Gage sig file - 512 bytes total *)
chartversion="";
chartversion=ReadString[infile,14];
dummy=BinaryRead[infile,"Integer16"];
dummy=ReadString[infile,9];
dummy=BinaryRead[infile,"Integer16"];
dummy=ReadString[infile,256];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartsamplerateindex="";
chartsamplerateindex=BinaryRead[infile,"Integer16"];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartsamples="";
chartsamples=BinaryRead[infile,"Integer32"];
chartinputrangeindex="";
chartinputrangeindex=BinaryRead[infile,"Integer16",ByteOrdering->-1];
chartinputcoupling="";
chartinputcoupling=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"Integer32"];
charttime="";
charttime=BinaryRead[infile,"UnsignedInteger16"];
chartdate="";
chartdate=BinaryRead[infile,"UnsignedInteger16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartboardtype="";
chartboardtype=BinaryRead[infile,"UnsignedInteger16"];
chartresolution="";
chartresolution=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
chartsampleoffset="";
chartsampleoffset=BinaryRead[infile,"Integer16"];
chartsampleresolution="";
chartsampleresolution=BinaryRead[infile,"Integer16"];
chartsamplebits="";
chartsamplebits=BinaryRead[infile,"Integer16"];
chartextendedtriggertime="";
chartextendedtriggertime=BinaryRead[infile,"UnsignedInteger32"];
chartinputimedance="";
chartinputimedance=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Integer16"];
dummy=BinaryRead[infile,"Real32"];
dummy=BinaryRead[infile,"Real32"];
dummy=BinaryRead[infile,"Integer32"];
dummy=BinaryRead[infile,"UnsignedInteger16"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=BinaryRead[infile,"UnsignedInteger32"];
dummy=ReadString[infile,127];
(* End of header data *)

(* Read all signal data from one file *)
signal=Table[0.,{i,1,chartsamples}];
Do[
signal[[i]]=BinaryRead[infile,"Integer16"]
,{i,1,chartsamples}
];

(* Perform FFT analysis of signal data for one file *)
Do[
{
fsignal=Abs[Fourier[signal[[fftbins*j-fftbins+1;;fftbins*j]]]]^2/50; (* power spectrum in uW into 50 ohms *)
fsignal = fsignal+dcoffset; (* add DC offset to spectrum prior to conversion to dBm *)
fsignal=fsignal/1000; (* convert uW to mW *)
fsignal=10*Log10[fsignal/1]; (* convert mW to dBm *)
fsignal = fsignal+postdcoffset;  (* bring floor above 0 dBm prior to conversion to SPS data *)
spssignal=Round2[fsignal[[lofelement;;hifelement]]*adcsens+spsoffset]; (* convert to SPS data and add DC offset *)

(* Write one FFT sweep to the SPS file *)
Do[
{
BinaryWrite[outfile,spssignal[[i]],"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,outputbincount}
];
(* Write EOR marker *)
(* The EOR marker is two consecutive 0xFE characters, ASCII "þ" *)
(* This is seens as 65278 decimal when read as a 2-byte integer  *)
BinaryWrite[outfile,65278,"UnsignedInteger16",ByteOrdering->+1];

}
,{j,1,fftsweeps}
];

(* Add one set of inter-burst padding sweeps to SPS file *)
Do[
{
(* Write one FFT sweep to the SPS file *)
Do[
{
BinaryWrite[outfile,0,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,outputbincount}
];
(* Write EOR marker *)
(* The EOR marker is two consecutive 0xFE characters, ASCII "þ" *)
(* This is seens as 65278 decimal when read as a 2-byte integer  *)
BinaryWrite[outfile,65278,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,SPSpaddingsweeps}
];

Close[infile];

}
,{k,1,filecount}
];

Close[outfile];

(* Output the analysis report *)
txtoutput[[41]]="DC offset per FFT element zero: "<>ToString[fsignal[[1]]]<>" μW (last FFT sweep of last data file)"; 
txtoutput[[3]]="Analysis complete: "<>DateString[{"Day"," ","MonthNameShort"," ","Year"," ","Hour24",":","Minute",":","Second"," Local"}];

panelout=Panel[Column[{
Column[txtoutput[[10;;11]]],
Column[txtoutput[[12;;15]],Alignment->":"],
Column[txtoutput[[16;;54]],Alignment->":"]
}],
Column[txtoutput[[1;;4]]]];

Export[PDFoutfilename<>".pdf",panelout];

Print["Processing ended "<>DateString[]]

(* Create output array                                          *)
txtoutput[[13]]="        RSS version: "<>chartversionF;
txtoutput[[14]]="        Chart start: "<>chartstartF;
txtoutput[[15]]="          Chart end: "<>chartendF;

txtoutput[[17]]="         Local name: "<>chartlocalname;
txtoutput[[18]]="     Chart location: "<>chartlocation;
txtoutput[[19]]="           Latitude: "<>ToString[PaddedForm[chartlat,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(S) -","(N) +"}]]<>" degrees";
txtoutput[[20]]="          Longitude: "<>ToString[PaddedForm[chartlong,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(E) +","(W) -"}]]<>" degrees";
txtoutput[[21]]="          Time zone: UTC+"<>IntegerString[chartTZ];
txtoutput[[22]]="           Observer: "<>chartauthor;

txtoutput[[24]]="      Data duration: "<>ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]<>" s";
txtoutput[[25]]="        Sweep Start: "<>ToString[PaddedForm[freqhi,{10,6},NumberPadding->{"","0"}]]<>" MHz";
txtoutput[[26]]="          Sweep End: "<>ToString[PaddedForm[freqlo,{10,6},NumberPadding->{"","0"}]]<>" MHz";
txtoutput[[27]]="             Sweeps: "<>ToString[sweeps];
txtoutput[[28]]="         Sweep Rate: "<>ToString[PaddedForm[chartsweeprate,{6,3},NumberPadding->{"","0"}]]<> "Hz";
txtoutput[[29]]="           Channels: "<>IntegerString[chartchannels];
txtoutput[[30]]=" Channel separation: "<>ToString[PaddedForm[channelsep,{6,3},NumberPadding->{"","0"}]]<>" kHz";
txtoutput[[31]]=" Channel Dwell Time: "<>ToString[channeldwelltime]<>" μs";

txtoutput[[33]]="           Data Max: "<>ToString[chartmaxF]<>" ADC Counts (raw data before averaging)";
txtoutput[[34]]="           Data Min: "<>ToString[chartminF]<>" ADC Counts (raw data before averaging)";

txtoutput[[36]]=Style["Continuum Plot: Average of "<>ToString[sweeps]<>" Sweeps, "<>IntegerString[3*(k-1)]<>" dB Attenuation"<>" = "<>IntegerString[Round2[srctemps[[k]]]]<>" kK at "<>calplane,Bold];

txtoutput[[38]]=Rasterize[plots[[k]],ImageSize->650];

(* Create channel frequency list                                *)
dataFM=Table[ToString[PaddedForm[N[freqhi-i*channelsep*10^-3],{5,3},NumberPadding->{"","0"}]],{i,1,chartchannels}];

(* Prepend column labeles and position numbers                  *)
A=Prepend[avgA," ADC (counts) "];
Frqs=Prepend[dataFM,"  Freq (MHz)  "];
Ch=Prepend[Table[i,{i,1,chartchannels}],"Channel"];
Dummy=Table["        ",{i,1,chartchannels+1}];

txtoutput[[40]]=TableForm[Transpose[{Dummy,Ch,Frqs,A}],TableAlignments->Center];

(* Output the result                                            *)
panels[[k]]=Panel[Column[{
Column[txtoutput[[10;;11]]],
Column[txtoutput[[12;;37]],Alignment->":"],
txtoutput[[38]],
"",
Style["    ADC Output Averaged Over "<>IntegerString[sweeps]<>" Sweeps",Black],
txtoutput[[40]]
}],
Column[txtoutput[[1;;3]]]];

Export[rptfilename<>IntegerString[(k-1)*stepsize,10,2]<>" dB Attenuation Step.pdf",panels[[k]]];

DateString[(SPSchartstart-2)*86400,{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute",":","Second"," Local"}]

Close[outfile]

outfile

ToString[Round2[DateList[(SPSchartstart-2)*86400][[6]]]]

txtoutput[[34]]="FFT BW: "<>ToString[fftrbw*(fftbins/2)/10^6]<>" MHz"

Close

DPS Channel Rotator - Rotates channel 200 back down below channel 199 and writes new SPS file

(********************************************************************************************)
(*                                                                                          *)
(*  RSS and DPS data have channel 200 rolled over above channel 1.                          *)
(*  This routine puts channel 200 data back where it belongs.                               *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one file for input                                      *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130503\\130503085000.sps";
infile=OpenRead[filename,BinaryFormat->True];
outfile=OpenWrite[StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_r200.sps",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];
WriteString[outfile,chartversion];

chartstart="";
chartstart=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartstart,"Real64"];

chartend="";
chartend=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartend,"Real64"];

chartlat="";
chartlat=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlat,"Real64"];

chartlong="";
chartlong=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlong,"Real64"];

chartmax="";
chartmax=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmax,"Real64"];

chartmin="";
chartmin=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmin,"Real64"];

chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartTZ,"Integer16"];

chartsource="";
chartsource=ReadString[infile,10];
WriteString[outfile,chartsource];

chartauthor="";
chartauthor=ReadString[infile,20];
WriteString[outfile,chartauthor];

chartlocalname="";
chartlocalname=ReadString[infile,20];
WriteString[outfile,chartlocalname];

chartlocation="";
chartlocation=ReadString[infile,40];
WriteString[outfile,chartlocation];

chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartchannels,"Integer16"];

chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
BinaryWrite[outfile,chartnotelength,"Integer32"];

chartnote="";
chartnote=ReadString[infile,chartnotelength];
WriteString[outfile,chartnote];

(* Read number of frequency sweeps from neader notes            *)
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Create sweep arrays, then populate, rotate, write to file    *)
(* iterating over all sweeps in source file.                    *)
indataB=Table[0,{i,1,chartchannels}];
indataA=Table[0,{i,1,chartchannels}];
Do[
{
(* Read one frequency sweep                                  *)
Do[
{
indataB[[i]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
indataA[[i]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Rotate channel 200 back down to below channel 199         *)
indataB=RotateLeft[indataB];
indataA=RotateLeft[indataA];
(* Write data to output file                                 *)
Do[
{
BinaryWrite[outfile,indataB[[i]],"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataA[[i]],"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Write EOR marker                                          *)
indataB[[1]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataB[[1]],"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps}
];
Close[infile];
Close[outfile];
Print["Done."]