clear all;

%%%%%%%%%%  DEFINE INPUT PARAMETERS  %%%%%%%%%%%%%%%%%%%%%%%
vel_kmPhr = 40; %desired velocity of car km/hr %1 mile = 1.609344 km
mass = 1500; %kg -> increasing mass increases tyre friction
% accel = 0.5; %km/s^2 :NOT IMPLEMENTED YET
CSarea = 2; %meter^2 front cross-sectional area of car -> increases air-borne sound
TyreRad = 0.35; %meters radius of tyre->reduces rpm/velocity ratio
Nsecs = 4; %No of seconds, the window length of the sound
movingRtoL = 1; %1 ->car is moving Right to Left, 0 ->car moving Left to Right
%%%%%%%%%%%% DO NOT CHANGE ANYTHING BEYOND THIS POINT%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% DEFINE CONSTANTS %%%%%%%%%%%%%%%%%%%
vel_sound = 343;%meter/sec: velocity of sound in the air
dist_ears = 0.25;% meters :distance between two ears
Fs=100000;             % sampling rate
Max_RPM = 8000;
vel = vel_kmPhr/3.6;% in units of meter/sec. N meter/sec = N*3.6 km/hour
gearRatio=8;%Chosen to make some sense out of velocity/rpm numbers
% max_vel = (rpm/60)*2*pi*TyreRad/gearRatio; % N meter/sec = N*3.6 km/hour
rpm = vel*gearRatio*60/(2*pi*TyreRad);
if rpm > Max_RPM
    rpm = Max_RPM;
    vel = (rpm/60)*2*pi*TyreRad/gearRatio;
    vel_kmPhr = vel*3.6;
    disp(['The lazy car refuses to run faster than ', num2str(vel_kmPhr), ' KiloMeters/hour']);
    disp(['...so the speed will be set to ', num2str(vel_kmPhr), ' KiloMeters/hour']);
end
% Air resistance is approximated by the following formula (Fluid Mechanics by Landau and Lifshitz, [Beckham] chapter 6, [Zuvich])
%Fair =  0.5 * Cd * A * rho * v2 where Cd = 0.3(empirically), rho = 1.29.
Kair = 0.5*0.3*CSarea;
Ktyre = 30*Kair*mass/1500;%this relation of 30 was derived for mass = 1500 kg.
Max_s1 = 5; %maximum value of noise in stroke1 as a percentage to the overall maximum
Max_s2 = 1; %stroke 2 is compression, so not much sound generated
Max_s3 = 25; %explosion
Max_s4 = 20; %exhaust

Engine_resonance = 20; %Hz %Arbitrarily chosen, 
EngMass = 10;%Engines have very high stiffness/mass ratio - this is just the mass of the piston+cylinder
EngDamp = 40; %This should be independent of resonance frequency
EngSpring = (Engine_resonance*2*pi)^2*EngMass;
Body_resonance = 2;
BodyMass = 50;
BodyDamp = 30;
BodySpring = (Body_resonance*2*pi)^2*BodyMass;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%1 rpm = 60 rps = 30 EngCycles/sec
Tstroke = 30/(4*rpm); %period of one stroke
t = 0:1/Fs:Tstroke;
%Make a kernel with shape of positive-half sinusoid
sine_kernel = sin(pi*(1/Tstroke)*t)';
s1 = Max_s1*sine_kernel.*(rand(length(t),1)- 0.5);
s2 = Max_s2*sine_kernel.*(rand(length(t),1)- 0.5);

s3core = [1; 0.01*s2(1:end-1)];
% % s3norm = lsim(EngineMechLP,s3core,[]);
% % EngineMechLP_num = fir1(20,Engine_resonance*2/Fs);%DOES NOT WORK
% EngineMechLP_num = butter(3,Engine_resonance*2/Fs);%DOES NOT WORK
% s3norm = filter(EngineMechLP_num, [1], s3core); 
Engsys = tf([0.1],[EngMass, EngDamp, EngSpring]);
% figure(100); [h w] = freqs([1],[EngMass EngDamp EngSpring],[0:1:20000]*2*pi); loglog(w/(2*pi), abs(h));
% title('Engine mechanical freq response');
EngsD = c2d(Engsys, 1/Fs,'zoh');%discretize the model
s3norm = lsim(EngsD, s3core, []);
% s3norm = impulse(Engsys, t);
s3 = (Max_s3/max(s3norm))*s3norm; 
% figure(101); plot(s3core);

s4_kernel = [length(t)-1:-1:0]'*1/length(t);%ramp kernel
% s4_kernel = [length(t)-1:-1:0]'*1/length(t);%exponetial decay
s4 = Max_s4*s4_kernel.*(rand(length(t),1)-0.5);

oneCycle = [s1;s2;s3;s4];
% figure(105);plot(oneCycle,'.-');title('One cycle of Engine before application of muffler'); 
testH = repmat(oneCycle,round(rpm*Nsecs/30),1);
OverallLP_num = butter(30, Body_resonance*2/Fs);%The body that attenuates the sound before it reaches the muffler
MufflerIN = filter(OverallLP_num,[1],testH);

% Hmuff = tf([EngMass, 0.2*EngDamp, EngSpring],[EngMass, EngDamp, EngSpring]);%
rpm40kmph = (40/3.6)*gearRatio*60/(2*pi*TyreRad);
%design muffler with the central attenuation frequency equal to that of the sound generated at 40 Km/sec velocity
Hmuff_wn = 2*pi*rpm40kmph/2;
Hmuff_zeta1 = 0.1; Hmuff_zeta2 = 1; %zeta2/zeta1 = 10
Hmuff = tf([1, 2*Hmuff_zeta1*Hmuff_wn, Hmuff_wn^2], [1, 2*Hmuff_zeta2*Hmuff_wn, Hmuff_wn^2]);
% figure(107); bode(Hmuff); title('Muffler response');
%muffD = c2d(Hmuff, 1/Fs,'zoh');%discretize the model:WHY? NO NEED
% sEngine = lsim(muffD, MufflerIN, []);%Sound from an engine
Tmuff   = [1:1:length(MufflerIN)]/Fs;
sEngine = lsim(Hmuff, MufflerIN, Tmuff);%Sound of engine out of the muffler

NormFactor_Air = 0.000001; %Factors to transform from force/friction to sound amplitude
NormFactor_Tyre = NormFactor_Air*0.0005; % chosen by trial and error method.

sAir = Kair*(vel*NormFactor_Air)^2*(rand(length(sEngine),1)- 0.5); %Air-borne sound is proportional to square of velocity
sTyre = Ktyre*vel*NormFactor_Tyre*(rand(length(sEngine),1)- 0.5);

% " The loudest component of the noise of passenger cars travelling over 40 km/hr is car tyre noise, and nobody knows why. 
% (Knowing why could be lucrative with 100 million car tyres being sold annually.) 
% Air pumping noise generation mechanisms are an active area of research, 
% and might include noise generated by partial singularities present due to the opening and closing
% of cavities as the tyre flexes." -  http://research.imb.uq.edu.au/~m.gagen/research/vel_sing.html
sstat = sEngine+ sAir+ sTyre; %stationary sound
%Now apply doppler effect on two halves of the sound 
% fD = vel_sound/(vel_sound - vel_source) OR fD = vel_sound/(vel_sound + vel_source)
stat1 = sstat(1:length(sstat)/2); % we will consider that vehicle is moving closer for this period
stat2 = sstat(length(stat1) + 1: end); % vehicle is moving farther

fD1 = vel_sound/(vel_sound - vel); 
fD2 = vel_sound/(vel_sound + vel);
fstat1 = fft(stat1);
fstat2 = fft(stat2);
F_step =  (Fs/2)*(1/(ceil((length(stat1)+1)/2)-1));%same for both halves
freq = 0:F_step:Fs/2;
dfreq1 = 0:F_step:ceil(Fs*fD1/2);
dfreq2 = 0:F_step:ceil(Fs*fD2/2);
dfst1 = zeros(length(dfreq1),1);
dfst2 = zeros(length(dfreq2),1);
for i = 1:length(freq)
    dfst1(round(i*fD1)) = fstat1(i);
    dfst2(round(i*fD2)) = fstat2(i);
end
idfst1 = real(ifft(dfst1(1:length(dfreq1))));
idfst2 = real(ifft(dfst2(1:length(dfreq2))));
testALL = [idfst1; idfst2];
% Now compute delay and attenuation to make the sound stereo
Ttotal = length(testALL)/Fs;
Dtotal = vel*Ttotal;
T1 = [0:1:length(idfst1)-1]/Fs;
T2 = [0:1:length(idfst2)-1]/Fs;
D1 = vel*T1(end) - vel*T1;
D2 = vel*T2;
min_dist = 0.2; %meters. If the distance is zero, sound = infinity.
%assuming motion from Left to right, if this is wrong the ears will be switched later 
dL1 = max(abs(D1-dist_ears/2),min_dist); dR1 = max(D1+dist_ears/2,min_dist);
% %     attL1 = D1./(D1-dL1); attR1 = D1./(D1-dR1);
    attL1 = 1./dL1; attR1 = 1./dR1;
    delL1 = floor(dL1*Fs/vel_sound); delR1 = floor(dR1*Fs/vel_sound);
    indx = 1:1:length(T1);
    earL1 = zeros(length(T1) + max(max(delL1,delR1)),1);
    earR1 = zeros(length(earL1),1);
    for i = 1:length(T1)
        earL1(indx(i) + delL1(i)) = earL1(indx(i) + delL1(i))+ idfst1(i)*attL1(i);
        earR1(indx(i) + delR1(i)) = earR1(indx(i) + delR1(i))+ idfst1(i)*attR1(i);
    end
    
    dL2 = max(D2+dist_ears/2,min_dist); dR2 = max(abs(D2-dist_ears/2),min_dist);
% %     attL2 = D2./(D2-dL2); attR2 = D2./(D2-dR2);
    attL2 = 1./dL2; attR2 = 1./dR2;
    delL2 = floor(dL2*Fs/vel_sound); delR2 = floor(dR2*Fs/vel_sound);
    indx = 1:1:length(T2);
    earL2 = zeros(length(T2) + max(max(delL2,delR2)),1);
    earR2 = zeros(length(earL2),1);
    for i = 1:length(T2)
        earL2(indx(i) + delL2(i)) = earL2(indx(i) + delL2(i))+ idfst2(i)*attL2(i);
        earR2(indx(i) + delR2(i)) = earR2(indx(i) + delR2(i))+ idfst2(i)*attR2(i);
    end 
    
earL = [nonzeros(earL1);(earL2)];
earR = [nonzeros(earR1);(earR2)];

%If motion is from right to left, swap the ear-signals
if (movingRtoL)
    temp = earL;
    earL = earR;
    earR = temp;
end
% Now attenuate all frequencies above 20KHz.
% LP20K = fir1(10, 20000*2/Fs, hann(11));
LP20K = butter(10, 20000*2/Fs);
testL = filter(LP20K,[1],earL);
testR = filter(LP20K,[1],earR);
maxFact = max(max(abs(testL)), max(abs(testR)));
testL = testL/maxFact;
testR = testR/maxFact;

if length(testL) >= length(testR)
    test = [testL(1:length(testR)), testR];
else
    test = [testL, testR(1:length(testL))];
end
%%%%%%%%%%% write wav file  %%%%%%%%%%%
wavname='test.wav';
wavwrite(test,Fs,wavname);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure(1)
% plot(test,'.-');
disp('Now playing...');
wavplay(test,Fs);
%Some ASSUMPTIONS:
%Tyres are never slipping
%Every motion is one dimensional
%Wind is always parallel and against car-motion (more valid for high velocities)

%Some NOTES: 
%car engines are made of aluminium alloys which have very high
%stiffness/weight ratio
%Wind Noise should be proportional to the frontal area, and velocity
%Tyre noise and wind are white noise sources with high amplitude at
%high velocities (they are not bandpass filters or even filters)
%Styre = K1*vel, Swind = K2*FrontArea*(vel)^2 
%Doppler effect does matter.

