function [tones,t_start,t_finish] = dtmf_det(x,D,M,Nw,win,fft_len,Nt,Np,Thresh,en_pt,tol,en_fi,en_di,f_name)
%tones = dtmf_det(x)
%return tones detected in vector x
%inputs:
%x        -> 8000 Hz audio vector read in from 16-bit linear .au file
%D        -> D sample-steps per window 
%M        -> M samples per window
%Nw       -> Nw consecutive windows to average via Welch method for psd
%win      -> window desired 0 -> rectangular, 1 -> Bartlett, 2 -> Hanning
%            3 -> Hamming
%fft_len  -> Length of fft to take
%Nt       -> Number of samples to define a valid tone detect length
%Np       -> Number of samples to seperate tone detect events
%Thresh   -> minimum threshold power to validate a detection
%en_pt    -> enable ploting of various stages of algorithm
%tol      -> tolerance on tone detect width
%en_fi    -> apply filter to data
%en_di    -> enable differential measurement for speech rejection
%f_name   -> input f_name that data tones were obtained from
%outputs:
%tones    -> tones detected from input audio vector as a string
%t_start  -> corresponding time in which tone was detected
%t_finish -> corresponding time in which tone stopped
% dtmf_det.m
% author: Nicholas Kottenstette
% date:   12/13/2004
% copyright: 2004 by Nicholas Kottenstette
  if en_pt
    fid = figure(2);
    set(gcf,'Name',f_name);
  end
  t0 = clock;%Start timing
  [m,n] = size(x);
  if n > m
    x = x';
  end
  %scale x to +/- 1 peak-to-peak if not doing differential detection
  %if ~en_di
  %  x = x*(1/(max(x)-min(x)));
  %end
  %rms_pow = mean(x.*x);
  %display(['rms_power = ', num2str(10*log10(rms_pow),3), ' dB']);
  Fs = 8000;
  dtl = 330;
  dth = 440;
  fl = [697 770 852 941];
  fh = [1209 1336 1477 1633];

  fv = [fl fh];           %tones to detect
 %fhv = fv*2;
  fhv = [fl/2 fh*2];      %second harmonics for high tones
                          %first harmonics for low tones
  
  fvmin = (1-tol)*fv;
  fvmax = (1+tol)*fv;
  fhvmin = (1-tol)*fhv;
  fhvmax = (1+tol)*fhv;
  Nfilt = 2^11;
  if en_fi
    status('Applying Tones bandpass',t0);
    [xf,h,G,fG]=bandpass(fvmin,fvmax,0,-40,x,Nfilt,Fs,-1,0);
    status('Tones Bandpass complete',t0);
    if en_di
      status('Applying harmonic bandpass',t0);
      [xhf,hh,Gh,fGh]=bandpass(fhvmin,fhvmax,0,-40,x,Nfilt,Fs,-1,0);
      status('Harmonic Bandpass complete',t0);
    end
    %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
    if en_pt
      NH = 2^ceil(log2(length(h))+4);
      H = dsp_fft(h,NH,[],[]);
      fH = [0:NH/2]*Fs/NH;
      lG = length(fG);
      subplot(2,2,1);
      plot(fH,20*log10(abs(H(1:1+NH/2))),fG(1:lG),20*log10(abs(G(1:lG))), ...
           '.');
      title([f_name(1:min(length(f_name),10)), ': Tone Bandpass Filter Response']);
      xlabel('frequency Hz');
      ylabel('H(f) dB');
      if en_di
        Hh = dsp_fft(hh,NH,[],[]);
        subplot(2,2,2);
        plot(fH,20*log10(abs(Hh(1:1+NH/2))),fGh(1:lG),20*log10(abs(Gh(1:lG))), ...
             '.');
        title([f_name(1:min(length(f_name),10)), ': Harmonics Bandpass Filter Response']);
        xlabel('frequency Hz');
        ylabel('Hh(f) dB');
      end
    end
    %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%
  else
    xf = x;
    if en_di
      xhf = x;
    end
  end
  w = wind(M,win);
  status('Calculating PSD via welch',t0);
  E  = calc_psd_welch(xf,D,M,w,fft_len,Nw);
  %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
  if en_pt
    subplot(2,2,3);
    plot_periodgram(E,[],Fs,D,false);
    title([f_name(1:min(length(f_name),10)), ': PSD of filtered tones via welch']);
  end
  %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%
  status('Done calculating PSD via welch',t0);
  if en_fi & en_di
    status('Calculating Harmonic PSD via welch',t0);
    Eh = calc_psd_welch(xhf,D,M,w,fft_len,Nw);
    
    %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
    if en_pt & en_di
      subplot(2,2,4);
      plot_periodgram(Eh,[],Fs,D,false);
      title('PSD of filtered harmonics via welch');
    end
    %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%
    status('Done calculating Harmonic PSD via welch',t0);
  else
    if en_di
      Eh = E; %Without filtering wont change result of fft
    end
            %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
    if en_pt & en_di
      subplot(2,2,4);
      plot_periodgram(Eh,[],Fs,D,false);
      title('PSD of filtered harmonics via welch');
    end
    %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%
  end
  status('Calculating Tone Frequency Peaks',t0);
  Ep = freq_peaks(E,fvmin,fvmax,Fs);
  %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
  if en_pt
    fid = figure(fid + 1);
    set(gcf,'Name',f_name);
    subplot(2,2,1);
    plot_periodgram(Ep,fvmin,Fs,D,true);
    title('PSD freq-peaks of tones');
  end
  %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%
  status('Done calculating Tone Frequency Peaks',t0);
  if en_di
    status('Calculating Harmonic Tone Frequency Peaks',t0);
    Eph = freq_peaks(Eh,fhvmin,fhvmax,Fs);
    status('Done calculating Harmonic Tone Frequency Peaks',t0);
  end
  %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
  if en_pt & en_di
    subplot(2,2,2);
    plot_periodgram(Eph,fhvmin,Fs,D,true);
    title('PSD freq-peaks of harmonics');
  end
  %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%
  if en_di
    Edif = Ep - Eph;
  else
    Edif = Ep;% - rms_pow;
  end
  Edif = max(Edif,0);
  %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
  if en_pt
    subplot(2,2,3);
    plot_periodgram(Edif,fvmin,Fs,D,true);
    %msg = sprintf('PSD diff. en-di = %d; rms-power = %2.3f
    %dB',en_di,10*log10(rms_pow));
    msg = sprintf('PSD diff. en-di = %d',en_di);
    title(msg);
  end
  %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%
  status('tone_detect',t0);
  [tones,t_start,t_finish]=tone_detect(Edif,Nt,Np,Thresh,en_pt,Fs,D,fv,f_name);
  status('tone_detect complete',t0);
  
function [tones,t_start,t_finish]=tone_detect(E,Nt,Np,Thresh,en_pt,Fs,D,fv,f_name)
  tones='';
  t_start=[];
  t_finish=[];
  TC = ['1' '2' '3' 'A';...
        '4' '5' '6' 'B';...
        '7' '8' '9' 'C';...
        '*' '0' '#' 'D'];
  %identify tones using ML detection
  %assume length of pause is as likely as length of tone
  %assume tones are uniformly distributed unless doing dial tone
  %detection once a valid tone other than a dial tone is detected
  %it is not possible to detect an additional dial tone
  %First we determine if any valid tone levels are present
  %a valid level is one which exceeds a given threshold for
  %Nt consecutive blocks after a valid pause length has been detected
  %once valid tones are detected we then take the most likely by
  %taken the one with the greatest power
  %the pause length is defined as the time between valid tones Np
  %the dial-tone detect does not requir a valid pause time
  dt = D/Fs;
  ni = 1;
  [me,ne] = size(E);
  
  Edl = E(1:4,1:ne);
  Edh = E(5:8,1:ne);
  
  Edl = Edl(:,ni:ne);
  Edh = Edh(:,ni:ne);
  Etl = zeros(size(Edl));
  Eth = zeros(size(Edh));
  %Apply threshold detection
  for i = 1:4
    for j = 1:4
      Idlh = find( (Edl(i,:) >= Thresh) & (Edh(j,:) >= Thresh) );
      %Determine index for all samples which are in Nt consecutive order
      for k = 1:Nt-1
        if length(Idlh) <= 1
          Idlh = [];
          break;
        end
        dIdlh = Idlh(2:length(Idlh)) - Idlh(1:length(Idlh)-1);
        %Get index for all consecutive threshold exceeding levels
        Idlh = Idlh(1 + find( dIdlh == 1 ));
      end
      %We now have indexes for valid low and high tone pairs (i,j) without
      %considering the pause time between tones but meets consecutive
      %threshold limits fill in temp E matrix with appropriate values
      if length(Idlh) ~= 0
        %Next need to restore the indexes which were "chomped off"

        tope = [Idlh(2:length(Idlh)),Idlh(length(Idlh))+2];
        bote = Idlh;
        topb = [Idlh(1)-2,Idlh(1:length(Idlh)-1)];
        botb = Idlh;
        ends = find( (tope-bote) > 1 );
        begs = find( (topb-botb) < -1 ); 
        Iends = Idlh(ends);
        Ibegs = max(Idlh(begs)-(Nt-1),1);
        nIdlh = zeros(1,length(Idlh) + length(ends)*(Nt-1));
        nee = 0;
        for ii = 1:length(Iends)
          inds = Ibegs(ii);
          inde = Iends(ii);
          dN = inde - inds;
          se = nee + 1;
          nee = se + dN;
          try,
            res = [Ibegs(ii):Iends(ii)];
            nIdlh(se:nee) = res;
          catch,
            display(se)
            display(nee)
            display(res)
            display(nIdlh)
            display(inds)
            display(inde)
            display(Idlh)
            display(dN)
            error('dtmf_det1:debug','%s',lasterr);
          end
        end        

        Idlhold = Idlh;
        Idlh = nIdlh;
        try,
          Etl(i,Idlh) = Edl(i,Idlh);
          Eth(j,Idlh) = Edh(j,Idlh);
        catch,
          display(Idlhold)
          display(nIdlh)
          error('dtmf_det1:debug','%s',lasterr);
        end
      end
    end
  end
  %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
  if en_pt
    subplot(2,2,4);
    plot_periodgram([Etl;Eth],fv,Fs,D,true);
    title([f_name(1:min(length(f_name),10)), ': Threshold Rejection PSD']);
  end
  %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%%%%%%
  %Update Edl and Edh with ML seperation
  Edl = ml_detect(Etl,Nt);
  Edh = ml_detect(Eth,Nt);
  [rl,cl]=find(Edl>0);
  [rh,ch]=find(Edh>0);
  %display(sum(cl-ch))
  %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%
  if en_pt
    figure(gcf + 1);
    set(gcf,'Name',f_name);
    subplot(2,2,1);stem(ch*dt,rh);
    ylabel('high-index');xlabel('time (seconds)');
    title([f_name(1:min(length(f_name),10)),': ml-detection for high-tones']);
    subplot(2,2,3);stem(cl*dt,rl);
    ylabel('low-index');xlabel('time (seconds)');
    title([f_name(1:min(length(f_name),10)), ': ml-detection for low-tones']);
  end
  %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%%%%%%
  %There is now only one unique column for every tone pair
  %Now just need to check for Np spacing to differentiate tones
  [m,n]=size(Edl);
  Esl = zeros(m,1);
  Esh = zeros(m,1);
  flag = false;
  if length(cl)
    ne = cl(length(cl));
    dn = n - ne;
    if dn > Np
      flag = true;
    end
    ns = cl(1);
  end
  for i = 2:length(cl)
    ci = cl(i);
    cip = cl(i-1);
    dn = ci-cip;
    if (dn > Np)
      il = find(Esl == max(Esl));
      ih = find(Esh == max(Esh));
      if((length(il)== 1) & (length(ih) == 1))
        tones    = [tones TC(il,ih)];
        t_start  = [t_start, ns*dt];
        t_finish = [t_finish, cip*dt];
        %very very dangerous for plotting only
        temp = Edl(il,ns:cip);
        Edl(:,ns:cip) = 0;
        Edl(il,ns:cip) = temp;
        temp = Edh(ih,ns:cip);
        Edh(:,ns:cip) = 0;
        Edh(ih,ns:cip) = temp;
        %very very dangerous for plotting only
      else
        %very very dangerous for plotting only
        Edl(:,ns:cip) = 0;
        Edh(:,ns:cip) = 0;
        %very very dangerous for plotting only
      end
      Esl = Edl(:,ci);
      Esh = Edh(:,ci);
      ns = ci;
    elseif dn == 1
      Esl = Esl + Edl(:,ci);
      Esh = Esh + Edh(:,ci);
    else %Possibly New tone other tone lost
      il = find(Esl == max(Esl));
      ih = find(Esh == max(Esh));
      inl = find(Edl(:,ci) == max(Edl(:,ci)));
      inh = find(Edh(:,ci) == max(Edh(:,ci)));
      if (inl ~= il) & (inh ~= ih)
        Esl = Edl(:,ci);
        Esh = Edh(:,ci);
        ns = ci;
      end
    end
  end
  if flag
    il = find(Edl(:,ns) == max(Edl(:,ns)));
    ih = find(Edh(:,ns) == max(Edh(:,ns)));
    tones = [tones, TC(il,ih)];
    t_start  = [t_start, ns*dt];
    t_finish = [t_finish,ne*dt]; 
  end
  %%%%%%%%%%%%%%%% Plotting %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if en_pt
    subplot(2,2,2);
    plot_periodgram(Edh,fv(5:8)*(1-.015),Fs,D,true);
    title([f_name(1:min(length(f_name),10)), ': Space filtering PSD High Freq']);
    subplot(2,2,4);
    plot_periodgram(Edl,fv(1:4)*(1-.015),Fs,D,true);
    title([f_name(1:min(length(f_name),10)), ': Space filtering PSD Low Freq']);
  end
  %%%%%%%%%%%%%%%% Done Plotting %%%%%%%%%%%%%%%%%%%%%%%
  
  
function Ep = freq_peaks(E,fmin,fmax,Fs)
  [m,n] = size(E);
  df = Fs/m;
  m_min = floor(fmin/df);
  m_max = ceil(fmax/df);
  Ep = zeros(length(fmin),n);
  for i = 1:length(fmin)
    Ep(i,:) = max(E(m_min(i):m_max(i),:),[],1);
  end
  
function E = calc_psd_welch(x,D,M,w,fft_len,Nw)
  [m,n]=size(x);
  if m < n
    x = x';
  end
  L = (length(x) - (M-1))/D + 1;
  L = ceil(L);
  %don't pad with zeros will cause wacky-detection with noisy data
  %pad with the last repeated frame of D samples
  lx = length(x);
  dlx = lx;
  lL = (L-1)*D + M -1 + 1;
  if lL > lx
    x =  [x;zeros(lL-lx,1)];
  end
  i = 1;
  off = 0;
  while(dlx < lL)
    xis = (i-1)*D + dlx + 1;
    xie = xis + D - 1;
    dlx = dlx + D;
    if (dlx > lL)
      off = dlx - lL;
      dlx = lL;
    end
    x(xis:xie-off) = x(xis-D:xie-D-off);
  end
  E = zeros(fft_len,L);
  U = mean(w.*w);
  MU = M*U;
  try,
    fp = load('fft_coef/coef.mat');
    c = fp.c;
    fp = load('fft_coef/rind.mat');
    rind = fp.rit{floor(log2(fft_len))};
  catch,
    c = [];
    rind = [];
    display('dtmf_det.m Warning unable to load coeficients and reverse index');
  end
  for i = 1:L
    xs = (i-1)*D + 1;
    xe = xs + M - 1;
    xv = x(xs:xe,1).*w;
    E(:,i) = (1/MU)*(abs(dsp_fft(xv,fft_len,c,rind)).^2);
  end
  sub = floor((Nw-1)/2);
  add = ceil((Nw-1)/2);
  for i = 1:L
    ms = max(1,i-sub);
    me = min(L,i+add);
    E(:,i) = mean(E(:,ms:me),2);
  end
  
function w = wind(Lw,win)
  n = [0:Lw-1]';
  switch win
   case 0 %Rectangular
    w = ones(Lw,1);
   case 1 %Bartlett
    w = 1 - 2*abs(n - (Lw-1)/2)/(Lw-1);
   case 2 %Hanning
    w = (1/2)*(1 - cos(2*pi*n/(Lw-1)));
   case 3 %Hamming
    w = 0.54 - 0.46*cos(2*pi*n/(Lw-1));
   otherwise %Rectangular
    w = ones(Lw,1);
  end

function status(msg,t0)
  display([msg ' ' num2str(etime(clock,t0)) ' seconds elapsed since dtmf_det1 called']);
  
function plot_periodgram(E,f,Fs,D,sp)
  [m,n]=size(E);
  t=[0:n-1]*(D/Fs);
  if length(f) == 0
    df = Fs/m;
    f = [0:m/2]*df;
    E = E(1:(m/2)+1,:);
  else
    [m,n] = size(f);
    if m > n
      f = f';
    end
    %E = [E;E(length(f),:)];
    %f = [f f(length(f))*1.015];
  end
  if sp
    %Add freq spacing between the E points in order to look "seperated"
    [m,n] = size(E);
    Et = zeros(2*m,n);
    ft = zeros(1,2*m);
    Et(1:2:2*m,:) = E;
    ft(1,1:2:2*m) = f;
    df = f*0.015;
    ft(1,2:2:2*m) = f*(1.015);
    E = Et;
    f = ft;
  end
  mesh(t,f,10*log10(abs(E)+1e-6));view(2); ...
      axis('tight');
  xlabel('time (seconds)');
  ylabel('frequency (Hz)');

%Slide a window of width Nt and detect ML tones
%ML detection
function E = ml_detect(E,Nt)
  [r,c] = find(E > 0);
  %Remove the least likely tones using a window of width Nt
  %to determine ML tone
  sub = floor((Nt-1)/2);
  add = ceil((Nt-1)/2); 
  i=1;
  [m,n]=size(E);
  while i < length(c)
    ind = find(c(i)==c);
    %If Multiple likely tones found
    if length(ind) > 1
      ns = max(c(i)-sub,1);
      ne = min(c(i)+add,n);
      Ewin = mean(E(:,ns:ne),2);
      [emax,tone]=max(Ewin);
      Ekeep = E(tone(1),c(i));
      E(:,c(i))=0*E(:,c(i));
      %Keep only the most likely
      E(tone(1),c(i))=Ekeep;
    end
    i = i + length(ind);
  end
