MATLAB Examples

Contents

Example

We are loading acceleration data made of 288 time series assumed to be acceleration histories of a suspension bridge, recorded at mid-span. Only the natural frequency for the vertical DOF of the bridge is considered here.

Traffic induced vibration has been introduced every 25 samples, i.e. sample 25,50,...,275 near t = 300 sec.

Time series are supposed to be recorded between the 01/01/2016 at 00:00:00 and the 02/01/2016 at 23:50:00.

N.B.: Every data used in the present script is simulated.

clear all;close all;clc;
load data.mat
[Nsamples,N]=size(Az);

% to disable all the warnings here, because of the high number of samples
% used:
% warning('off')

Visualisation of acceleration data

The top figure shows wind induced vibrations only, whereas the bottom one shows a sudden peak indicating dominating traffic induced vibrations.

figure
subplot(211);
plot(t,Az(1,:))
ylabel('Az (m/s^2)')
title('wind induced vibrations only ')

subplot(212);
plot(t,Az(25,:))
ylabel('Az (m/s^2)')
xlabel('time (s)');
title('wind + traffic induced vibrations')
set(gcf,'color','w')

method 1: Days and Nights separation

To separate data dominated by wind from those dominated by traffic, we can separate diurnal and nocturnal data, by assuming that traffic-induced vibrations only occur at day time. Here day time is taken between 5 a.m. and 6 p.m This method is fast and efficient but too many samples are detected as "dominated by traffic".

% nocturnal and diurnal separation is taken as method = 1;
method = 1;


% optional paramters
param.tStart = '05:00:00'; % beginning of day time
param.tEnd = '18:00:00'; % end of day time

% Each sample is tested individually in a loop
y = zeros(1,Nsamples);
for ii=1:Nsamples,
param.time = time(ii);% mandatory parameter: information on time (datenum)
[y(ii)] = trafficDetection(Az(ii,:),method,param);
end
% find the sample detected as dominated by traffic
ind1 = find(y==0);
fprintf(['traffic introduced ',num2str(numel(ind0)),' times \n'])
fprintf(['traffic detected ',num2str(numel(ind1)),' times \n'])

% Illustration:
% green area are dominated by wind (nocturnal data)
% red area are dominated by traffic (diurnal data)
clf;close all;
figure
hold on
area(time,y,'FaceColor','g');
area(time,abs(y-1),'FaceColor','r');
legend('no traffic','traffic')
datetick('x','HH')
xlabel('time (HH)')
set(gcf,'color','w')
traffic introduced 11 times 
traffic detected 158 times 

method 2: fluctuating RMS of displacement

In this method (method = 2), the variation rate tau is evaluated, and broadly corresponds to the flucutating RMS of the bridge acceleration weighted by the wind velocity fluctuations. This method works pretty well, but depends on the initial threshold that has to be fixed. this threshold sepcifies the limit over which traffic is detected as non-negligible.

method = 2;
clear parameters


% optional paramters
% The threshold beyond which traffic is assumed to be dominating is fixed
param.Threshold = 0.7;
% number of subsegment to divide the accelerations data: 200 subsegments:
param.Nu = 200;



% Similarly to method 1, every sample is tested individually
y = zeros(1,Nsamples);
for ii=1:Nsamples,
param.u = u(ii,:); % the wind velocity flucutations is used as input
[y(ii)] = trafficDetection(Az(ii,:),method,param);
end

% check the number of sample detected as dominated by traffic
ind1 = find(y==0);
fprintf(['traffic introduced ',num2str(numel(ind0)),' times \n'])
fprintf(['traffic detected ',num2str(numel(ind1)),' times \n'])
disp([ind1(:)]);


% Illustration:
% top: no sudden jump of tau = only wind induced vibrations
% bottom: a jump of tau is detected near t = 300 sec, i.e. we have traffic
% induced vibrations

clf;close all;
T = reshape(t ,[],300);
U = reshape(u(15,:),[],300);
X = reshape(Az(15,:) ,[],300);
tau= [nan,abs(std(diff(X,1,2))./std(diff(U,1,2)))];
figure
plot(mean(T),tau,'ko','markerfaceColor','r');
set(gcf,'color','w')
xlabel('time (s)')
ylabel(' \tau ')
ylim([0,1])

U = reshape(u(25,:),[],300);
X = reshape(Az(25,:) ,[],300);
tau= [nan,abs(std(diff(X,1,2))./std(diff(U,1,2)))];
figure
hold on; box on;
plot(mean(T),tau,'ko','markerfaceColor','r');
% plot([300,300],[0,1],'k--'); % jump detected here
set(gcf,'color','w')
xlabel('time (s)')
ylabel(' \tau ')
ylim([0,1])
traffic introduced 11 times 
traffic detected 11 times 
    25
    50
    75
   100
   125
   150
   175
   200
   225
   250
   275

method 3: High frequency vs low frequency components

This method (method = 3) relies on the idea that traffic induced vibrations affects generally more the high frequency part of the response spectrum than the low frequency part. On the contrary, wind-induced vibrations are mainly concentrated on the low frequency part of the response spectrum. In the present exemple, the frequency threshold is taken as equal to 1Hz, i.e. below 1Hz we have the "low frequency domain" and over 1 Hz we have the "high frequency domain".

method = 3;
param.fs = 1/median(diff(t)); % sampling frequency
param.Threshold = 0.5; % frequency threshold

y = zeros(1,Nsamples);
for ii=1:Nsamples,
[y(ii)] = trafficDetection(Az(ii,:),method,param);
end
ind1 = find(y==0); % get samples detected as dominated by traffic

fprintf(['traffic introduced ',num2str(numel(ind0)),' times \n'])
fprintf(['traffic detected ',num2str(numel(ind1)),' times \n'])
fprintf('every sample with traffic has been detected: \n')
disp([ind1(:)]);

% Illustration of method 3:
% Low frequency part and high frequency part. The top figure shows the case
% of wind-induced vibrations only. We have a single peak located near 0.2
% Hz.The bottom picture shows multiple peaks at hgih frequencies that are
% introduced by traffic.


figure
title('Without traffic')
[PSD,freq]=pwelch(Az(15,:),1024,512,1024,param.fs);
indF = find(freq>=param.Threshold);
loglog(freq(indF),PSD(indF),'k',freq(1:indF(1)),PSD(1:indF(1)),'r')
legend('high frequency part','low frequency part','location','NorthWest')
axis tight
set(gcf,'color','w')
xlabel('freq (Hz)')
ylabel(' PSD (m^2s^{-1})')

figure
title('With traffic')
[PSD,freq]=pwelch(Az(25,:),1024,512,1024,param.fs);
indF = find(freq>=param.Threshold);
loglog(freq(indF),PSD(indF),'k',freq(1:indF(1)),PSD(1:indF(1)),'r')
legend('high frequency part','low frequency part','location','NorthWest')
axis tight
set(gcf,'color','w')
xlabel('freq (Hz)')
ylabel(' PSD (m^2s^{-1})')
traffic introduced 11 times 
traffic detected 11 times 
every sample with traffic has been detected: 
    25
    50
    75
   100
   125
   150
   175
   200
   225
   250
   275

Conclusions

Three simple methods to detect acceleration samples dominated by trafic induced vibrations are presented here. The first one is fast, but may be too conservative, and miss some important traffic induced vibrations. The second method works pretty well, but results depend on the threshold parameter that is not intuitive. The last method is easy to implement, works well too, and the threshold parameter is more intuitive.