Wow, I found this explanation really helpful, too. I am working on an accelerometer-data challenge: calculating approximate position from accelerometer data. It's a cyclic motion with vibration, and I've succeeded in applying your method to see the major frequencies in my data. I wonder if you could offer guidance on how I might approximate a "fit" to the major frequencies that could be double-integrated to compute approximate x, y, z positions of the device? It's an irregularly shaped rotor in a system with some translational motion. I've seen approaches to compute the translational component (eg. position of someone walking) by subtracting out the cyclical motion and vibration -- and these methods get pretty hairy. But in my case, the cyclical motion is what would be of interest (ie. fixed position on an irregularly shaped rotor) rather than the translational motion of the whole system. Any suggestions?

P.S. Data in same format as example above: {time, acceleration}

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:

1. tsamp, which is the time interval (in seconds) between samples
2. 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"}]