fs = 8192;
tonedur = 0.25;
pausedur = 0.05;

LBPF = 129;

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

% generate signal

xx=dtmfdial('7035569452', fs, tonedur, pausedur);
xx = xx + sigma*randn(size(xx));
soundsc(xx,fs)
subplot(1,1,1)
specgram(xx, 256, fs)

pause

% 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 );
end

% plot frequency response
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 by each of the seven 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

% decision logic
% decompose samples into periods for each number
Nnumbers = floor( length(xx)/(fs*(tonedur+pausedur)) );
NTonePlusPause = round(fs*(tonedur+pausedur));
NPause = round(fs*pausedur);
    
% score for each tone period: sum of squares in period
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')
set(gca, 'YTickLabel',...
    [ ' 691'; ' 770'; ' 852'; ' 941'; '1209'; '1336'; '1477'])

pause

% Decisions
% in each row of the score matrix find the biggest entry among the first
% four and final three columns
for nn=1:Nnumbers
    [ smax, imax_low(nn)] = max( score(nn, 1:4) );
    [ smax, 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; % between 1 and 12
rawnumbers = Digits( numbers );       % ASCII codes for dialed numbers
decodedNumber = char( rawnumbers )    % convert to string