Engee documentation
Notebook

Third-octave analysis for railway transport

Third-octave analysis is a method of frequency analysis of signals in which the spectrum is divided into bands, each of which is 1/3 octave. The use of third-octave analysis for trains makes it possible to study noise and vibrations and assess their impact on the environment, passengers and infrastructure. This method is especially important for:

  • Noise normalization (compliance with GOST, ISO, EN).

  • Diagnostics of the condition of the railway track and rolling stock.

  • Optimization of structures (for example, reduction of noise from wheelsets or rails).

Accelerometers**mounted on wheel sets, wagon frame, rails, and ground are often used for analysis (to assess the propagation of vibrations in buildings). For railway transport, the key third-octave bands usually lie in the range of 1-1000 Hz (vibration) and 20-5000 Hz (noise). Examples of important frequencies:

  • Low (1-80 Hz): Vibrations from wheel impacts on rail joints.

  • Medium (80-500 Hz): Rolling noise, engine hum.

  • High (500-5000 Hz): Creaking brakes, aerodynamic noise.

The construction of a one-third octave spectrum allows you to analyze peaks in certain bands that may indicate:

  • Wheel wear (increased level in the 63-250 Hz bands).

  • Rail defects (resonances at 31.5–125 Hz).

  • Suspension problems (2-20 Hz vibrations).

Let's consider the problem of analyzing real data from an accelerometer on railway tracks and constructing a one-third octave spectrum.

Data import and preprocessing

Connecting the necessary libraries:

In [ ]: 
using DSP, DelimitedFiles, Statistics, FFTW, Polynomials

Accelerometer data is stored in a file data.txt. We count them as a matrix in a variable datamat:

In [ ]: 
datamat = readdlm("data.txt", ',');

Let's extract from it separate vectors for time and "raw" vibration acceleration values (in m/s^2), and display them on a graph in the time domain.:

In [ ]: 
timevec = datamat[:,1];
data1 = datamat[:,2];
Plots.plot(timevec,data1,
            title = "Импортированные данные виброускорений",
            xlabel = "Время, с",
            ylabel = "Виброускорение, м/с^2",
            legend = false)
Out[0]:

We see that the recording of vibration acceleration values did not begin immediately at the initial moment of time, so it is logical to cut off part of the signal level close to zero. To do this, we will leave only the vector indexes corresponding to the signal values exceeding a certain threshold (in our case, 1).:

In [ ]: 
indices = findall(x -> x > 1, data1);
accel = data1[indices[1]:indices[end]];
t = timevec[1:length(accel)];
Plots.plot(t, accel,
            title = "Вырезанные значения виброускорений",
            xlabel = "Время, с",
            ylabel = "Виброускорение, м/с^2",
            legend = false)
Out[0]:

Transition to vibration speeds

To visualize vibration velocities, we need to integrate the vector of acceleration values. accel. Let's write a custom cumulative trapezoidal numerical integration function. cumtrapz, an analog of the one used in the environment MATLAB:

In [ ]: 
function cumtrapz(X::T, Y::T) where {T <: AbstractVector}
    @assert length(X) == length(Y)
    out = similar(X)
    out[1] = 0
        for i in 2:length(X)
            out[i] = out[i-1] + 0.5*(X[i] - X[i-1])*(Y[i] + Y[i-1])
        end
    out
  end
Out[0]:
cumtrapz (generic function with 1 method)

Let's apply a custom function to a vector accel, and then we will select a constant component (or trend) in the vibration velocity signal in order to remove it from the signal by element-wise subtraction, if necessary. To do this, we approximate the vibration velocity signal with a 12th-order polynomial and estimate its values at the same time points as the velocity signal.:

In [ ]: 
velocity_trend = cumtrapz(t,accel);
poly =  fit(t, velocity_trend, 12);
trend = poly.(t);
Plots.plot(t, velocity_trend, label = "Виброскорость",
            title = "Виброскорость и её тренд",
            xlabel = "Время, с")
Plots.plot!(t, trend, linewidth = 3, label = "Постоянная составляющая",
            ylabel = "Виброскорость, м/с")
Out[0]:

Spectral analysis by the Welch method

We will perform a spectral analysis of vibration accelerations using standard library packages. DSP.jl. First, let's determine the value of the sampling frequency fs the original signal in Hz, analyzing the step of the time vector dt:

In [ ]: 
dt = mean(diff(timevec));
fs = Int(floor(1/dt))
Out[0]:
3472

Let's set the parameters for spectrum analysis that were included in the terms of reference. In particular: the size of the "window" for the input of the FFT algorithm, the number of zeros for interpolating the power spectral density graph (SPM), an overlap of 50% for the Welch method, the output number of FFT samples (the size of the "window" + the number of zeros).

In [ ]: 
window_length = fs;
nzeros = 100000;
overlap = 0.5;
nFFT = window_length + nzeros;
noverlap = Int(floor(overlap*window_length));
stepsize = window_length - noverlap;

It is worth noting that the number of points of the input "window" is equal to the sampling frequency. This makes it possible to estimate the spectrum with a resolution of 1 Hz.

Also, according to the assignment, the analysis is performed only in the lower half of the spectrum (the first Nyquist zone) of positive frequencies. To do this, the signal is passed through a 12th-order Butterworth low-pass filter, and its normalized cutoff frequency is 0.5. The filtering result is recorded in a variable filtered:

In [ ]: 
Wp = 0.5;
myfilter = digitalfilter(Lowpass(Wp), Butterworth(12));
filtered = filtfilt(myfilter,accel);

Now let's use the function welch_pgram to estimate the spectral power density using the Welch method. The overlap parameters, the size of the input "window" and the output number of FFT points correspond to those previously set. Let's display the SPM charts in dB:

In [ ]: 
 pxx = welch_pgram(filtered, n = stepsize, noverlap = div(stepsize, 2); 
                onesided = true,
                nfft = nFFT, 
                fs = fs, 
                window = hanning);
pxx_nofilt = welch_pgram(accel, n = stepsize, noverlap = div(stepsize, 2); 
        onesided = true,
        nfft = nFFT, 
        fs = fs, 
        window = hanning);
Plots.plot(freq(pxx), pow2db.(power(pxx)), label = "После фильтра",
                title = "СПМ сигнала виброускорений",
                xlabel = "Частота, Гц")
Plots.plot!(freq(pxx_nofilt), pow2db.(power(pxx_nofilt)), label = "До фильтра",
                ylabel = "Спектральная плотность мощности, дБ/Гц")
Out[0]:

We will also display the amplitude spectrum in units of measurement of the signal. To do this, we need to determine the parameter of the equivalent noise band for the Hamming window function, which we used in the Welch method.:

In [ ]: 
mywin = hanning(window_length);
ENBW = fs * (sum(abs.(mywin).^2)) / (abs.(sum(mywin)).^2)
Out[0]:
1.5004321521175445

We get the amplitude value from the power value of the spectrum, taking into account the parameter ENBW:

In [ ]: 
amplitude = sqrt.(2*power(pxx)*ENBW);
Plots.plot(freq(pxx), amplitude, legend = false, 
            title = "Амплитудный спектр сигнала виброускорений",
            ylabel = "Амплитуда, м/с^2",
            xlabel = "Частота, Гц",
            xlim = (0,1000))
Out[0]:

The third octave spectrum

To construct a one-third octave spectrum, it is first necessary to determine the boundaries of the frequency ranges and their number. count:

In [ ]: 
nstart = -20;         
ncount = nstart;
count = 0;
fmax = 0;
while fmax*(2^(1/6)) <= fs/2
    fmax = (1000).*((2^(1/3)).^ncount);
    ncount = ncount + 1;
    count = count + 1;       
end

Let's calculate the average power in each of the ranges (vectors ff) in a loop and write the results into a vector octave_power:

In [ ]: 
ff = (1000).*((2^(1/3)).^(nstart:(ncount-1)));
octave_power = zeros(count);
psdx = power(pxx);
for i = 1:count
    f_low = ff[i]/(2^(1/6));
    f_high = ff[i]*(2^(1/6));
    idx = .&(freq(pxx) .> f_low, freq(pxx) .<= f_high);
    idx = vec(idx);    
    octave_power[i] = mean(psdx[idx])
end

Let's move on to the values of the average amplitude in the bands and display the third-octave spectrum with the function bar:

In [ ]: 
octave_amplitude = sqrt.(2*octave_power*ENBW);
freq_ticks = round.(ff, digits = 0);
BW = exp.(0.2:0.24:0.24*count);
bar(freq_ticks, octave_amplitude, bar_width = BW, xaxis=:log,
                 xticks = (ff, freq_ticks), color = :green,
                 title = "Третьоктавный спектр виброускорений",
                 xlabel = "Частота, Гц",
                 ylabel = "Амплитуда, м/с^2",
                 legend = false)
Out[0]:

Conclusion

Third-century analysis is a common method of diagnosing the web or wheel sets of rolling stock, which allows detecting, for example, uneven wear of rails or problems with wheel–rail contact. Spectrum analysis in a similar range is prescribed in GOST 31319-2009 and GOST 31296.1-2019.

We have reviewed the possibilities of the available tools of the Engee environment for conducting such an analysis.