Hi George,
I'll partition this problem into 2 pieces: (1) importing the data file; (2) taking PSD and plotting
Let's say your file looks like this:
time1, accel1
time2, accel2,
...
timeN, accelN
where from your mention of CSV, I assume the pair elements are comma-separated as shown, and each pair on a separate line.
So importing:
SetDirectory["C:\\Projects\\Accelerometer"]
myseries = N[Import["test.txt", "Table", "FieldSeparators" -> {","}]];
The enclosing N[] is to convert any integer numbers to machine precision (for compute efficiency later)mailto:wolfram-community@wolfram.com?subject=Wolfram%20Community%20Feedback
The terminating ; is to prevent the imported data from spewing to the screen
Dimensions[myseries]
should yield {n,2} and myseries is organized as nested Lists as so:
{{time1, accel1}, {time2, accel2},...{timeN, accelN}}
From here it's easy:
I'll define 3 psdFunc[] functions, the latter two providing you a PSD frequency-smoothing option. Evaluate all 3.
Consider a rate signal in units of quantity/sec, and (averaged-sample) duration t sec. Output later needs to be scaled by t to get PSD in units of rate^2/Hz.
psdFunc[seq_] := (Abs[Fourier[seq, FourierParameters -> {1, -1}]]^2)/Length[seq]
We can average in the frequency domain by taking the PSDs of non-overlapping blocks of samples, then averaging the PSDs:
psdFunc[seq_, nBlocks_] := (If[Mod[Length[seq], nBlocks ] != 0,
Print["Non-integer multiple nBlocks!"]; Return[]];
Apply[Plus, Table[psdFunc[seq[[Range[1 + (n - 1)*Length[seq]/nBlocks, (n)*Length[seq]/nBlocks]]]], {n, 1, nBlocks}]]/nBlocks);
Or better yet (to minimize sidelobe leakage), smooth the PSD in frequency domain by applying window therein (and subsampling).
psdFunc[seq_, nBlocks_] := (If[Mod[Length[seq], nBlocks ] != 0,
Print["Non-integer multiple nBlocks!"]; Return[]];
Partition[psdFunc[seq], nBlocks]).Table[1./nBlocks, {nBlocks}];
Now set two parameters in the expression below -- I've put in sample values, but you should re-check/correct for your data:
- tsamp, which is the time interval (in seconds) between samples
- fbin, the number of raw frequency bins (of frequency resolution 1/bigT) to average
You may also want to truncate the series length to a convenient multiple of 2 as the length of myseries:
myseries=myseries[[1;;512>, All]];
tsamp = myseries[[2,1]]-myseries[[1,1]];
nrec = Length[myseries];
bigT = nrec*tsamp;
fbin = 2;
ListPlot[tsamp*psdFunc[myseries[[All, 2]], fbin][[1 ;; nrec/2/fbin]],
Joined -> True, PlotRange -> All, DataRange -> {0, 1/tsamp/2},
GridLines -> Automatic, Frame -> True,
FrameLabel -> {"Frequency (Hz)", "G^2/Hz", "Accelerometer PSD"},
PlotStyle -> {{Blue, Thickness[0.005`]}, {Brown, Thickness[0.005`]}},
BaseStyle -> {FontSize -> 16, FontFamily -> "Arial"}]