disp(' ');
disp(' digitize_SPL.m    ver 1.2   May 21, 2009');
disp(' ');
disp(' By Tom Irvine   Email: tomirvine@aol.com');
%
%  This program is a modified version of a similar program written by 
%  Kent Hardy and Ryan Edwards.
%
%  This function allows the user to pick points off an image file.  The
%  intended use is to be able to scan a hard copy of an SPL plot. 
%  
%  Variables:
%    "image_file" is an image file which you have already loaded (perhaps
%  with the "imread" function).  The image must be upright when viewed with
%  the "imagesc" function, so you may need to rotate the data somewhat.  I
%  recommend that you try the function "rot90" to accomplish this.
%  
%
%  PROGRAM USE:
%    When you execute the function, the program will prompt you to pick
%  the lower left corner of the data plot, and then enter the frequency
%  and SRS value at that point.  Then, you will pick the lower right
%  corner of the data plot and input the frequency at that point.  Third,
%  you will pick the upper left corner of the data plot and input the
%  SRS magnitude at that point.  Finally, you will pick the points on the
%  visible SPL curve.  Remember to hit RETURN when you are done picking.
%  The program will then plot what you have chosen on top of the image so
%  you can see if you like how it fits together.
%    Once you have given the OK, the code converts the data points to
%  frequency/SPL format.  If you pick a point near the lower or upper
%  frequency bounds, the code will assumme you meant that the point should
%  be at the maximum/minimum frequency of the plot.  This should help
%  later when you try to put everything to a common frequency resolution.
%
clear bspl;
clear f;
clear fc;
clear freq;
clear B;
clear nn;
clear xx;
clear yy;
clear xxyy;
clear xxyy3;
clear srs_out;
clear spl_out;
clear figure(1);
%
disp(' ')
disp(' Select file input method ');
disp('   1=external image file ');
disp('   2=image preloaded into Matlab ');
file_choice = input('');
%
%% if(file_choice==1)
%%    disp(' Enter the input image filename: ');
%%    filename = input(' ','s');
%%    image_file = imread(filename);
%%else
%%    image_file = input(' Enter the matrix name:  ');
%%end
%%%%%%%%%%%%%%%%%%
if(file_choice==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
%      
    image_file = imread(filename);
%
else
    image_file = input(' Enter the image name:  ');
end
%%%%%%%%%%%%%%%%%%
%
  disp(' ');
% disp(' Select graph type ');
% disp(' 1=SRS  2=spl  3=other');
% itype = input(' ');
itype=3;
%
yesno='No';
while yesno(1)~='Y';
  clf;
  imagesc(image_file);
  axis image off; % Turn off the image axes
  set(gca,'Position',[0 0 1 1]);
  set(gcf,'colormap',[[1 0 0]; [0 0 0]]);
 
  zoommsg = {'Zoom in on the portion of the plot that you want to digitize.',...
        'Take care to keep the axes limits visible.',...
        'Press OK when you want to proceed'};
  zoom on;
  q0 = msgbox(zoommsg);
  waitfor(q0);
  
  q1=questdlg('Pick lower left corner of plot. ','','Proceed','Cancel','Proceed');
  if q1(1)=='P'; x1y1=ginput(1); else; error('Shutting Down'); end
%
  if(itype==1)
    txt2=cellstr(strvcat('Natural Frequency (Hz)','SRS (G)'));
  else
    txt2=cellstr(strvcat('Frequency (Hz)','SPL (dB)'));   
  end
%
  q1=inputdlg(txt2);
  q1a=char(q1);
  f1=str2double(q1a(1,:));
  s1=str2double(q1a(2,:));
  %
  q1=questdlg('Pick lower right corner of plot. ','','Proceed','Cancel','Proceed');
  if q1(1)=='P'; x2y2=ginput(1); else; error('Shutting Down'); end
% 
  if(itype==1)
     txt2=cellstr('Natural Frequency (Hz)');
  else
     txt2=cellstr('Frequency (Hz)');     
  end  
%
  q1=inputdlg(txt2);
  q1a=char(q1);
  f2=str2double(q1a(1,:));
  %
  q1=questdlg('Pick upper left corner of plot. ','','Proceed','Cancel','Proceed');
  if q1(1)=='P'; x3y3=ginput(1); else; error('Shutting Down'); end
%
  if(itype==1)
%     txt2=cellstr('SRS (G)');
  else
     txt2=cellstr('SPL (dB)');   
  end
%
  q1=inputdlg(txt2);
  q1a=char(q1);
  s2=str2double(q1a(1,:));
  %
  %  Plot the lines which were chosen and give the user a chance to change.
  %
  x1=[x1y1(1) x2y2(1)]; y1=[x1y1(2) x2y2(2)];
  q1=line(x1,y1); set(q1,'color',[1 1 0]);
  x1=[x1y1(1) x3y3(1)]; y1=[x1y1(2) x3y3(2)];
  q1=line(x1,y1); set(q1,'color',[1 1 0]);
  yesno=questdlg('Are the axes acceptable?');
end
%
yesno='No';
while yesno(1)~='Y';
%
if(itype==1)
  txt=[
    'Pick the points of the SRS curve.'
    'Choose as many points as desired.'
    'Press RETURN when done choosing. '];
else
   txt=[
    'Pick the points of the SPL curve.'
    'Choose as many points as desired.'
    'Press RETURN when done choosing. '];   
end
%
  q1=questdlg(txt,'','Proceed','Cancel','Proceed');
  if q1(1)=='P'; [xx,yy]=ginput; else; error('Shutting Down'); end
  %
  %  Overlay chosen points on image.
  %
  numpts=length(xx);
  for k1=1:numpts-1
    x1=[xx(k1) xx(k1+1)];
    y1=[yy(k1) yy(k1+1)];
    q1=line(x1,y1);
    set(q1,'color',[0 1 1]);
  end
  yesno=questdlg('Are the points acceptable?');
  if strcmpi(yesno,'No')
      q1 = findall(gca,'Color',[0 1 1]);
      delete(q1);
  end
end
xxyy=[xx yy]';
%
%  Translate the data into actual SRS numbers assumming a
%  loglog plot.  Start by rotating the data to a directly
%  x-y coordinate system (this to account for the fact
%  that the scanner often rotates the page slightly).
%
thet1=atan((x1y1(2)-x2y2(2))/(x2y2(1)-x1y1(1)));
thet2=atan((x3y3(1)-x1y1(1))/(x3y3(2)-x1y1(2)));
theta=(thet1+thet2)/2;
%  Create coordinate transformation matrix.
trmat=[cos(theta) -sin(theta)
  sin(theta) cos(theta)];
%  Rotate the axes and the curve.
xxyy2=trmat*xxyy;
xxyy2=xxyy2';
x1y1a=trmat*x1y1';
x2y2a=trmat*x2y2';
x3y3a=trmat*x3y3';
%
%  Convert to log data.
%
freq_decs=log10(f2/f1);
srs_decs=log10(s2/s1);
freq_length=(x2y2a(1)-x1y1a(1))/freq_decs;
% srs_length=(x1y1a(2)-x3y3a(2))/srs_decs;
srs_length=(x1y1a(2)-x3y3a(2));
for k1=1:numpts
  xxyy3(k1,1)=f1*10^((xxyy2(k1,1)-x1y1a(1))/freq_length);
  %  If frequency point is near the upper or lower limit, then
  %  assumme that the user intended to pick a point at the limit.
  fred1=abs((xxyy3(k1,1)-f1)/f1);
  fred2=abs((xxyy3(k1,1)-f2)/f2);
  if fred1<0.02; xxyy3(k1,1)=f1; end
  if fred2<0.02; xxyy3(k1,1)=f2; end
  xxyy3(k1,2)=s1 + (s2-s1)*((-xxyy2(k1,2)+x1y1a(2))/srs_length);
end
%    
%  Output 
%
if(itype==1)
    srs_out=xxyy3;
    disp(' ')
    disp(' The output filename is "srs_out" ')
    nn=length(srs_out(:,1));
    for(i=1:nn-1)
        if(srs_out(i,1)>=srs_out(i+1,1))
            disp(' Frequency error. Sorting data to correct error.');
        end
    end
    B = sortrows(srs_out);
    srs_out=B;
else
    spl_out=xxyy3;
% 
    nn=length(spl_out(:,1));   
    for(i=1:nn-1)
        if(spl_out(i,1)>=spl_out(i+1,1))
            disp(' Frequency error.  Sorting data to correct error.');
        end
    end 
    B = sortrows(spl_out);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear ref;
ref =1;
%
clear f;
clear aspl;
clear summ;
summ=0.;
%
f=B(:,1);
aspl=B(:,2);
%
disp(' ');
disp(' Convert to one-third octave format');
disp(' ');
%
	fc(1)=10.;
	fc(2)=12.5; 
	fc(3)=16.;
	fc(4)=20.;
	fc(5)=25.;
	fc(6)=31.5;
	fc(7)=40.;
	fc(8)=50.;
    fc(9)=63.;
	fc(10)=80.;
	fc(11)=100.;
	fc(12)=125.;
	fc(13)=160.;
	fc(14)=200.;
	fc(15)=250.;
	fc(16)=315.;
    fc(17)=400.;
	fc(18)=500.;
	fc(19)=630.;
	fc(20)=800.;
	fc(21)=1000.;
	fc(22)=1250.;
	fc(23)=1600.;
	fc(24)=2000.;
    fc(25)=2500.;
	fc(26)=3150.;
	fc(27)=4000.;
	fc(28)=5000.;
	fc(29)=6300.;
	fc(30)=8000.;
	fc(31)=10000.;
	fc(32)=12500.;
	fc(33)=16000.;
	fc(34)=20000.;
	fc(35)=25000.;
	fc(36)=31500.;
%
    num=max(size(f));
    k=1;
	for(i=1:36)
%
		for(j=1:(num-1))
%
			if(fc(i)>f(j) && fc(i) < f(j+1) )
%
				x1=f(j);
				x2=f(j+1);
%                
				x=fc(i);
%
				y1=aspl(j);
				y2=aspl(j+1);
%
				LXL=log10(x2)-log10(x1);
%
				if(LXL <=1.0e-20)
					disp(' LXL error ');
                end
%
				c2=(log10(x)-log10(x1))/LXL;
				c1=1-c2;
%
				bspl(k)=c1*y1+c2*y2;
%
				db=bspl(k);
                freq(k)=fc(i);
                k=k+1;             
				ms = ref*10^(db/20);
				summ=summ+(ms^2);
%
				break;
            end
			if(fc(i)==f(j))
				bspl(k)=aspl(j);
				db=bspl(k);
                freq(k)=fc(i);
                k=k+1;
				ms = ref*10.^(db/20);
				summ=summ+(ms^2);
				break;
            end
     		if(fc(i)==f(j+1))
				bspl(k)=aspl(j+1);
				db=bspl(k);
                freq(k)=fc(i);
                k=k+1;
				ms = ref*10.^(db/20);
				summ=summ+(ms^2);
				break;
            end       
        end
    end
%
	summ=sqrt(summ);
%
	rms=20.*log10(summ/ref);
%
	out1=sprintf('\n\n  Overall Level = %8.4g dB\n',rms);
    disp(out1);
%
    spl_out=[freq' bspl']
    disp(' ')
    disp(' The output filename is "spl_out" ')  
%
figure(2);
plot(freq,bspl);
%
ylabel('SPL (dB)');
xlabel('Center Frequency (Hz)');
%
out7=sprintf('One-Third Octave SPL   OASPL=%7.4g dB  Ref = 20 micro Pa',rms);
title(out7);
%
grid on;
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','lin');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp('Output ASCII text file?')
choice=input(' 1=yes  2=no  ' );
disp(' ')
%
if choice == 1 
   disp(' Enter output filename ');
   TH_filename = input(' ','s');
%
   fid = fopen(TH_filename,'w');
end
%
if choice == 1 
%
    for(i=1:max(size(freq)))
        fprintf(fid,'%12.5g %12.5g \n',freq(i),bspl(i));
    end
%
   disp(' ');
   disp(' Output filename: ');
   out1=sprintf('  %s',TH_filename);
   disp(out1);
   fclose(fid);
end