% Simulated accelerometer data from elevator in Taipei 101 office tower 
% (509 m, 101 floors).
% Accelerates to 16.8 m/s in 16s; reaches observation deck (91st floor) 
% after 37s from fifth floor.
% References: en.wikipedia.org/wiki/elevator, en.wikipedia.org/wiki/Taipei_101, 
% accessed 2014-Feb-04.

function [a,t,dt, a_true] = accel(rms_err)

% time array
dt = 0.02; % sampling period, seconds.
tmax = 37; % seconds, time at end of journey.
trest = 5; % seconds at rest before and after journey.
t = -trest:dt:tmax+trest;  % time in seconds.


% acceleration array
accel_interval = 16; % s
max_speed = 16.8;  % m/s
g = 9.81;  % acceleration of gravity, m/s^2
a = g * ones(size(t));
ss = find( (t<accel_interval) & (t>0) );
a(ss) = a(ss) + max_speed/accel_interval;
ss = find( (t>(tmax-accel_interval)) & (t<tmax) );
a(ss) = a(ss) - max_speed/accel_interval;

% smooth the acceleration profile
ramp_up = 0:dt:1;   % 1s up-ramp of triangular smoothing kernel
ramp_down = (1-dt):(-dt):0;  % 1s down-ramp of triangular smoothing kernel
filter = zeros(size(a)); 
filter(1:numel(ramp_down)) = ramp_down;
filter( (end+1-numel(ramp_up)):end ) = ramp_up;
filter_f = fft(filter);
filter_f = filter_f / filter_f(1);  % filter normalization
a_f = fft(a) .* filter_f.^3;  % 3 applications of the filter.
a_true = real( ifft(a_f) ); % the theoretical acceleration curve, before we apply noise.

% noise
% noise = normrnd(0, rms_err, size(a));
noise = rms_err * randn( size(a) );
a = a_true + noise;

% discretization
da = 0.03;
a = da * round(a/da);