% Shock Response Spectrum calculation (Absolute Acceleration)
% Written by Bob Light, June 1999
% Based on impulse Invariant digital recursive filtering 
% Ref. Tom Irvine Introduction to the Shock Response Spectrum 
% See http://www.vibrationdata.com/SRS.htm

% Example input & basic plotting added by Dave Corben, Philips

% calculation is based on time response
dt=0.001; tpulse=0.030; tmax=0.050;
npulse=tpulse/dt; nmax=tmax/dt+1;
t=0 : dt : tmax;
y  = zeros(1,nmax);

%y(1:npulse)= ones(1,npulse); % square wave input
y(1:npulse) = sin( pi*t(1:npulse)/tpulse ); % half sine input

% srs_filt.m
%
% expects time vector in t and data in y

% initialize some values
delta_t=abs(t(2)-t(1));
zeta=0.0;
fmin=20.0;
fmax=(1.0/delta_t)/10.0;
d_octaves=1/6;
n_octaves=log(fmax/fmin)/log(2);
k_octaves=2.0^d_octaves;
f_stop=ceil(n_octaves/d_octaves);

% calculate frequencies
f_n(1)=fmin;
w_n(1)=2.0*pi*fmin;
for k=2:f_stop,
	w_n(k)=w_n(k-1)*k_octaves;
	f_n(k)=f_n(k-1)*k_octaves;
end
w_d=w_n*sqrt(1.0-zeta^2);

% initialize filter coefficients
a1=2.0*exp(-1.0*zeta*w_n*delta_t).* cos(w_d*delta_t);
a2=-1.0*exp(-2.0*zeta*w_n*delta_t);
b1=2.0*zeta*w_n*delta_t;
b2=w_n*delta_t.*exp(-1.0*zeta*w_n*delta_t);
b2=b2.*((w_n./w_d)*(1.0-2.0*zeta^2).*sin(w_d*delta_t)-2.0*zeta.*cos(w_d*delta_t));

%calculate response
x_max=zeros(size(w_n));x_min=x_max;
%%resp=zeros(1,size(t));
resp=zeros(1,length(t));
for k=1:f_stop,
   A=[1,-a1(k),-a2(k)];
   B=[b1(k),b2(k),0];
	resp=filter(B,A,y);
	x_max(k)=max(resp);
	x_min(k)=min(resp);
end;
x_mm=[x_max;abs(x_min)];
x_mm=max(x_mm);

% plot
figure
plot(f_n,x_mm,'k-x')
axis([fmin fmax 0 3])
xlabel('frequency (Hz)'); ylabel('Acceleration')

% clean up
clear a1;clear a2;clear b1;clear b2;clear A;clear B;
clear k;clear zeta;clear fmin;clear fmax;clear resp;
clear ans;clear w_n;clear w_d;clear delta_t;clear f_stop;
clear d_octaves;clear k_octaves;clear n_octaves;
