0
|
21813 Views
|
3 Replies
|
0 Total Likes
View groups...
Share

# Recorded accelerometer data to PSD plot?

Posted 10 years ago
 Hi, 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 George
3 Replies
Sort By:
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 10 years ago
 Frank,thank you so much for the very detailed reply. I will be trying it shortly.thanks again,George
Posted 10 years ago
 Hi George, I'll partition this problem into 2 pieces: (1) importing the data file; (2) taking PSD and plottingLet'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 screenDimensions[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"}]