function identifycellular(filename);
% IDENTIFYCELLULAR: Identify cellular touchtone phone number.
%
%   IDENTIFYCELLULAR("soundfile.au") loads the recording of a
%   cellular touchtone number, and attempts to identify the
%   telephone number by correlating against known cellular
%   tones.

% Author:
%   Alexander Perlis
%   SouthWest Regional Institute in the Mathematical Sciences
%   c/o Department of Mathematics
%   The University of Arizona
%   Tucson, AZ 85721
%   (520) 621-2137
%   aprl@math.arizona.edu
%
% Obtain the latest version of this software from:
%   http://www.math.arizona.edu/~rims/workshops/fourier/softwareMATLAB
%
% Bugs/problems/ideas:
%   None.
%
% Change history:
%   May 1997  Original version
%   Dec 1997  Modified to also work with Linux on Pentium with sound card.

%----------------------------------------------------------------------------


%% This program attempts to identify the cellular phone number
%% stored as an .au file in filename.

y = auread(filename);

%% Load prerecorded cellular tones
currentpath = pwd;
cd ~/matlab/cellular;
T1 = auread('cell1.au');
T2 = auread('cell2.au');
T3 = auread('cell3.au');
T4 = auread('cell4.au');
T5 = auread('cell5.au');
T6 = auread('cell6.au');
T7 = auread('cell7.au');
T8 = auread('cell8.au');
T9 = auread('cell9.au');
Tstar = auread('cellstar.au');
T0 = auread('cell0.au');
Tpound = auread('cellpound.au');
cd(currentpath);

%% Now try to recognize:

y = y/max(abs(y));
figure;
plot(y);

windowsize = 200;
pw = 2*windowsize;  %% the width of the patch that moves across
h = patch([1,1+pw,1+pw,1],[-1,-1,1,1],'r');
set(h,'EraseMode','xor');
drawnow;

T1 = T1(1:windowsize);
T2 = T2(1:windowsize);
T3 = T3(1:windowsize);
T4 = T4(1:windowsize);
T5 = T5(1:windowsize);
T6 = T6(1:windowsize);
T7 = T7(1:windowsize);
T8 = T8(1:windowsize);
T9 = T9(1:windowsize);
Tstar = Tstar(1:windowsize);
T0 = T0(1:windowsize);
Tpound = Tpound(1:windowsize);

T1 = T1/norm(T1);
T2 = T2/norm(T2);
T3 = T3/norm(T3);
T4 = T4/norm(T4);
T5 = T5/norm(T5);
T6 = T6/norm(T6);
T7 = T7/norm(T7);
T8 = T8/norm(T8);
T9 = T9/norm(T9);
Tstar = Tstar/norm(Tstar);
T0 = T0/norm(T0);
Tpound = Tpound/norm(Tpound);

x = 1;
len = length(y) - windowsize;
%disp('Length is')
%disp(len)
grabthreshold=0.4;

while x < len,
  % Advance window until no more noise
  w = y(x:x+windowsize);
  while (x < len) & (max(abs(w)) <= grabthreshold),
    x = x + windowsize;
    set(h,'XData',[x,x+pw,x+pw,x]);drawnow;
%    disp('skipping, x is')
%    disp(x)
    if (x < len)
      w = y(x:x+windowsize);
    end
  end
  % Advance window until it covers tone
  while (x < len) & (max(abs(y(x))) <= grabthreshold),
    x = x + 1;
    set(h,'XData',[x,x+pw,x+pw,x]);drawnow;
%    disp('covering tone, x is')
%    disp(x)
  end
  %% Now try to correlate
  while x < len,
    w = y(x:x+windowsize);
    %% In case window has nothing in it, get out of this loop.
    %% Can't happen the first time through, but could happen
    %% if we don't correlate and move the window.
    if max(abs(w)) <= grabthreshold,
      break;
    end
    w = w/norm(w);
    corr = [max(abs(xcorr(w,T1))), max(abs(xcorr(w,T2))),...
          max(abs(xcorr(w,T3))), max(abs(xcorr(w,T4))),...
          max(abs(xcorr(w,T5))), max(abs(xcorr(w,T6))),...
          max(abs(xcorr(w,T7))), max(abs(xcorr(w,T8))),...
          max(abs(xcorr(w,T9))), max(abs(xcorr(w,Tstar))),...
          max(abs(xcorr(w,T0))), max(abs(xcorr(w,Tpound)))];
%    disp(corr)
    if max(corr) > .65,
%      disp('Partial correlation')
      corr = corr > .65;
      if sum(corr) == 1,
        %% Okay, a single guy correlated above 0.1, so print him
        digits = ['1';'2';'3';'4';'5';'6';'7';'8';'9';'*';'0';'#'];
        sound(w);
        disp(char(corr*digits))
        %% We're correlated, so move window off the tone
        w = y(x:x+windowsize);
        while (x < len) & (max(abs(w)) >= grabthreshold),
          x = x + windowsize;
          set(h,'XData',[x,x+pw,x+pw,x]);drawnow;
          if x < len,
            w = y(x:x+windowsize);
          end
        end;
        %% Now break out of TryToCorrelate while loop
        break;
      else
        disp('Partial correlation, but didn''t correlate.')
      end
    end
  %% We didn't correlate, so move window, try again
  %% To move past garbage, we skip more than 1, in fact 20 at a time.
  x = x + 20;
  set(h,'XData',[x,x+pw,x+pw,x]);drawnow;
%  disp('didnt correlate, try again, x is')
%  disp(x)
  end
%disp('Now at position: ')
%disp(x)
end

% [eof]