%《數字信號處理教程——MATLAB釋義與實現》第九章例9.4.1程序hc941
% DTMF雙頻撥號信號的生成和檢測程序
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
%
d=input('鍵入一位電話號碼= ','s');  % 輸入一個字符
symbol=abs(d);                      % 求它的ASCII碼
tm=[49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68];   % 16個ASCII碼
for p=1:4;
    for q=1:4;
        if tm(p,q)==abs(d); break,end       % 檢測碼相符的列號q
    end
        if tm(p,q)==abs(d); break,end      % 檢測碼相符的行號p
end
f1=[697,770,852,941];                   % 行頻率向量
f2=[1209,1336,1477,1633];               % 列頻率向量
n=0:2040;                               % 為了發聲，加長序列
x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 構成雙頻信號
sound(x);                               % 發出聲音
disp('雙頻信號已經生成并發出'),pause
% 接收檢測端的程序
%k = [18 20 22 24 31 34 38 42];          % 要求的DFT樣本序號
N=205;k=round([f1,f2]/8000*N);
for m = 1:8;
    X(m)=gfft(x(1:N),k(m));              % 用Goertzel算法計算八點DFT樣本
end
val = abs(X);                           % 列出八點DFT向量
stem(k,val,'.');grid;xlabel('k');ylabel('|X(k)|') % 畫出DFT(k)幅度
set(gcf,'color','w')                    % 置圖形背景色為白
shg,disp('圖上顯示的是檢測到的八個近似基頻的DFT幅度');pause % 
limit = 80;                 %
for s=5:8;
    if val(s) > limit, break, end       % 查找列號
end
for r=1:4;
    if val(r) > limit, break, end       % 查找行號
end
disp(['接收端檢測到的號碼為',setstr(tm(r,s-4))])  % 顯示接收到的字符