Message Boards Message Boards

Recorded accelerometer data to PSD plot?

Posted 10 years ago


I would like to take data recorded from an accelerometer during a random vibration test and generate a PSD plot from it.

It would be a log=-log plot with the vertical axis in G^2/Hz and the horizontal axis in Hz.

The data is a CSV file of pairs: {time, accel}

I'm new to Mathematica and I am sure this is not that hard but after hours of searching and contacting tech support multiple times I have not been successful. I think I am close using the Periodgram or the PowerSpectralDensity, but it never seems to quite work out. The horizontal axis seems to be in radians and not really modifiable. I did this about 10 yrs ago in MathLab and it was a snap.

Has anyone been successful in doing this and have some insight on how to do this?

thanks in advance


POSTED BY: George Thiel
3 Replies
Posted 9 years ago

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}

POSTED BY: Caitlin Ramsey
Posted 10 years ago


thank you so much for the very detailed reply. I will be trying it shortly.

thanks again,


POSTED BY: George Thiel

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:

  myseries = N[Import["test.txt", "Table", "FieldSeparators" -> {","}]];

The enclosing N[] is to convert any integer numbers to machine precision (for compute efficiency later) 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"}]
POSTED BY: Frank Iannarilli
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract