%% Dual-Tone Multi-Frequency (DTMF) Signal Detection
% This demonstration showcases the use of the Goertzel function as a part
% of a DFT-based DTMF detection algorithm.
% 
% Dual-tone Multi-Frequency (DTMF) signaling is the basis for voice
% communications control and is widely used worldwide in modern telephony
% to dial numbers and configure switchboards. It is also used in systems
% such as in voice mail, electronic mail and telephone banking.
%
% <<phonepad.png>>

%   Copyright 1988-2005 The MathWorks, Inc.
%   $Revision: 1.1.6.5 $  $Date: 2005/11/18 14:22:43 $

%% Generating DTMF Tones
% A DTMF signal consists of the sum of two sinusoids - or tones - with
% frequencies taken from two mutually exclusive groups. These frequencies
% were chosen to prevent any harmonics from being incorrectly detected
% by the receiver as some other DTMF frequency. Each pair of tones
% contains one frequency of the low group (697 Hz, 770 Hz, 852 Hz, 941 Hz)
% and one frequency of the high group (1209 Hz, 1336 Hz, 1477Hz) and
% represents a unique symbol. The frequencies allocated to the push-buttons
% of the telephone pad are shown below:
%
%                          1209 Hz   1336 Hz   1477 Hz   
%                         _ _ _ _ _ _ _ _ _ _ _ _ _ _ _    
%                        |         |         |         |  
%                        |         |   ABC   |   DEF   |  
%                697 Hz  |    1    |    2    |    3    |     
%                        |_ _ _ _ _|_ _ _ _ _|_ _ _ _ _|  
%                        |         |         |         | 
%                        |   GHI   |   JKL   |   MNO   |  
%                770 Hz  |    4    |    5    |    6    |       
%                        |_ _ _ _ _|_ _ _ _ _|_ _ _ _ _|   
%                        |         |         |         |  
%                        |   PRS   |   TUV   |   WXY   |  
%                852 Hz  |    7    |    8    |    9    |         
%                        |_ _ _ _ _|_ _ _ _ _|_ _ _ _ _|   
%                        |         |         |         |  
%                        |         |         |         |  
%                941 Hz  |    *    |    0    |    #    |       
%                        |_ _ _ _ _|_ _ _ _ _|_ _ _ _ _|  
%                                                           
% First, let's generate the twelve frequency pairs
clear all;

errores=[];
eje_x=10:-10:-80; %rango de valores de SNR
% SNR=-20;
for SNR=eje_x,
error_SNR=[];
    for contador=1:10, %numero de iteraciones para cada SNR

symbol = {'1','2','3','4','5','6','7','8','9','*','0','#'};
lfg = [697 770 852 941]; % Low frequency group
hfg = [1209 1336 1477];  % High frequency group
f0  = [];
f   = [];
for c=1:4,
    for r=1:3,
        f0 = [ f0 [lfg(c);hfg(r)] ];
    end;
end;
anhadir=[[lfg(1);hfg(1)] [lfg(1);hfg(2)] [lfg(1);hfg(3)] [lfg(2);hfg(1)]];
f = [f0 f0 f0 f0 f0 f0 f0 f0 anhadir];
[fil col]=size(f);
%%
% Next, let's generate and visualize the DTMF tones
Fs  = 8000;       % Sampling frequency 8 kHz
N = 800;          % Tones of 100 ms
t   = (0:N-1)/Fs; % 800 samples at Fs
pit = 2*pi*t;

tones = zeros(N,size(f,2));
for toneChoice=1:col,
    % Generate tone
    tones(:,toneChoice) = sum(sin(f(:,toneChoice)*pit))';
    % Plot tone
%     subplot(4,3,toneChoice),plot(t*1e3,tones(:,toneChoice));
%     title(['Symbol "', symbol{toneChoice},'": [',num2str(f(1,toneChoice)),',',num2str(f(2,toneChoice)),']'])
%     set(gca, 'Xlim', [0 25]);
%     ylabel('Amplitude');
%     if toneChoice>9, xlabel('Time (ms)'); end
end
%  set(gcf, 'Color', [1 1 1], 'Position', [1 1 1280 1024])
% annotation(gcf,'textbox', 'Position',[0.38 0.96 0.45 0.026],...
%     'EdgeColor',[1 1 1],...
%     'String', '\bf Time response of each tone of the telephone pad', ...
%     'FitHeightToText','on');

%% Playing DTMF Tones
% Let's play the tones of phone number 508 647 7000 for example. Notice
% that the "0" symbol corresponds to the 11th tone.
 %for i=[5 11 8 6 4 7 7 11 11 11],
 %    p = audioplayer(tones(:,i),Fs);
 %    play(p)
 %    pause(0.5)
 %end
 
%% Estimating DTMF Tones with the Goertzel Algorithm
% The minimum duration of a DTMF signal defined by the ITU standard is 40
% ms. Therefore, there are at most 0.04 x 8000 = 320 samples available for
% estimation and detection. The DTMF decoder needs to estimate the
% frequencies contained in these short signals. 
%
% One common approach to this estimation problem is to compute the
% Discrete-Time Fourier Transform (DFT) samples close to the seven
% fundamental tones. For a DFT-based solution, it has been shown that using
% 205 samples in the frequency domain minimizes the error between the
% original frequencies and the points at which the DFT is estimated.
Nt = 205;
original_f = [lfg(:);hfg(:)];  % Original frequencies
%%
k = round(original_f/Fs*Nt);  % Indices of the DFT
estim_f = round(k*Fs/Nt);      % Frequencies at which the DFT is estimated
%%
% To minimize the error between the original frequencies and the points at
% which the DFT is estimated, we truncate the tones, keeping only 205
% samples or 25.6 ms for further processing.
tones = tones(1:205,:);

%anhadir ruido a la sehal
%%%%%%%%%% prueba con rand
% x=rand(size(tones));
% tones =sqrt(x)+tones;
%%%%%%%%%%%%%%%%%%%%

 x1=auread ('male_voice.au');
 x2=auread ('female_voice.au');
for i=1:col,
tones(:,i)=tones(:,i)+x2(1:205);
end;

tones= awgn (tones,SNR); % ruido anhadido 
%%
% At this point we could use the Fast Fourier Transform (FFT) algorithm to
% calculate the DFT. However, the popularity of the Goertzel algorithm in this
% context lies in the small number of points at which the DFT is estimated.
% In this case, the Goertzel algorithm is more efficient than the FFT algorithm.
%
% Plot Goertzel's DFT magnitude estimate of each tone on a grid
% corresponding to the telephone pad.
% figure,
numeros=[];
num=0;
% validity=zeros(1,12);
for toneChoice=1:col,
    % Select tone
    tone=tones(:,toneChoice);
    % Estimate DFT using Goertzel
    ydft(:,toneChoice) = goertzel(tone,k+1); % Goertzel use 1-based indexing
    % Plot magnitude of the DFT
%     subplot(4,3,toneChoice),stem(estim_f,abs(ydft(:,toneChoice)));
%     title(['Symbol "', symbol{toneChoice},'": [',num2str(f(1,toneChoice)),',',num2str(f(2,toneChoice)),']'])
%     set(gca, 'XTick', estim_f, 'XTickLabel', estim_f, 'Xlim', [650 1550]);
%     ylabel('DFT Magnitude');
%     if toneChoice>9, xlabel('Frequency (Hz)'); end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%validity checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% One_tone=abs(ydft(:,toneChoice));
% [valor1,max1]=max(One_tone);
% One_tone(max1)=0;
% M1=estim_f(max1);
% [valor2,max2]=max(One_tone);
% One_tone(max2)=0;
% M2=estim_f (max2);
% [valor3,max3]=max(One_tone);
% M3=estim_f (max3);
% cond1=0;cond2=0;cond3=0;T1=30;T2=.6;
% if ((M1 > 1000) & (M2 <1000))|((M2 > 1000) & (M1 <1000)),cond1=1;end;
% if (abs(M2))> T1, cond2=1; end;
% if (M2/M3)>T2, cond3=1; end;
% if (cond1==1)&(cond2==1)&(cond3==1),validity(toneChoice)=1; end;
%end;
% validity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Detecting DTMF Tones
% The digital tone detection can be achieved by measuring the energy 
% present at the seven frequencies estimated above. Each symbol can be
% separated by simply taking the component of maximum energy in the lower
% and upper frequency groups.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
One_tone=abs(ydft(:,toneChoice));
[valor1,max1]=max(One_tone);
One_tone(max1)=0;
M1=estim_f(max1);
[valor2,max2]=max(One_tone);
One_tone(max2)=0;
M2=estim_f (max2);
if ((M1==702)&(M2==1210))|((M2==702)&(M1==1210)), num=1;
elseif ((M1==702)&(M2==1327))|((M2==702)&(M1==1327)), num=2;
elseif ((M1==702)&(M2==1483))|((M2==702)&(M1==1483)), num=3;
elseif ((M1==780)&(M2==1210))|((M2==780)&(M1==1210)), num=4;
elseif ((M1==780)&(M2==1327))|((M2==780)&(M1==1327)), num=5;
elseif ((M1==780)&(M2==1483))|((M2==780)&(M1==1483)), num=6;
elseif ((M1==859)&(M2==1210))|((M2==859)&(M1==1210)), num=7;
elseif ((M1==859)&(M2==1327))|((M2==859)&(M1==1327)), num=8;
elseif ((M1==859)&(M2==1483))|((M2==859)&(M1==1483)), num=9;
elseif ((M1==937)&(M2==1210))|((M2==937)&(M1==1210)), num=10;
elseif ((M1==937)&(M2==1327))|((M2==937)&(M1==1327)), num=11;
elseif ((M1==937)&(M2==1483))|((M2==937)&(M1==1483)), num=12;
else num=20; 
end;

numeros=[numeros num];
end;

numeros;
a= [1 2 3 4 5 6 7 8 9 10 11 12];
correctos=[a a a a a a a a 1 2 3 4];
error=length(find(correctos-numeros)); %errores en una secuencia de 100 tonos
error_SNR= [error_SNR error]; %errores en 10 iteraciones para una SNR dada
end;
% error_grafica
promedio_error_SNR= sum(error_SNR(1:end))/(length(error_SNR));
errores=[errores promedio_error_SNR];
end;

figure;
% plot (eje_x,error_grafica); 
plot(eje_x,errores); grid;
title('Probability of false detection');
ylabel('% error');
xlabel ('SNR (dB)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set(gcf, 'Color', [1 1 1], 'Position', [1 1 1280 1024])
% annotation(gcf,'textbox', 'Position',[0.28 0.96 0.45 0.026],...
%     'EdgeColor',[1 1 1],...
%     'String', '\bf Estimation of the frequencies contained in each tone of the telephone pad using Goertzel', ...
%     'FitHeightToText','on');