Max velocity on 1/3 octave band from Acceleration-Time data

I have acceleration time data taken by accelerometer. I need to find the peak velocity of each frequency on the 1/3 octave band range (4, 5, 6.3, 8 etc). im confused whether should i do since im new into matlab. ive done several script (help from GPT), but weird results.

[Mathieu NOE]

hello

can you provide the data file as well ?

[Gerald]

sorry for the late response. here is the sample data

Answer 1

[Mathieu NOE]

hello again

so I decided to make already some modifications in your code , and test it with synthetic data (at least you know what the result should be so it's easier to chase down the most evident mistakes)

as usual with fft , you have to pay attention how to normalize your fft output so for example a sinewave of amplitude 1 (peak amplitude) should give you the same amplitude once fft'ed. that was the major hurdle in your code .

the others mods are less critical , for example when doing the time integration, there is always the risk that the output is not zero-mean even though you have applied a bandpass filter on the incoming signal. Beware that transients can cause drift . So I prefer to simply do the bandpass filtering after the cumtrapz operation, so we are on the safe side and we know that now the output signal is balanced. In a further refinement you could even remove some of the signal that correspond to the transient time response of the bandpass filter.

with the example below a unitary amplitude (sine) acceleration signal at 10 Hz is used. Converted to velocity the amplitude becomes 1/(2*pi*10) = 0.016 (whatever the units are) which is what we obtain in both fft and 1/3 octave bands

%% 1. Load Acceleration Data
% filename = 'Imperial_Valley.asc'; % Change to your actual filename
% data = readmatrix(filename, 'FileType', 'text');
% acc = data(:,1) * 10; % Convert cm/s² to mm/s²

% some synthetic data to test the code
N = 1e4;
Fs = 200; % Sampling frequency (adjust if needed)
dt = 1 / Fs; % Time step
t = (0:N-1) * dt;
acc = sin(2*pi*10*t); % a unitary amplitude, 10 Hz sine wave

%% 2. Define Parameters
Fs = 200; % Sampling frequency (adjust if needed)
N = length(acc);
dt = 1 / Fs; % Time step
t = (0:N-1) * dt;
freqs = [4, 5, 6.3, 8, 10, 12.5, 16, 20, 25, 31.5, 40, 50, 63]; % 1/3-octave bands

%% 3. Apply High-Pass Filter (Remove Low-Frequency Drift)
d_hp = designfilt('highpassiir', 'FilterOrder', 4, ...
    'HalfPowerFrequency', 0.5, 'SampleRate', Fs);
acc_filtered = filtfilt(d_hp, acc);

%% 4. Compute FFT-Based Velocity Spectrum
N_pad = 2^nextpow2(N); % Zero-padding for better frequency resolution
Y = abs(fft(acc, N_pad))/N; % /!\ normalized fft !!
select = (1:N_pad/2);
Y = Y(select)*2; % Single-sided spectrum
f = (select - 1)*Fs/N_pad;
% remove freq = 0 data
ind = f>0;
f = f(ind);
Y = Y(ind);
% % debug plot
% figure, plot(f, abs(Y))

% Convert acceleration to velocity in frequency domain
omega = 2 * pi * f;
V_fft = Y ./ omega;

%% 5. Compute 1/3-Octave Band Velocity Levels
velocity_levels = zeros(size(freqs));

for i = 1:numel(freqs)
    % Define 1/3 octave band edges
    lower_freq = freqs(i) / (2^(1/6));
    upper_freq = freqs(i) * (2^(1/6));
    
    % Design bandpass filter for the current band
    d = designfilt('bandpassiir', 'FilterOrder', 4, ...
        'HalfPowerFrequency1', lower_freq, ...
        'HalfPowerFrequency2', upper_freq, ...
        'SampleRate', Fs);

%     % Apply filter
%     filtered_signal = filtfilt(d, acc_filtered);
    
    % Convert acceleration to velocity (integrate)
    velocity_signal = cumtrapz(t, acc_filtered); % Integration
    
    % Apply filter after integration (to remove DC  / drift mostly) 
    velocity_signal = filtfilt(d, velocity_signal);

    % Get peak velocity for this band
    velocity_levels(i) = max(abs(velocity_signal));
end

%% 6. Plot FFT Velocity Spectrum
figure;
plot(f, V_fft, 'b');
xlabel('Frequency (Hz)');
ylabel('Velocity (mm/s)');
title('FFT-Based Velocity Spectrum');
grid on;
xlim([0, max(freqs) + 5]);

%% 7. Plot 1/3-Octave Band Velocity Spectrum
figure;
semilogx(freqs, velocity_levels, 'o-r'); % Log scale for better visibility
xlabel('Frequency (Hz)');
ylabel('Peak Velocity (mm/s)');
title('1/3-Octave Band Velocity Spectrum');
grid on;

[William Rose]

@Gerald,

@Mathieu NOE provided a very thoughtful and good answer that took a real effort. It is nice to "accept" an answer in that case, so I am accepting his answer on your behalf.

[Mathieu NOE]

hello @William Rose

that's very kind , I really appreciate, - but you too has spend time on the topic so I hope @Gerald will at least vote for your contribution !

all the best for both of you

Answer 2

[William Rose]

@Gerald,

I have been doing the same tyhing as @Mathieu NOE b ut not as quicjkly.

One thing I noticed is that the line

V_fft = Y ./ (1i * omega);

doesn't do what you want. It makes a 2000x2000 array. Do

V_fft = Y ./ (1i * omega');

instead.

I made a simulated velocity waveform with three pure sinusoids, each with amplitude=1, at frequencies 10, 20, and 40 Hz. Remember that this is assumed to be velocity in cm/s.

Fs=200; dt=1/Fs; t=(0:2000)'*dt;
v=cos(2*pi*10*t)+cos(2*pi*20*t)+cos(2*pi*40*t);

I computed the acceleration and saved the result as a file.

a=diff(v)/dt;
save('Imperial_Valley0.asc','a','-ascii')

I imported the simulated data file with load() instead of readmatrix() since I got funny results with readmatrix().

I removed the zero-padding because it does not increase the true frequency resolution. Zero-padding the time-domain signal results in a spectrum which is, in effect, interpolated to a higher sampling rate (in the frequency domain) from the spectrum of the non-zero-padded signal. But this seeming improvement in frequency reoslution is only due to interpolation, and therefore gives a misleading indicaiton of the actual frequency resolution. To get a true improvement in frequency resolution, you need to sample for a longer period of time.

The use of filtfilt() with a high-pass filter, to remove low frequency drift, gives bad results, at least for my simulated data file. See plot below of the simulated acceleration signal, before and after filfilt().

Therefore I recommend using

acc_filtered=detrend(acc,2);

to remove quadratic trend (or 3, to remove cubic trend, etc) from the signal. If you plot the acceleration signal before and after detrend(), you will see that the unwanted start and end transients, which occur with filtfilt(), do not occur with detrend().

More later.

[Mathieu NOE]

also I wonder why we use a time approach (integration) to do the ocatve band analysis as we have already performed a fft

there is a more direct path to convert fft to octave (or Nth) bands (and do the integration in frequency domain must faster)

for example :

NarrowToNthOctave - File Exchange - MATLAB Central

[Gerald]

@Mathieu NOE

are you trying to say that, obtaining the max velocity on 1/3 octave bands shouldve been simpler than this?

[Mathieu NOE]

@Gerald

yep , but maybe this is a debate if you want the rms value in a band or if you really want the "true" peak velocity . In the older days of analog octave band analysers , the signal was split by fixed analog filters and then the rms value was derived , but the ratio between rms and peak depends of the nature of the signal. so it's easier (and definitively more direct) if you only want rms . The question is what does the "exact" peak velociy bring as additionnal info vs. the rms value ?

[Gerald]

since i need it for impact analysis for a structure. i need the true peak velocity. and what we've been doing and disccussed are for the true peak right? since rms should be averaging each freq inside the band for each

[William Rose]

@Gerald,

When you integrate acceleration to estimate velocity, and then use max(abs(...)) to find the peak velocity, you can get misleading results. Consider the following example: a 10 Hz sinusoidal velocity signal, with amplitude=1.

Fs=200; dt=1/Fs; t=(0:50)'*dt;
v=cos(2*pi*10*t);       % velocity (amplitude=1)
a=diff(v)/dt;           % acceleration
v1=cumtrapz(dt,a);      % velocity estimated from acceleration
% Plot velocity, estimated velocity, and acceleration.
figure; subplot(211), plot(t,v,'-r.',t(2:end),v1,'-b.')
grid on; legend('original v','v from a'); ylabel('Velocity')
subplot(212), plot(t(2:end),a,'-r.')
grid on; ylabel('Acceleration'); xlabel('Time (s)')

Estimate peak velocity from the velocity signal derived by integration:

vpeak=max(abs(v1));
fprintf('Estimated max velocity=%.2f.\n',vpeak)

Estimated max velocity=1.95.

The peak velocity of the original signal is +-1. Therefore the approach above overestimates the peak by a factor of 2 in this case, because the reconstructed velocity trace varies from 0 to -2, while the original trace varies from +1 to -1. This kind of error can happen in your code. Since you are reconstructing signals in frequency bands that do not extend to f=0, you may remove the mean value of the signal you get from cumtrapz():

v1=cumtrapz(dt,a);  % this line is same as above
v1=v1-mean(v1);     % remove the mean value

With this change, you will get a better estimate for the peak velocity.

vpeak=max(abs(v1));
fprintf('Estimated max velocity=%.2f.\n',vpeak)

Estimated max velocity=0.98.

OK.

[William Rose]

@Gerald,

I always learn from the code and comments of @Mathieu NOE.

In the version attached here, I keep the original order which you had: band-pass filter the acceleration, then reconstruct velocity for that band, by integrating the band-passed acceleration. Then remove the mean value of the velocity signal, as discussed in my previous post. It gives the results shown below.

The original velocity signal, in this example, has three cosine waves, at frequencies 10, 20, 40 Hz, with amplitude=1 cm/s for each component. Therefore the expected amplitude is 10 mm/s, at f=10, 20, 40 Hz, and zero at other frequencies. The plot aboive shows amplitude ~=10 at f=10, 20, 40, as expected. The amplitude is considerably lower, as expected, at intermediate bands with center frequencies above 10 Hz. For bands with center frequencies < 10 Hz, the estimated peak velocity is higher than expected. This is probably due to a combination of non-ideal band-pass filters and imperfect trend removal, etc.

By the way, if you want to estimate amplitudes for bands with center frequencies at 50 Hz and higher, I recommend using a sampling rate several times greater than 200 Hz. Consider a sinusoid with f=63 Hz. If Fs=200 Hz, this sine wave is only sampled about three times per period, which is probably not enough to get an accuarate estimate of its amplitude. Theory says you'll be OK as long as your signal frequency is below the Nyquist frequency, but in practice, oversampling is useful.

I changed the data file extension from .asc to .txt, since Matlab Answers won't let me attach .asc.

[Gerald]

so is it normal that there is some drift from expected value (10 mm/s) due to combination of filters etc? is there a way to fix the problem? Also many thanks for guiding me and giving me suggestion. Super appreciated it @Mathieu NOE @William Rose

[Mathieu NOE]

hello again

that's why I put the (only) HP filter after the cumptrapz integration , so it removes DC & drifts that are both in the input signal and also occasionnaly caused during the integration process

[Mathieu NOE]

@Gerald

since i need it for impact analysis for a structure. i need the true peak velocity.

are you saying you are doing hammer tests ? FRF ?

this may interest you in that case : How to make an FRF with my accelerometer Data and FFT? - MATLAB Answers - MATLAB Central

of course you can adapt the solution to have a FRF in velocity vs force format (instead of acceleration vs force)

[Gerald]

not really a hammer test. its a train generated vibration. the standard for this assesment in my region is using the peak velocity on the 1/3 bands

[Gerald]

Hello again both.

i attached my almost finished codes for this issue.

i tried @Mathieu NOE earlier suggestion which is inserting the filter after cumtrapz. it seems that the result are more or less the same.

also, in my code which to visualise the velocity-time graph. i used the original acc instead of acc_detrend because it becomes parabolic (up). is it really like that?

this is the original acc

this is the acc_detrend

[Gerald]

Oh here is the data i used. actually im using simple sine provided by @William Rose but changed the file name.

txt because matlab says no to asc

[Mathieu NOE]

ok , good to know

just as I mentionned above the octave analysis can be done with filters (as you do) , which is the most precise method , but more and more the octave is simply derived from fft (see for example : Octave Analysis in ObserVIEW - FFT and Filter-based - Vibration Research)

to come back to your code , and to test it with 3 tones like @William Rose , my simple suggestion here is to simply discard the transients introduced by the filters (here I arbitrary choose to remove 20% at both ends)

NB that the ocatve plot is no more "poluted" by low DC / drift

all the best

%% 1. Load Acceleration Data
% filename = 'Imperial_Valley.asc'; % Change to your actual filename
% data = readmatrix(filename, 'FileType', 'text');
% acc = data(:,1) * 10; % Convert cm/s² to mm/s²

% some synthetic data to test the code
N = 1e4;
Fs = 200; % Sampling frequency (adjust if needed)
dt = 1 / Fs; % Time step
t = (0:N-1)' * dt;
% acc = sin(2*pi*10*t); % a unitary amplitude, 10 Hz sine wave

acc = sin(2*pi*10*t)+sin(2*pi*20*t)+sin(2*pi*40*t); % unitary amplitude sine waves

%% 2. Define Parameters
Fs = 200; % Sampling frequency (adjust if needed)
N = length(acc);
dt = 1 / Fs; % Time step
t = (0:N-1) * dt;
freqs = [4, 5, 6.3, 8, 10, 12.5, 16, 20, 25, 31.5, 40, 50, 63]; % 1/3-octave bands

%% 3. Apply High-Pass Filter (Remove Low-Frequency Drift)
d_hp = designfilt('highpassiir', 'FilterOrder', 4, ...
    'HalfPowerFrequency', 0.5, 'SampleRate', Fs);
acc_filtered = filtfilt(d_hp, acc);

%% 4. Compute FFT-Based Velocity Spectrum
N_pad = 2^nextpow2(N); % Zero-padding for better frequency resolution
Y = abs(fft(acc_filtered, N_pad))/N; % /!\ normalized fft !!
select = (1:N_pad/2);
Y = Y(select)*2; % Single-sided spectrum
f = (select - 1)*Fs/N_pad;
% remove freq <1 Hz data (up to you to change the cut off freq)
ind = f>1;
f = f(ind);
Y = Y(ind);

% Convert acceleration to velocity in frequency domain
omega = 2 * pi * f;
V_fft = Y(:) ./ omega(:);

%% 5. Compute 1/3-Octave Band Velocity Levels
velocity_levels = zeros(size(freqs));

for i = 1:numel(freqs)
    % Define 1/3 octave band edges
    lower_freq = freqs(i) / (2^(1/6));
    upper_freq = freqs(i) * (2^(1/6));
    
    % Design bandpass filter for the current band
    d = designfilt('bandpassiir', 'FilterOrder', 4, ...
        'HalfPowerFrequency1', lower_freq, ...
        'HalfPowerFrequency2', upper_freq, ...
        'SampleRate', Fs);

    % Apply filter
    filtered_signal = filtfilt(d, acc_filtered);
    
    % Convert acceleration to velocity (integrate)
    velocity_signal = cumtrapz(t, filtered_signal); % Integration
    
    % Apply HP filter after integration (to remove DC  / drift mostly) 
    velocity_signal = filtfilt(d_hp, velocity_signal);

    % Get peak velocity for this band
    % remove the initial and end transients from observation 
    ind1 = round(0.2*N); % discard the first 20% samples 
    ind2 = round(0.8*N); % discard the last 20% samples
    velocity_levels(i) = max(abs(velocity_signal(ind1:ind2)));
end

%% 6. Plot FFT Velocity Spectrum
figure;
plot(f, V_fft, 'b');
xlabel('Frequency (Hz)');
ylabel('Velocity (mm/s)');
title('FFT-Based Velocity Spectrum');
grid on;
xlim([0, max(freqs) + 5]);

%% 7. Plot 1/3-Octave Band Velocity Spectrum
figure;
semilogx(freqs, velocity_levels, 'o-r'); % Log scale for better visibility
xlabel('Frequency (Hz)');
ylabel('Peak Velocity (mm/s)');
title('1/3-Octave Band Velocity Spectrum');
grid on;

[Mathieu NOE]

@Gerald

hello again

so finally, problem solved ?