%        ***** DEMONSTRATION OF EXPONENTIAL WINDOW METHOD  ****
%                 1-dof system subjected to step load

pi2=2*pi;		% 6.28318...

% *** FFT data
nn=128;			% number of points in time domain
nf=nn/2+1;		% number of points in frequency domain (0 to Nyquist)
tp=5;			% period of response
dt=tp/nn;		% time step
df=1/tp;		% frequency step, in Hz
dw=pi2*df;		% frequency step, in rad/s
nyq=0.5/dt;		% nyquist frequency, in Hz

% ****** Oscillator data
fn=1.2;         	% natural frequency, in Hz
wn=pi2*fn;		%    "        "    , in rad/s
mass=5.;		% oscillator mass
stiff=mass*wn^2;	% spring stiffness
flex=1/stiff;		% flexibility or compliance
d=0.01;			% fraction of damping
p=1000;			% amplitude of load

% *** Define forcing function
time=(0:dt:tp-dt);	% time axis vector
freq=(0:df:nyq);	% frequency axis vector (in Hz)
om=(0:dw:pi2*nyq);	%     "       "    "    (in rad/s)
force=p*(time>1);	% load amplitude * step function = p*step(t-1)
plot(time,force)	% plot forcing function
title('Forcing function')
xlabel('Time (s)')
grid on
pause

% ******     a) Solution using standard FFT method   *******

r=om/wn;		% tuning ratio
r2=r.^2;		% square of tuning ratio
h=flex./(1-r2+2*i*d*r);	% transfer function
F=fft(force*dt);	% FFT of forcing function
plot(freq,abs(h));	% plot transfer function, absolute value
title('Standard transfer function')
xlabel('Frequency (Hz)')
grid on
pause;
plot(freq,abs(F(1:nf)));% plot FT of forcing function
title('Fourier transform of forcing function')
xlabel('Frequency (Hz)')
grid on
pause;
H=[h,conj(h(nf-1:-1:2))];	% creating mirror image for FFT 
u=real(ifft(F.*H*df));	% Inverse FFT
plot(time,u)		% plot response
title('Standard (periodic) response function')
xlabel('Time (s)')
grid on
pause

% ******     b) Solution using exponential window method   *******
a=dw;			% imaginary component of frequency
E=exp(-time*a);		% decaying exponential window
E1=1./E;		% rising exponential window
force=force.*E;		% apply window to forcing function
r=(om-i*a)/wn;		% complex tuning ratio

r2=r.^2;		% square of complex tuning ratio
h=flex./(1-r2+2*i*d*r);	% transfer function
F=fft(force*dt);	% FFT of forcing function
plot(freq,abs(h))	% plot transfer function, absolute value
title('EWM Transfer function')
xlabel('Frequency (Hz)')
grid on
pause
plot(freq,abs(F(1:nf))); % plot FT of forcing function
title('EWM Fourier Transform of forcing function')
xlabel('Frequency (Hz)')
grid on
pause;
H=[h,conj(h(nf-1:-1:2))];	% creating mirror image for FFT 
u=real(ifft(F.*H*df));	% Inverse FFT
u=E1.*u;		% apply rising exponential window
plot(time,u)		% plot response
title('EWM response function')
xlabel('Time (s)')
grid on