function tones = dtmfdial( string, fs, tonedur, pausedur)
% dtmfdial - create DTMF (TouchTone) signal
%
% usage:
%   tones = dtmfdial( string)
%   tones = dtmfdial( string, fs)
%   tones = dtmfdial( string, fs, tonedur, pausedur)
%
% inputs:
%   string -  numbers to dial (e.g. '1234567890*#')
%   fs      - sampling rate (optional, defaults to 8192Hz)
%   tonedur - duration of tone (optional, defaults to 0.5s)
%   pausedur- duration of pause between tones (optional, defaults to 0.1s)
%

%% check inputs
if ( nargin < 4 )
    pausedur = 0.1;
    
    if ( nargin < 3 )
        tonedur = 0.5;
        
        if ( nargin < 2 )
            fs = 8192;
        end
    end
end

if ( nargin < 1 )
    error( 'dtmfdial must have at least one input' );
end

%% lookup table to translate digits string into numbers
% the following produces the ASCII code for all buttons on the phone pad
Digits = double('123456789*0#');

% then assign each of these ASCII codes a number between 1 and 12; this
% number indicates a button's position on the dial pad
InverseDigits = zeros(1,max(Digits) );
for kk=1:12
    InverseDigits( Digits(kk) ) = kk;
end

%% Parse the diall string
% turn dial string into ASCII code ...
RawNumbers = double( string );
% ... then map ASCII codes into button numbers between 1 and 12
numbers = InverseDigits( RawNumbers );

% ensure numbers are integers between 1 and 12
numbers = round( numbers ); % silently discard fractional part
if ( min( numbers ) < 1 || max( numbers ) > 12 )
    error( 'invalid dial string!' );
end

%% construct signal
% convert durations to number of samples
Ntone  = round( fs*tonedur );
Npause = round( fs*pausedur);

% figure out how long the output signal will be
Nnumbers = length( numbers );
Nsamples = Nnumbers*(Ntone + Npause); 

tones = zeros(1, Nsamples );
pause = zeros(1, Npause);

% associate numbers with DTMF pairs (12-by-2); record normalized frequencies!
dtmfpairs = ...
    [ 697 697 697 770 770 770 852 852 852 941 941 941;
      1209 1336 1477 1209 1336 1477 1209 1336 1477 1209 1336 1477 ]/fs;

% loop over all numbers
for kk = 1:length(numbers)
    Start = (kk-1)*(Ntone + Npause) + 1;
    End   = kk*(Ntone + Npause);
    
    freqs = dtmfpairs( :, numbers(kk) );
    currtone = 0.5* cos( 2*pi*freqs(1)*(0:Ntone-1) ) + ...
        0.5*cos( 2*pi*freqs(2)*(0:Ntone-1) );
    tones(Start:End) = [ currtone pause ];
end