%dtmfMain - script for demonstrating DTMF dial tones
%
% This script demonstrates dual-tone multiple-frequency (DTMF) dial tones.
% It illustrates how these tones can be generated in MATLAB and shows a
% possible strategy for decoding these signals. 
%
% This script is simplified by assuming that tone durations are known and
% all identical. 

%% set parameters
fs = 8192;
tonedur = 0.25;
pausedur = 0.05;

LBPF = 150; % length of bandpass filters

SNR = 60; % signal-to noise ratio in dB
sigma = sqrt(0.5/10^(SNR/10));

numberToDial = '7039934810';

%% generate signal
% dtmfdial parses the string containing the phone number and generates the
% signals
xx = dtmfdial(numberToDial, fs, tonedur, pausedur);
xx = xx + sigma*randn(size(xx)); % add noise
soundsc(xx,fs)
subplot(1,1,1)
spectrogram(xx, 256, 128, 256, fs,'yaxis');

pause

%% Recovering the dial string
% the strategy for decoding the signal and to recover the dial string is as
% follows. A bandpass filter is designed for each of the seven DTMF
% frequencies. Then, the signal is passed through each of the seven
% filters. To determine if a frequency is present during a given dialing
% period the strength (energy) at the output of the filters is measured. The
% pair of frequencies (one from the low group and one from the high group)
% with the highest powers indicates which button was pushed. 

%% design BP filters
DTMFFreqs = [ 691 770 852 941 1209 1336 1477 ];

hh = zeros(length(DTMFFreqs), LBPF);

for kk=1:length(DTMFFreqs)
    %hh(kk,:) = makeBPF( DTMFFreqs(kk)/fs, LBPF );
    hh(kk,:) = makeBPF2( DTMFFreqs(kk)/fs, LBPF-1 );
end

% plot frequency responses
ff = 0 : 0.001 : 0.5;
for kk = 1:length(DTMFFreqs)
    subplot( 7, 1, kk )
    plot( ff*fs, abs( FreqResp( hh(kk,:), ff) ), 'b-', ...
        DTMFFreqs, abs( FreqResp( hh(kk,:), DTMFFreqs/fs) ), 'bo');
    grid
    ylabel( '|H(f)|' )
    %xlabel( 'Frequency (Hz) ')
    axis([0 fs/2 0 1.1])
end

pause

%% filter the signal with each of the seven bandpass filters
yy = zeros(length(DTMFFreqs), length(xx)+length(hh(1,:))-1);

for kk=1:length(DTMFFreqs)
    yy(kk,:) = conv( hh(kk,:), xx );
end

tt = 0 : 1/fs : (length(yy(1,:))-1)/fs;
for kk = 1:length(DTMFFreqs)
    subplot( 7, 1, kk )
    plot( tt, yy(kk,:) );
    grid
    ylabel( 'y(t)' )
    %xlabel( 'Frequency (Hz) ')
    axis([0 (length(yy(1,:))-1)/fs -1 1])
end

pause

% divide samples into periods for each number
Nnumbers = floor( length(xx)/(fs*(tonedur+pausedur)) );
NTonePlusPause = round(fs*(tonedur+pausedur));
NPause = round(fs*pausedur);

%% measure energy at output of bandass filters
% score for each tone period: sum of squares in each period (energy)
score = zeros(Nnumbers, length(DTMFFreqs));
for nn=1:Nnumbers
    Startnn = (nn-1)*NTonePlusPause + 1 + floor(LBPF/2);
    Endnn = nn*NTonePlusPause - NPause - floor(LBPF/2);
    for kk = 1:length(DTMFFreqs)
        score(nn,kk) = sum( yy(kk,Startnn:Endnn).^2 );
    end
end

subplot(1,1,1)
bar3( score' );
xlabel('Digit period')
ylabel('Tone')
zlabel('Score (Energy)')
set(gca, 'YTickLabel',...
    [ ' 691'; ' 770'; ' 852'; ' 941'; '1209'; '1336'; '1477'])

pause

%% Decisions: find frequencies with highest power
% in each row of the score matrix find the biggest entry among the first
% four and final three columns
imax_low = zeros(1, Nnumbers); % column number
imax_high = zeros(1, Nnumbers);% row number on dial pad
for nn=1:Nnumbers
    [ ~, imax_low(nn)] = max( score(nn, 1:4) );
    [ ~, imax_high(nn)] = max( score(nn, 5:7) );
end

%% decode
% lookup table to translate numbers string into numbers
Digits = double('123456789*0#');      % table of ASCII codes for dial-pad

numbers = (imax_low-1)*3 + imax_high; % map row/column to button numbers between 1 and 12
rawnumbers = Digits( numbers );       % Look up ASCII codes for dialed numbers
decodedNumber = char( rawnumbers );   % convert ASCII codes to string

fprintf('The number is: %s\n', decodedNumber);