Search code examples
signal-processingfftoctavevibration

Car vibration evaluated by Fast Fourier Transform in Octave


I must evaluate vibration of a car. For this trial I used accelerometer. Collected data are depended on the time.

I need to convert the data from time domain to frequency domain by FFT. Unfortunately, I am not familiar with coding and FTT very well, however I found and used the code below.

What is strange for me, that maximum high point has 0Hz. Please, see attached image. Anyway, is there a way how to do a graph more obvious? For example, cut a series in x-axis and shows only data under 200Hz.

enter image description here

clc
A=xlsread('50_dirt_road.xlsx');
t=A(:,9);
s=A(:,8);
Ts = mean(diff(t));                                     % Sampling Interval
Fs = 1/Ts;                                              % Sampling Frequency
Fn = Fs/2;                                              % Nyquist Frequency
L = numel(t);                                           % Signal Length
sm = s - mean(s);                                       % Mean-Corrected Signal (Eliminates 0 Hz Offset)
FTs = fft(sm)/L;                                        % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:numel(Fv);                                       % Index Vector
[MaxV,idx] = max(abs(FTs(Iv))*2);                       % Maximum V & Index
Freq = Fv(idx);                                         % Frequency Of Maximum V
figure
plot(Fv, abs(FTs(Iv))*2)
grid
text(Freq, MaxV, sprintf('\\leftarrow %.4f G, %.0f Hz', MaxV, Freq), 'HorizontalAlignment','left')
xlabel('Frequency (Hz)')
ylabel('Amplitude')

Could you double check it, please? My variables are defined as follows:

  • s: measured "G" value. In total 3395 measured values.

  • t: time. Every single value was recorded after 0.001s, in total 3.395s.


Solution

  • The maximum high point you are referring to is normal. Without having frequency cut offs, the beginning, (nearest zero), of your fourier transform will always have the most spectral energy because it is resolving sampling rates tending towards zero which as you can imagine would have a ton of representation in the time series. Obviously 0 hz is impossible, so you should choose a lower bound frequency that makes sense in the domain of car vibrations. So I recommend a band pass filter, one that does forward and backward filtering (filt filt) to preserve the original time series and then running the result through this analysis. I would start with a butterworth filter and experiment with othes if desired:

    https://octave.sourceforge.io/signal/function/butter.html

    EDIT:

    I was using t but the signal is in s. Here is the whole thing:

    clc
    A=xlsread('50_dirt_road.xlsx');
    t=A(:,9);
    s=A(:,8);
    Ts = mean(diff(t));                                     % Sampling Interval
    Fs = 1/Ts;                                              % Sampling Frequency
    Fn = Fs/2;                                              % Nyquist Frequency
    L = numel(t);                                           % Signal Length
    
    %1st order butterworth filter with a band pass of 1hz to 200hz in radians
    %forward and reverse filtered
    [b,a] = butter(1, [1/(L/2), 200/(L/2)]);
    filtered_s = filtfilt(b,a,s);
    
    sm = filtered_s - mean(filtered_s);                     % Mean-Corrected Signal (Eliminates 0 Hz Offset)
    FTs = fft(sm)/L;                                        % Fourier Transform
    Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
    Iv = 1:numel(Fv);                                       % Index Vector
    [MaxV,idx] = max(abs(FTs(Iv))*2);                       % Maximum V & Index
    Freq = Fv(idx);                                         % Frequency Of Maximum V
    figure
    plot(Fv, abs(FTs(Iv))*2)
    grid
    text(Freq, MaxV, sprintf('\\leftarrow %.4f G, %.0f Hz', MaxV, Freq), 'HorizontalAlignment','left')
    xlabel('Frequency (Hz)')
    ylabel('Amplitude')
    

    LAST EDIT:

    Okay, so I think I got you an intro in to the world of filtering when I think all you were looking to do is constrain the axis of your FFT. Same arguments for the band pass filter above, 1hz and 200hz. The above code should work but the following code is probably what you were looking for originally:

    clc
    A=xlsread('50_dirt_road.xlsx');
    t=A(:,9);
    s=A(:,8);
    Ts = mean(diff(t));                                     % Sampling Interval
    Fs = 1/Ts;                                              % Sampling Frequency
    Fn = Fs/2;                                              % Nyquist Frequency
    L = numel(t);                                           % Signal Length
    
    sm = s - mean(s);                     % Mean-Corrected Signal (Eliminates 0 Hz Offset)
    FTs = fft(sm)/L;                                        % Fourier Transform
    Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
    freqMask = (Fv > 1) & (Fv < 200);
    Fv = Fv(freqMask);
    FTs = FTs(freqMask);
    Iv = 1:numel(Fv);                                       % Index Vector
    [MaxV,idx] = max(abs(FTs(Iv))*2);                       % Maximum V & Index
    Freq = Fv(idx);                                         % Frequency Of Maximum V
    figure
    plot(Fv, abs(FTs(Iv))*2)
    grid
    text(Freq, MaxV, sprintf('\\leftarrow %.4f G, %.0f Hz', MaxV, Freq), 'HorizontalAlignment','left')
    xlabel('Frequency (Hz)')
    ylabel('Amplitude')