disp(' ')
disp(' blast.m  version 1.2    February 14, 2013 ')
disp(' By Tom Irvine   Email:  tom@vibrationdata.com')
disp(' ')
disp(' This program calculates a response time history from');
disp(' an input time history and a transfer function where');
disp(' the transfer function is a complex Fourier transform.');
%
close all hidden;
%
clear H;
clear frf;
clear impulse_resp;
clear t;
clear ttt;
clear y;
clear t_imp;
clear length;
clear cc;
%
fig_num=1;
%
disp(' ');
disp(' Enter the time history');
disp(' ');
[t,y,dt,sr,tmx,tmi,nt,ncontinue]=enter_time_history();
%
%
figure(fig_num);
fig_num=fig_num+1;
%
plot(t,y);
title('Time History');
grid on;
xlabel('Time(sec)');
%
tmin=t(1);
t=t-tmin;
%
disp(' ');
[H,df,nf,ncontinue]=enter_H_complex();
%
figure(fig_num);
fig_num=fig_num+1;
%
plot(H(:,1),abs(H(:,2)));
title('Transfer Function Magnitude');
grid on;
xlabel('Frequency(Hz)');
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
%
disp(' ');
disp(' Enter transfer type: ');
disp('  1=half-sided (extends to Nyquist freq)');
disp('  2=full-sided (extends to sample rate freq)');
%
itype=input(' ');
%
if(itype==1)
    frf=zeros(2*nf,2);
    frf(1:nf,:)=H(1:nf,:);
%
    for i=1:nf
        frf(i+nf,1)=(i+nf)*df;
    end
%
    aa=H(:,2);
    bb= flipud(aa);
%
    for i=1:nf
        frf(i+nf,2)=real(bb(i))-(1i)*imag(bb(i));
    end
%
    nf=2*nf;
%
else
    frf(:,2)=H(:,2);
end
%
frf_amp=frf(:,2);
%
frf_amp=flipud(frf_amp);
%
disp(' ');
disp(' inverse Fourier transform ');
[t_real,t_imag]=inverse_fourier_transform(frf_amp,nf);
%
impulse_resp=( t_real+(1i)*t_imag );
%
impulse_resp=real(impulse_resp);
%
impulse_dt=1/max(frf(:,1));
tmax=(nf-1)*impulse_dt;
%
t_imp=linspace(0,tmax,nf);
%
%
[ttt,impulse_resp]=linear_interpolation(t_imp,impulse_resp,dt);
impulse_resp(1)=0.;
%
nt_th=nt;
nt_imp=length(impulse_resp);
%
%
impulse_resp=impulse_resp/nt_imp;   % normalization factor applied
%
cc=zeros((nt_th+nt_imp-1),1);
%
%  Convolution
%
progressbar
for i=1:nt_th
    progressbar(i/nt_th);
    cc((1:nt_imp)+i-1)=cc((1:nt_imp)+i-1)+y(i)*impulse_resp(1:nt_imp);
end
progressbar(1);
%
ntc=length(cc);
%
tc=linspace(0,(ntc-1),ntc);
%
figure(fig_num);
fig_num=fig_num+1;
plot(ttt,impulse_resp);
title('Impulse Response Function');
grid on;
xlabel('Time(sec)');
grid on;
xlabel('Time(sec)')
%
figure(fig_num);
fig_num=fig_num+1;
plot(tc,cc);
title('Response');
grid on;
xlabel('Time(sec)');
%
sz=size(cc);
if(sz(2)>sz(1))
    cc=cc';
end
%
sz=size(tc);
if(sz(2)>sz(1))
    tc=tc';
end
response=[tc cc];
%
disp(' Output file: ');
disp('   response ');
out1=sprintf('\n max=%8.4g   min=%8.4g ',max(cc),min(cc));
disp(out1);