disp(' ')
disp(' frf_from_th.m ')
disp(' ver 1.7  January 22, 2013 ')
disp(' by Tom Irvine  Email: tom@irvinemail.org ')
disp(' ')
disp(' This program calculates the frequency response function ')
disp(' and of a force and response time history pair. ')
disp(' It is best suited for an impulse force input.')
disp(' It processes the data as a single record. ');
%
close all;
%
clear THM;
clear t;
clear f;
clear amp;
clear length;
clear fff;
clear Z;
%
clear H1;
clear H2;
clear COH;
clear frf;
%
fig_num=1;
%
tpi=2.*pi;
%
disp(' ');
disp(' Select units:  1=English  2=metric ');
iu=input(' ');
%
disp(' ');
disp(' Select response amplitude type: ');
disp('  1=displacement ');
disp('  2=velocity ');
disp('  3=acceleration ');
iamp=input(' ');
%
disp(' ');
disp(' The input array must have three columns: ');
if(iu==1)
    if(iamp==1)
        disp('  time  force(lbf)  disp(in) ');
    end
    if(iamp==2)
        disp('  time  force(lbf)  vel(in/sec) ');
    end
    if(iamp==3)
        disp('  time  force(lbf)  accel(G) ');        
    end    
else
    if(iamp==1)
        disp('  time  force(lbf)  disp(m) ');
    end
    if(iamp==2)
        disp('  time  force(lbf)  vel(m/sec) ');
    end
    if(iamp==3)
        disp('  time  force(lbf)  accel(m/sec^2) ');        
    end      
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select input method: ');
disp('  1=Signals are in two separate arrays.  Each with time(sec) & amp ');
disp('  2=Signals are in a common array.  time(sec) & force & response ');
im=input(' ');
%
iflag=0;
%
if(im==1)
    disp(' ');
    disp(' Force ');
    disp(' ');
    [t1,A,dt1,sr1,tmx1,tmi1,n1,~]=enter_time_history();
    disp(' ');
    disp(' Response ');
    disp(' ');    
    [t2,B,dt2,sr2,tmx2,tmi2,n2,ncontinue]=enter_time_history();
%
    if(tmx1 ~= tmx2 || tmi1 ~= tmi2 || n1 ~= n2 )
       disp(' Error:  time columns must be the same ');
       iflag=1;
    else
       dt=dt1;
       sr=sr1;
       t=t1;
       n=n1;   
       tmx=tmx1;
       tmi=tmi1;
    end
%
else
    [t,A,B,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history_three();
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(ncontinue==1 && iflag ~=1)
%
    disp(' ');
%
    F=A;
    amp=B;
%
    nhalf=floor(n/2);
%
    dur=t(n)-t(1);
    df=1/dur;
    fmax=(1/dt);
%
%
%%
    for i=1:50
        if(n<2^i)
            nq=2^i;
            break;
        end
    end
%
    freq=linspace(0,fmax,nq);
    freq=fix_size(freq);
    omega=2*pi*freq;
    omega2=omega.^2;
%
    disp(' Calculate FFTs ');
%
    force_ft=fft(F,nq);  % potential zero-padding
    amp_ft=fft(amp,nq);
%
    dH1=(conj(force_ft)).*force_ft;
     H1=((conj(force_ft)).*amp_ft)./dH1;
%
    dH2=(conj(amp_ft)).*force_ft;
     H2=((conj(amp_ft)).*amp_ft)./dH2;    
%
    disp(' ');
    ipf=input(' Apply FRF peak filter?  1=yes  2=no ');
%
    disp(' ');
    if(ipf==1)
        disp(' Filter H1 ');
        [H1]=frf_median_filter(H1);
        disp(' Filter H2 ');    
        [H2]=frf_median_filter(H2);
    end
%
    disp(' Calculate Coherence ');
%    
    COH=zeros(nhalf,1);
    for i=1:nhalf
        COH(i)=abs(H1(i)/H2(i));
    end
%
    frf=H1;
%%
%
%  Impulse Response Function
%
    clear Y;
    clear amp_real;
    clear amp_imag;
    Y = ifft(frf,nq,'symmetric');
    if(iamp==3)
        Y(1)=0;
    end
    amp_real=real(Y);
    amp_imag=imag(Y);
%
%%
%
    clear fff;
    clear Z
    clear ph
%
    fff=freq(1:nhalf);
%
    Z =abs(frf(1:nhalf));
%
    ph =(180/pi)*atan2(imag(frf(1:nhalf)), real(frf(1:nhalf)));               
%%
%%%
    t=fix_size(t);
    t_real=fix_size(amp_real);
    t_imag=fix_size(amp_imag);
    fff=fix_size(fff);
%
    Z=fix_size(Z);
    ph=fix_size(ph);
%
%%%
%%
    disp(' ');
    fmin=input(' Enter plot lower frequency(Hz) ');
    fmax=input(' Enter plot upper frequency(Hz) ');
%
    clear temp_z;
    clear temp_f;
    clear temp_p;
    clear temp_c;
%
    for i=1:nhalf
        if(fff(i)>fmax)
            break;
        else
            temp_f(i)=fff(i);
            temp_z(i)=Z(i);
            temp_p(i)=ph(i);
            temp_c(i)=COH(i);
        end
    end
%
    clear fff;
    clear Z;
    clear ph;
    clear COH;
%
    fff=temp_f;
    Z=temp_z;
    ph=temp_p;
    COH=temp_c;
%
    figure(fig_num);
    fig_num=fig_num+1;
    subplot(3,1,1);
    plot(fff,ph);
%
    if(iamp==1)
        axis([fmin,fmax,-180,0])
        set(gca,'MinorGridLineStyle',':','GridLineStyle',':',...
                       'XScale','log','YScale','lin','ytick',[-180,-90,0]);    
    end
    if(iamp==2)
        axis([fmin,fmax,-180,180])
        set(gca,'MinorGridLineStyle',':','GridLineStyle',':',...
                         'XScale','log','YScale','lin','ytick',[-180,-90,0,90,180]);    
    end
    if(iamp==3)
        axis([fmin,fmax,0,180])
        set(gca,'MinorGridLineStyle',':','GridLineStyle',':',...
                         'XScale','log','YScale','lin','ytick',[0,90,180]);    
    end
%
    ylabel('Phase (deg)');
    if(iamp==1)
        title('Receptance Frequency Response Function');
    end
    if(iamp==2)
        title('Mobility Frequency Response Function');
    end
    if(iamp==3)
        title('Acceleration Frequency Response Function');         
    end
    grid on;
%
    subplot(3,1,[2 3]);
    plot(fff,Z);
%
    ymax=10^(ceil(log10(max(Z))));
%    
    ymin=ymax;
    clear length;
    for i=1:length(Z)
        if(Z(i)<ymin && fff(i)>fmin && fff(i)<fmax )
            ymin=Z(i);
        end
    end
%
    ymin=10^(floor(log10(ymin)));
%
    if(ymin<ymax/10000)
        ymin=ymax/10000;
    end
    axis([fmin,fmax,ymin,ymax]);
    set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale',...
                                                      'log','YScale','log');
    if(iamp==1)
        if iu==1
            ylabel('Disp/Force (in/lbf)');
        else
            ylabel('Disp/Force (m/N)');   
        end
    end
    if(iamp==2)
        if iu==1
            ylabel('Vel/Force (in/sec/lbf)');
        else
            ylabel('Vel/Force (m/sec/N)');   
        end
    end
    if(iamp==3)
        if iu==1
            ylabel('Accel/Force (G/lbf)');
        else
            ylabel('Accel/Force (m/sec^2/N)');   
        end
    end
    xlabel('Frequency (Hz)');
    grid on;
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(fff,COH);
    grid on;
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
         'YScale','lin');
    xlabel('Frequency(Hz)'); 
    ylabel('(\gamma_x_y)^2'); 
    title('Coherence');
    ymin=0.;
    ymax=1.2;
    axis([fmin,fmax,ymin,ymax]); 
%
    ijk=1;
    while(ijk==1)
        tmin=input(' Enter plot start time(sec) ');
        tmax=input(' Enter plot end time(sec) ');
        figure(fig_num);
        fig_num=fig_num+1;
        disp(' ');
        plot(t(1:n),t_real(1:n));
        ymin=min(ylim);
        ymax=max(ylim);
        axis([tmin,tmax,ymin,ymax]);
        if(iamp==1)
            title('Receptance Impulse Response Function');
        end
        if(iamp==2)
            title('Mobility Impulse Response Function');
        end
        if(iamp==3)
            title('Acceleration Impulse Response Function');         
        end
        grid on;
        if(iamp==1)
            if iu==1
                ylabel('Disp/Force (in/lbf)');
            else
                ylabel('Disp/Force (m/N)');   
            end
        end
        if(iamp==2)
            if iu==1
                ylabel('Vel/Force (in/sec/lbf)');
            else
                ylabel('Vel/Force (m/sec/N)');   
            end
        end
        if(iamp==3)
            if iu==1
                ylabel('Accel/Force (G/lbf)');
            else
                ylabel('Accel/Force (m/sec^2/N)');   
            end
        end
        xlabel('Time (sec)');
        grid on;
%
%%%
        disp(' ');
        disp(' Replot impulse response function with different limits?');
        disp('  1=yes  2=no ');
        ijk=input(' ');
    end
%
%
    disp(' ');
    disp(' Output arrays ');
    if(iamp==1)
            clear receptance_frf;
            clear receptance_irf;
            receptance_frf=[fff Z];
            receptance_irf=[t(1:n) t_real(1:n)];
            disp(' receptance_frf & receptance_irf ');
    end
    if(iamp==2)
            clear mobility_frf;
            clear mobility_irf;
            mobility_frf=[fff Z];
            mobility_irf  =[t(1:n) t_real(1:n)];
            disp(' mobility_frf & mobility_irf ');            
    end
    if(iamp==3)
            clear accelerance_frf;  
            clear accelerance_irf;
            accelerance_frf=[fff Z];
            accelerance_irf  =[t(1:n) t_real(1:n)];
            disp(' accelerance_frf & accelerance_irf ');             
    end
end