Hi Neil, thank you again very much for your great help!

Hi John,

What I want to do is to develop a random vibe testing PSD profile based on the actual accelerometer data to realistically shake my devices. PSD plot is a power spectrum density plot.

Two examples of Matlab/Python script generated PSD plots from the data I shared for 1) sps = 8192, df=1.22Hz, dof = 2, and 2) sps=1024, df=9.77Hz, dof=18, are shown below (someone made them for me as reference). In general, the number of spike should be reduced on PSD so we need to select the right sps (sample per segment) from the data. It looks like sps 1024 case is better than sps 8192 case for this dataset. I am still working on Neil's codes to see if I can reproduce identical plots by Mathematica. Thanks!

Delan

Delan,

Yes, fs is just the inverse of tsamp. I used the fs notation because I took the simple example from your file, "Vibration_Analysis_Examples.m". The use for fs or tsamp is to figure out Nyquist to know what the highest frequency on the plot is numerically (your FFT runs from 0 hz to 1/(2*fs) -- Nyquist, if you plot just the first half of the data; or 0hz to 1/fs if you plot the repeating second half (as I did). In MATLAB the code was written to show only one half (which makes sense) and plots it from 0 to Nyquist and multiplies the magnitude by 2 so the magnitude agrees with the coefficient of the Sin waves. You can duplicate this in MMA easily.

For varied length/rate cases you should adjust the fs so that you correctly compute the maximum value on the x axis. As you probably already know, The resolution of your FFT is based on the length of the time waveform (the longer the waveform, the higher the resolution)

By running the example from the MATLAB code, you can make sure your MMA code gives the same answer before you run it on your live data. (They make certain assumptions and use certain conventions so you might want to match that).

Lastly, if you download MATLink for MMA, you can always run MATLAB code from MMA if that helps. See this Thread for downloading it as the website is temporarily down

Regards,

Neil

Here is the FFT of the Sample data (Note the factor of 2 from the Matlab code)

And the Power (note the power is the FFT mag squared)

Delan,

Your code is close to correct.

1. There is no need in Mathematica (MMA) to use a power of 2 for the FFT (especially such a short FFT as you have)

2. Your time data in the CSV has rounding issues so tsamp. I would round it so it is 0.0001 (Maybe best to take the last value and the first and divide by the length or use Round)

3. For this type of analysis your Fourier options should be FourierParameters -> {-1, 1}]

4. The power spectrum is calculated

   Abs[Fourier[xx, FourierParameters -> {-1, 1}]^2]
   

5. This version is sample rate/length invariant. You could use the built in function PowerSpectralDensity[xx, w, FourierParameters -> {-1, 1}] to get a continuous function of w but MMA's version would need to be scaled for length and the computation gets long for long data sets.

6. Note that there is a factor of 2 between the Matlab and the MMA FFT because in Matlab they truncated it and mapped it to remove the mirrored side of the spectrum.

To verify your code, run the example you posted:

f1 = 60;
f2 = 22;
f3 = 100;
fs = 500;
xx = Table[
   2*Sin[2*Pi*f1*t] + 1*Sin[2*Pi*f2*t] + 1.5*Sin[2*Pi*f3*t], {t, 0, 
    0.5 - 1/fs, 1/fs}];
ListLinePlot[xx]
(* FFT *)
ListPlot[Abs[Fourier[xx, FourierParameters -> {-1, 1}]], 
 PlotRange -> All, ImageSize -> Large, Filling -> Axis, 
 DataRange -> {0, fs}]
(* Power Spectral Density *)
ListPlot[Abs[Fourier[xx, FourierParameters -> {-1, 1}]^2], 
 PlotRange -> All, ImageSize -> Large, Filling -> Axis, 
 DataRange -> {0, fs}]

I hope this helps,

Regards,

Neil

Hi Delan,

Have you looked at Periodogram? Can you post the Matlab / Python code and results?

Rohit