数字信号处理实验2

一、实验目的

1、产生数字8的双音频信号的离散信号，其是对双音频信号时长数为500ms的连续信号采样获得。

2、用N=205点的戈泽尔算法，识别出“8”这个符号。

（注意，要找好需查找的 k1 和 k2 )

码音时长：250~500ms (注：静音时，信号为0)

静音	号码1	静音	号码2				静音	号码8	静音

静音时长：200ms~400ms

单号码音时长：250~500ms (注：静音时，信号为0)

用N=205点的戈泽尔算法，识别出电话号码。

4、在理想的拨号音上叠加一个50Hz的工频干扰信号。

用N=205点的戈泽尔算法，识别出电话号码。

5、在理想的拨号音上叠加一个白噪声信号。

在拨号信号持续期内的信噪比分别为10dB，5dB，3dB，2dB，1dB，0dB的情况下，用N=205点的戈泽尔算法，识别出电话号码。

二、实验原理

DTMF是按键式电话信令中的一般名称，它等效于在贝尔系统内部正在使用的按纽式拨号系统。DTMF在电子邮件系统和电话银行系统中也得到广泛的应用；后者，用户可以从电话送出DTMF信号根据一个菜单来任选一项。

DTMF信号很容易用软件产生并用数字滤波器检测到，也能用软件来实现，做成是对8个频率单音调谐的。通常，DTMF信号经由一个编码/解码芯片或用线性A/D和D/A转换器与模拟环境接口。编解码芯片包含了一个双向数字/模拟接口全部必要的A/D和D/A，采样和滤波电路。

DTMF音可以用数学方法产生，也可以用查表方式产生。在一种硬件实现（例如在数字信号处理器）中，两个正弦波的数字样本用数学方法产生、加权并相加在一起。其和经对数压缩后送到编解器转换到一模拟信号。在8kHz的采样率下，硬件必须每隔0.125ms输出一个样本。在这种情况下，不用正弦函数查表的办法，这是由于正弦波的值能很快计算出而勿需像查表那样要求占用大的数据存储量。在MATLAB中为了仿真和研究的目的，查表法可能是一个好的办法。

在接收端接收到这个来自编解码器的经对数压缩的8bit数字数据字后，将它们对数扩展到16bit的线性格式，然后检测到这些单音并判断出传递的数字。检测算法可以是利用FFT算法的一种DFT实现，或者是一种滤波器柜的实现。

三、实验条件

1．使用MATLAB及其函数进行实验

四、实验步骤

1．分析问题

2．MATLAB运行

3．分析实验结果

五、实验代码

主代码：

clc
clear all
f0 = 8000;
T = 1/f0;

% 1、产生数字8的双音频信号的离散信号，其是对双音频信号时长数为500ms的连续信号采样获得。
% 每隔0.125ms，时长数为500ms的连续信号，共4000个
t = 0:1:3999;
[h, l] = mfind('8');
q1 = sin(2*pi*h*t*T)+sin(2*pi*l*t*T);
% sound (q1,f0)

% 2、用N=205点的戈泽尔算法，识别出“8”这个符号。（注意，要找好需查找的k1和k2）
s2 = shibie(q1)

% 3、产生一理想的8位数的电话拨号音信号（例如82317231）,用N=205点的戈泽尔算法，识别出电话号码。
% dail = '#*D0CA12';     %设置8位数的电话拨号音信号
% dail = '56478521';     %设置8位数的电话拨号音信号
% dail = '12367889';     %设置8位数的电话拨号音信号
dail = '*0#D7BA1';     %设置8位数的电话拨号音信号
q3 = [];
for i = 1:1:8
    %200ms~400ms对应1600~3200，250~500ms对应2000~4000
    [h, l] = mfind(dail(i));
    t1 = round(1600 + 1600*rand(1));
    t20 = round(2000 + 2000*rand(1));
    t2 = 0:1:t20;
    q3 = [q3,zeros(1,t1),sin(2*pi*h*t*T)+sin(2*pi*l*t*T)];
end
t1 = round(1600 + 1600*rand(1));
q3 = [q3,zeros(1,t1)];
% sound (q3,f0)
i = 1;
s3 = liushibie(q3)

% 4、在理想的拨号音上叠加一个50Hz的工频干扰信号。用N=205点的戈泽尔算法，识别出电话号码。
l = length(q3);
t = 1:1:l; 
q4 = q3 + cos(2*pi*50*t*T);
s4 = liushibie(q4)

% 5、在理想的拨号音上叠加一个白噪声信号。
% 在拨号信号持续期内的信噪比分别为10dB，5dB，3dB，2dB，1dB，0dB的情况下，用N=205点的戈泽尔算法，识别出电话号码。
% （分析哪种情况不能识别?）
q510 = awgn(q3, 10);
s510 = liushibie(q510)
q55 = awgn(q3, 5);
s55 = liushibie(q55)
q53 = awgn(q3, 3);
s53 = liushibie(q53)
q52 = awgn(q3, 2);
a52 = filter(lbq, q52);
s52 = liushibie(a52)
q51 = awgn(q3, 1);
a51 = filter(lbq, q51);
s51 = liushibie(a51)
q50 = awgn(q3, 0);
a50 = filter(lbq, q50);
s50 = liushibie(a50)

figure(1)
subplot(3,1,1)
plot(q1)
title('数字8的双音频信号的离散信号')
hold on
subplot(3,1,2)
plot(q3)
title('理想电话拨号音信号')
hold on
subplot(3,1,3)
plot(q4)
title('理想的拨号音上叠加一个50Hz的工频干扰信号')
hold on

figure(2)
subplot(6,1,1)
plot(q510)
title('叠加信噪比分别为10dB白噪声信号的拨号音')
hold on
subplot(6,1,2)
plot(q55)
title('叠加信噪比分别为5dB白噪声信号的拨号音')
hold on
subplot(6,1,3)
plot(q53)
title('叠加信噪比分别为3dB白噪声信号的拨号音')
hold on
subplot(6,1,4)
plot(q52)
title('叠加信噪比分别为2dB白噪声信号的拨号音')
hold on
subplot(6,1,5)
plot(q51)
title('叠加信噪比分别为1dB白噪声信号的拨号音')
hold on
subplot(6,1,6)
plot(q50)
title('叠加信噪比分别为0dB白噪声信号的拨号音')
hold on

mfind函数：

function [h, l] = mfind(x)
    if (x == '1')
        h = 697;l = 1209;
    elseif x == '2'
        h = 697;l = 1336;
    elseif x == '3'
        h = 697;l = 1477;
    elseif x == '4'
        h = 770;l = 1209;
    elseif x == '5'
        h = 770;l = 1336;
    elseif x == '6'
        h = 770;l = 1477;
    elseif x == '7'
        h = 852;l = 1209;
    elseif (x == '8')
        h = 852;l = 1336;
    elseif x == '9'
        h = 852;l = 1477;
    elseif x == '0'
        h = 941;l = 1336;
    elseif x == '#'
        h = 941;l = 1209;
    elseif x == '*'
        h = 941;l = 1477;
    elseif x == 'A'
        h = 697;l = 1633;
    elseif x == 'B'
        h = 770;l = 1633;
    elseif x == 'C'
        h = 852;l = 1633;
    elseif x == 'D'
        h = 941;l = 1633;
    end
end

shibie函数

function [s] = shibie(y)
    V(1)=DFT(18,y);
    V(2)=DFT(20,y);
    V(3)=DFT(22,y);
    V(4)=DFT(24,y);
    V(5)=DFT(31,y);
    V(6)=DFT(34,y);
    V(7)=DFT(38,y);
    V(8)=DFT(42,y);
    if V(1)>=80&&V(5)>=80
        s='1';
    elseif V(1)>=80&&V(6)>=80
        s='2';
    elseif V(1)>=80&&V(7)>=80
        s='3';
    elseif V(1)>=80&&V(8)>=80
        s='A';
    elseif V(2)>=80&&V(5)>=80
        s='4';
    elseif V(2)>=80&&V(6)>=80
        s='5';
    elseif V(2)>=80&&V(7)>=80
        s='6';
    elseif V(2)>=80&&V(8)>=80
        s='B';
    elseif V(3)>=80&&V(5)>=80
        s='7';
    elseif V(3)>=80&&V(6)>=80
        s='8';
    elseif V(3)>=80&&V(7)>=80
        s='9';
    elseif V(3)>=80&&V(8)>=80
        s='C';
    elseif V(4)>=80&&V(5)>=80
        s='#';
    elseif V(4)>=80&&V(6)>=80
        s='0';
    elseif V(4)>=80&&V(7)>=80
        s='*';
    elseif V(4)>=80&&V(8)>=80
        s='D';
    else
        s='E';
    end
end

liushibie函数

function [S] = liushibie(x)
%     clear
    m=length(x);
%     global flag
    l=1;
    b=0;
    sum=0;
    for i=1:1:100
        sum=sum+x(i)^2;
    end
    s(1)=sum/100;
    for i=2:1:(m-99)
        s(i)=(100*s(i-1)-x(i-1)^2+x(i+99)^2)/100;
    end
     if s(1)<0.65
         h=0;
     else
         h=1;
     end
    for i=1:1:(m-99)
        if s(i)<0.65
            h=0;
        end
        if (h==0)&&(s(i)>0.8)&&(s(i)>=s(i-1)&&s(i)>=s(i+1))
            S(l)=Sampling(x(i+100:1:i+304)); %富裕100个点取205个点到304，即100-304共205个点
            if S(l)=='E'||l>8
                S='已不能识别';
%                 flag=1;
                break
            end
            l=l+1;
            h=1;
        end
    end
%     if strcmp(s(1),'E')
%         S='已不能识别';

DFT函数

function [v] = DFT(k,y)
v=0;
    for i=0:1:204
        v=v+y(i+1)*exp(-j*2*pi/205*k*i);
    end
v=abs(v);
end

Sampling函数

function [s] = Sampling(y)
    V(1)=DFT(18,y);
    V(2)=DFT(20,y);
    V(3)=DFT(22,y);
    V(4)=DFT(24,y);
    V(5)=DFT(31,y);
    V(6)=DFT(34,y);
    V(7)=DFT(38,y);
    V(8)=DFT(42,y);
    if V(1)>=80&&V(5)>=80
        s='1';
    elseif V(1)>=80&&V(6)>=80
        s='2';
    elseif V(1)>=80&&V(7)>=80
        s='3';
    elseif V(1)>=80&&V(8)>=80
        s='A';
    elseif V(2)>=80&&V(5)>=80
        s='4';
    elseif V(2)>=80&&V(6)>=80
        s='5';
    elseif V(2)>=80&&V(7)>=80
        s='6';
    elseif V(2)>=80&&V(8)>=80
        s='B';
    elseif V(3)>=80&&V(5)>=80
        s='7';
    elseif V(3)>=80&&V(6)>=80
        s='8';
    elseif V(3)>=80&&V(7)>=80
        s='9';
    elseif V(3)>=80&&V(8)>=80
        s='C';
    elseif V(4)>=80&&V(5)>=80
        s='#';
    elseif V(4)>=80&&V(6)>=80
        s='0';
    elseif V(4)>=80&&V(7)>=80
        s='*';
    elseif V(4)>=80&&V(8)>=80
        s='D';
    else
        s='E';
    end
end

lbq函数

function Hd = lbq
%LBQ 返回离散时间滤波器对象。

% MATLAB Code
% Generated by MATLAB(R) 9.9 and DSP System Toolbox 9.11.
% Generated on: 18-Dec-2021 19:59:49

% Equiripple Lowpass filter designed using the FIRPM function.

% All frequency values are in Hz.
Fs = 8000;  % Sampling Frequency

Fpass = 1650;            % Passband Frequency
Fstop = 1700;            % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

% [EOF]

六、实验结果

s2 =
‘8’
s3 =
‘*0#D7BA1’
s4 =
‘*0#D7BA1’
s510 =
‘*0#D7BA1’
s55 =
‘*0#D7BA1’
s53 =
‘*0#D7BA1’
s52 =
‘*0#D7BA1’
s51 =
‘*0#D7BA1’
s50 =
‘*0#D7BA1’

七、结果分析

1、可以通过这种方式稳定的传输信噪比大于3dB以上的信号。

2、运算过程很快，适合实际中使用。