function dtmfdec(wavfn)
%sivann 2008
%input: wavfn.wav file
%output: time + detected letter combination in wavfn.txt
%calls dtmfparse to convert dtmf digits to NMEA
%
% for video input: c:\winap\ffmpeg.exe -i c:\temp\ilisia2fotis.divx.pcm.avi  -ab 16k out.wav
%C:\winap\ffmpeg>ffmpeg.exe -i c:\temp\ilisia2fotis.divx.pcm.avi  -ar 16000 -ac 1 c:\temp\out.wav
%
% TODO: wavread in parts to avoid 'out of memory' for large wavs
%[samples channels]=wavread('c:/temp/out44.wav','size');
%[Y,FS,NBITS]=wavread('c:/temp/out44.wav',[1 10000000]);
%



[Y,FS,NBITS]=wavread(wavfn); % USAGE: dtmfdec('filename.wav');
[pat,wfn,ext,ver]=fileparts(wavfn);
txtfn=fullfile(pat,[wfn '.txt']);
fid=fopen(txtfn,'w'); disp(['truncating ' txtfn]);

f = [697 770 852 941   1209 1336 1477 1633]; %dtmf freqs
letters=[ '1' '2' '3' 'A' ;'4' '5' '6' 'B' ;'7' '8' '9' 'C';'*' '0' '#' 'D' ];

Y=sum(Y',1); %add both audio channels

N = 205;%# of samples per goertzel (fft) (bin size) http://ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html
Fs=16000;%convert to this freq first
Y=Y(1:(FS/Fs):end); %subsamle (sucks)

y=Y(1,:); %stereo->mono

% Indices of the DFT 
k = round(f/Fs*N); 
% Frequencies at which DFT is actually estimated  
ef = round(k*Fs/N); 

%create a band-pass filter
[b,a]=cheby1(5,0.5,[500/Fs*2 1800/Fs*2]');
%[H,w] = freqz(b,a,512); plot(w*Fs/(2*pi),abs(H));grid;  %plot filter

%filter signal
yf=filter(b,a,y);
%yf=y;
silence=0;

%find max value ymax
ymax=[];
for i=1:N:(length(yf)-N)
  y1=yf(i:i+N-1);
  ydft = goertzel(y1,k+1); 
  ymax=max([max(abs(ydft)) ymax]);
end  

%init variables
detected=[];fletters=[];seq=[];seqidx=1;isvalid=0;

%loop through blocks of N wave samples
for i=160000:N:length(yf)-N
  y1=yf(i:i+N-1);
  tstart=i.*1/Fs ; tend=(i+N-1).*1/Fs ;%seconds
  ydft = goertzel(y1,k+1); 
  %stem(ef,abs(ydft));  AXIS([f(1) f(end) 0 50 ]);pause(0.02) %plot values
  xt=abs(ydft)>(ymax/3); % values above thresh (binary)
  
  if sum(xt) < 2 silence=1; isvalid=0;   letter=[' ' ' '];
  
  elseif sum(xt) == 2 % two maxima 
	isvalid=isvalid+1;
	dtl=find (xt(1:4)); %f 1-4 represent low freqs
	dth=find (xt(5:end));
	%[dtl dth isvalid]
	letter(isvalid)=letters(dtl,dth);
  end

  if (isvalid==2) && (letter(1) ~= letter(2))isvalid=0;  letter=[' ' ' '];end; %wrong letter detected
  
  if (isvalid == 2) %occurs on second detection. following detections will be ignored (isvalid>2)
	detected=[detected; letter(1)  tstart]; %
	fletters=[fletters ;letter(1)]; %found letters
	%sprintf('%c',fletters)
	
    if (seqidx>60) || (letter(1)=='A') seqidx=1;end; %at start
    
	seq(seqidx)=letter(1);seqidx=seqidx+1;


	if (seqidx>40) && (letter(1)=='B') 
	  str=dtmfparse(seq,tstart); % parse the seq string
      fprintf(fid,'%s\n',str);
      seq=[];seqidx=1;

	end; %at end
  end
end

fclose(fid);

