%
disp('  ');
disp('  half_power_bandwidth.m  ver 1.1  April 29, 2013 ');
disp(' ');
disp('  by Tom Irvine ');
%
disp('  This method performs a curve-fit to determine the damping ratio ');
disp('  for a system excited by base excitation or an applied force. ');
disp(' ');
disp('  This method is equivalent to the half-power bandwidth method.');
disp('  ');
disp('  The input data must be a frequency response function. ');
disp('  ');
%
close all;
%
clear length;
clear f;
clear a;
clear fr;
clear ar;
clear omega;
%
fig_num=1;
%
disp(' ');
disp(' The response data array must have two columns:  freq (Hz) & amplitude ');
disp(' ');
%
disp(' Select data input method ')
disp('  1=data preloaded in matlab ')  
disp('  2=external ASCII file ');
disp('  3=Excel file ');
%
method=input(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(method==1)
    THF = input(' Enter the input filename:  ');
    sz=size(THF);
    f=double(THF(:,1));        
    a=double(THF(:,2));      
end
if(method==2)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
%      
    fid = fopen(filename,'r');
%        THF = fscanf(fid,'%g %g',[2 inf]);
%        THF=THF';
    THF=load(filename);
    sz=size(THF);
    f=double(THF(:,1));
    a=double(THF(:,2));
end
if(method==3)
    [filename, pathname] = uigetfile('*.*');
    xfile = fullfile(pathname, filename);
%        
    THF = xlsread(xfile);
%        
    f=double(THF(:,1));
    a=double(THF(:,2));  
end
%
n=length(f);    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Enter excitation type ');
disp(' 1=base acceleration  2=applied force ');
iex=input(' ');
%
disp(' ');
disp(' Enter response type ');
disp(' 1=acceleration  2=velocity  3=displacement ');
iresp=input(' ');
%
disp(' Enter response amplitude type ');
%
if(iresp==1)
  disp('  1= acceleration ');
  disp('  2= acceleration^2 ');
  irat=input(' ');
  if(irat==1)
      tlabel='Acceleration';
%    
      if(iex==1)
         MFRF=@(omega,rho,damp)(sqrt((1+(2*damp*rho)^2)/((1-rho^2)^2+(2*damp*rho)^2)));      
      else
         MFRF=@(omega,rho,damp)(sqrt(omega^2/((1-rho^2)^2+(2*damp*rho)^2)));             
      end
%      
  else
      tlabel='Acceleration^2';
%    
      if(iex==1)
         MFRF=@(omega,rho,damp)((1+(2*damp*rho)^2)/((1-rho^2)^2+(2*damp*rho)^2));        
      else
         MFRF=@(omega,rho,damp)((omega^2/((1-rho^2)^2+(2*damp*rho)^2)));          
      end
%       
  end
end
%
if(iresp==2)
  disp('  1= velocity ');
  disp('  2= velocity^2 ');
  irat=input(' ');
  if(irat==1)
      tlabel='Velocity';
%    
      if(iex==1)
         MFRF=@(omega,rho,damp)(sqrt(((1+(2*damp*rho)^2)/omega)/((1-rho^2)^2+(2*damp*rho)^2)));       
      else
         MFRF=@(omega,rho,damp)(sqrt(omega/((1-rho^2)^2+(2*damp*rho)^2)));              
      end
%      
  else
      tlabel='Velocity^2';
%    
      if(iex==1)
         MFRF=@(omega,rho,damp)(((1+(2*damp*rho)^2)/omega)/((1-rho^2)^2+(2*damp*rho)^2));       
      else
         MFRF=@(omega,rho,damp)((omega/((1-rho^2)^2+(2*damp*rho)^2)));          
      end
%       
  end   
end
%
if(iresp==3)
  disp('  1= displacement ');
  disp('  2= displacement^2 ');
  irat=input(' ');
  if(irat==1)
      tlabel='Displacement';
 %    
      if(iex==1)
         MFRF=@(omega,rho,damp)(sqrt(((1+(2*damp*rho)^2)/omega^2)/((1-rho^2)^2+(2*damp*rho)^2)));        
      else
         MFRF=@(omega,rho,damp)(sqrt(1/((1-rho^2)^2+(2*damp*rho)^2)));            
      end
%      
  else
      tlabel='Displacement^2';
%    
      if(iex==1)
         MFRF=@(omega,rho,damp)(((1+(2*damp*rho)^2)/omega^2)/((1-rho^2)^2+(2*damp*rho)^2));        
      else
         MFRF=@(omega,rho,damp)((1/((1-rho^2)^2+(2*damp*rho)^2)));            
      end
%       
  end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
plot(f,a);
grid on;
xlabel('Frequency (Hz)');
%
out1=sprintf('%s Response Function',tlabel);
title(out1);
%
fig_num=fig_num+1;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
fstart=input(' Enter starting frequency(Hz) for curve-fit ');
  fend=input(' Enter   ending frequency(Hz) for curve-fit ');
%
if(fstart<min(f))
    fstart=min(f);
end
if(fend>max(f))
    fend=max(f);
end
%
disp(' ');
nt=input(' Enter number of trials (suggest >= 10000)  ');
%
for i=1:n
    if(f(i)>=fstart)
        i1=i;
        break;
    end
end
%
for i=1:n
    if(f(i)>fend)
        i2=i-1;
        break;
    end
end
%
maxa=0;
%
for i=i1:i2
    if(a(i)>maxa)
        maxa=a(i);
        maxf=f(i);
    end
end
%
tpi=2*pi;
%
df=fend-fstart;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
error_max=1.0e+90;
%
omega=tpi*f;
%
Ar=0;
Axr=0;
fnr=0;
dampr=0;
%
disp(' ');
disp('  i     A     fn     damp ');
for i=1:nt
%    
    if(i<5 || rand()<0.5)
        A=(0.99+0.02*rand());
        damp=0.5*(rand())^2;
        fn=maxf;
%
    else
        if(rand()<0.4)
              fn=fnr*(0.995+0.01*rand());  
          damp=dampr*(0.99+0.02*rand());       
                A=Ar*(0.995+0.01*rand());
        else
            if(rand()<0.5)
                   fn=fnr*(0.995+0.01*rand());
            else
               damp=dampr*(0.995+0.01*rand());                  
            end
        end
    end
%
    if(fn<fstart || fn>fend)
        fn=maxf;
    end
%
    if(A<0.2 || A>5)
        A=maxa*(0.8+0.4*rand());
    end
%
    if(damp>0.5)
        damp=0.5*(rand())^2;        
    end
%
    error=0;
%
    omegan=tpi*fn;
%
    c=zeros(i2,1);
    for j=i1:i2
        omega=tpi*f(j);
        rho=omega/omegan;
        c(j)=MFRF(omega,rho,damp);
    end
%    
    AX=A*maxa/max(c);
    c=AX*c;
%
    for j=i1:i2
        bbb=abs(log10(a(j)/c(j)));
        error=error+bbb;
    end
%
    if(error<error_max)
%        
        error_max=error;
             dampr=damp;
              fnr=fn;
               Ar=A;
              Axr=AX;
              cr=c;
%
        out1=sprintf(' %d  %8.4g  %8.4g  %8.4g',i,A,fn,damp);
        disp(out1);
%
        k=1;
%
        for j=i1:i2
%        
            fr(k)=f(j);
            ar(k)=cr(j);
            k=k+1;
%
        end
%
        figure(fig_num);
        plot(f(i1:i2),a(i1:i2),fr,ar);
        grid on;
        legend('Measured Data','Curve-fit',1);
        xlabel('Frequency (Hz)');
        out1=sprintf('%s Response Function',tlabel);
        title(out1);
%
    end
%
end
%
Qr=1/(2*dampr);
%
disp(' ');
disp('         Results ');
disp(' ');
out1=sprintf('            fn =%8.4g Hz \n damping ratio = %8.4g \n             Q =%8.4g \n',fnr,dampr,Qr);
disp(out1);
%