% Script for producing a Fourier transform of vehicle acceleration
% over a test road for final project.
% Chris Gerdes
% ME 106/227
% Spring 2002
% 5/24/2002 

% Do a 10 second simulation with a sampling rate of 1 kHz

t = 0:0.001:10;

% Load in the road profile for our road.  This corresponds to
% the time vector above and a speed of 20 m/s.  So the vector
% is just the height encountered every millisecond by a vehicle
% moving over the road at 20 m/s.

load zroad

% Set the corresponding points along the road for the front axle

zf = zroad(1001:11001);

% Offset the rear axle by a time equal to L/V which is just
% offsetting the vector by (L/V)*1000 steps.
% For this test, V = 20m/s and we'll take L = 2.7m from
% the data in Assignment 7.  Thus we want 1000*(2.7/20) or
% 135 points offset.  You should change this for the wheelbase
% you have in your vehicle for the final project

zr = zroad(1001-135:11001-135);

% Set your input vector equal to the road inputs for lsim

u = [zf; zr];

% Here is where you stick in your 3 mass vehicle model and 
% simulate with appropriate A, B and C matrices to get the
% pitch and bounce accelerations and the driver's seat 
% vertical acceleration.  To show a nice Fourier transform
% of your data, take the following steps with your
% acceleration vectors.  The example assumes you have an
% acceleration vector of length 10001 called 'acc'.

% Set up the frequency axis for plotting the data.  This 
% should go up to half of the sampling frequency and have
% half as many points as you have time points minus 1 (we
% will drop the last time point when we take the fft.

freq = (0:4999)*(500/5000);

% Take the Fourier Transform of the acceleration data
% and use 10000 points.

accfft = fft(acc,10000);

% Scale the magnitude of the FFT appropriately based on our
% class discussion.  This gives it some physical meaning,
% although our main objective is to be able to compare two
% different designs.  Drop the second half of the points
% which are simply a mirror of the first.

Magaccfft = 2*abs(accfft(1:5000))/5000;

% Plot the Fourier transform magnitude.  You should do this
% for the sample Mercedes-like parameters and for your
% design and compare. 

plot(freq,Magaccfft);