MATLAB Examples

## 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.