数字信号处理课程设计

《数字信号处理》课程设计实验报告-PAGE34-中南大学数字信号处理课程设计报告题目数字信号处理课程设计学生姓名阮光武指导教师李宏学院信息科学与工程学院专业班级电子信息工程0704班

目录第一章实验概述 41.1DTF在双音多频信号检测中的应用 41.2频率分布 4第二章总体设计及关键技术 52.1工作过程及关键技术 52.2信号滤波演示 52.2.1工作流程 52.2.2设计要求 6第三章程序实现 73.1DFT在双音多频信号检测中的应用的程序实现 73.1.1代码说明 73.1.2部分代码 73.1.3结果展示 93.2信号滤波演示的程序实现 113.2.1关键部分的概述 113.2.2关键部分的设计 153.2.3结果分析 32第四章结束语 33参考文献 34

第一章 实验概述在日常生活中我们当然会遇到许多的问题,也许外面看上去会很复杂.但是你仔细思考后就会发现原来他们是这么的简单.就如同我们熟悉的电话一样,你还不太明白这电话号码是如何拨号的.但是你看了这个实验之后也就会发现原来是我们所学过的.DFT在双音多频信号（DTMF）检测中的应用及信号滤波演示。你肯定也见过那些触摸屏拨号的,很奇怪他们是如何工作的.下面的实验我们就会了解这些问题.当然在实验中我们肯定会遇到许多的难题,但是我们要努力解决他们,成功就会在眼前.1.1DFT在双音多频信号检测中的应用双音多频信号是电话系统中电话机与交换机之间的一种用户信令，通常用于发送被叫号码。在使用双音多频信号之前，电话系统中使用一连串的断续脉冲来传送被叫号码，称为脉冲拨号。电话机向交换机同时发送两个频率的信号，发送按键信息。一共有8个频率的音频信号，分为2组，每组4个，两两组合共可以代表16个按键，分别代表0-9、#、*等按键。1.2频率分布频率组合（标称值）：行频率：697Hz、770Hz、852Hz、941Hz；列频率：1209Hz、1336Hz、1477Hz、1633Hz行频率和列频率组成16种双频数字信号。电话通信系统中拨号频率的分配如下图所示：图1-1电话的频率分配图一个DTMF信号是由一个高频段和一个低频段的音频信号混合而成，可以直接由两个频率的正弦序列相加来实现。通过对两个两个频率的正弦序列相加进行DFT变换，我们可以清晰的看见，在变换后的的幅度中，有两个值远远大于另外的值。那么，根据这个性质，接受端就可以通过技术上的处理来获得原来信息，这也正是DTMF在电话通信系统中的工作原理。DTMF在通信中有一定的指标。CCIT对DTMF信号所规定的指标是：传送/接收率为每秒10个数字，即每个数字100ms。每个数字传送过程中，信号必须存在至少45ms，且不得多于55ms，100ms里的其余时间是静音。实验中取8000HZ的抽样频率，DTMF信号持续50ms,按照CCIT的指标DFT的N=8000*0.05=400。但为了尽量减少泄漏，N=205时最为理想。则有下表：表一8个基音频率及其N=205时对应的DFT序号

第二章 总体设计及关键技术2.1工作过程及关键技术图2-1DTMF信号产生及检测流程从图上可以看出,整个工作流程是从电话号码拨出,然后检测得到的不同电话号码,生成不同的信号.总体来说分为两步:（1）拨号端根据输入的不同电话键值，生成相应的DTMF信号，展示该信号时/频域波形并播放该信号；（2）局端根据接收到的DTMF信号检测出拨号端发送的11位电话号码。根据要求,每播一个号就将其时/频域波形展示出来。当拨号完成后，取出每个号码的目的7个k点(图2.2中的A,B,C,D不考虑在内)并都存在一设定好的文件中。当发送端点击发送后，检测端读取事先约定的文件信息。将文件中的信息读取出来后，连续7个归为一组，通过值的比较便可一一检测到所发送的信息。2.2.信号滤波演示2.2.1工作流程图2.3信号滤波演示系统工作流程2.2.2设计要求信号发生器：能实现信号的导入。有动态导入与静态导入两种。①静态型：直接输入(或从文件读取)测试信号序列；②动态型：输入由多个不同频率正弦信号叠加组合而成的模拟信号公式100sin(2πf1t)+100sin(2πf2t)+…+100sin(2πfnt)、采样频率（Hz）以及采样点数，动态生成该信号的采样序列，作为测试信号。DFT频谱分析—使用DFT对产生的测试信号进行频谱分析并展示其幅频、相频特性，指定需要滤除或保留的频带，通过选择滤波器类型（IIR/FIR），确定对应的滤波器（低通、高通、带通、带阻）技术指标。滤波器设计：根据IIR/FIR数字滤波器技术指标设计滤波器，生成相应的滤波器系数，并展示对应的滤波器幅频（衰减）、相频特性。(1)IIRDF设计：使用双线性变换法，可选择滤波器基型（巴特沃斯或切比雪夫型）；(2)FIRDF设计：使用窗口法，可选择窗口类型。数字滤波：根据设计的滤波器系数，对测试信号进行滤波。(1)IIRDF：要求通过差分方程迭代实现滤波，未知初值置零处理；(2)FIRDF：要求通过快速卷积实现滤波。可以选择使用重叠相加或重叠保留法进行卷积运算，并动态展示卷积运算的详细过程。输出信号分析：展示滤波后信号的幅频与相频特性，分析是否满足滤波要求。对同一滤波要求，根据输出信号频谱，对比分析各类滤波器的差异。

第三章 程序实现3.1.DFT在双音多频信号检测中的应用的程序实现3.1.1代码说明总休来说这个程序还是比较简单的,首先根据拨出的记录查找出拨出的电话号码.然后根据相应的行频率和列频率进行相加,求出相应的dft序列.最后输出号码。程序可以分为两部分。发送与检测。发送端对输入号码的相应频率进行演示。输出端对得到的信息进行分析。发送端：每拨一个号码，便通过函数调用Dtmf函数展示其DFT变换，并将以定义好的7个k值所对应DFT中的值记录在双方的接口缓存中。当码号拨好后，通过发送命令Call激活检测端，使检测端开始检测号码。接受端：当检测端被激活后，便进行译码工作。既将接口中的数据取出后进行重新分组并在每组中找到最大的两个值，通过与表一与图2.1的比较，便可以将7个值译为一个一位号码。3.1.2部分代码(1)转换函数functionDtmf(handles,f1,f2,num)holdon;globalcoun;a=get(handles.edit1,'String');ifsize(a)==[00]coun=0;fid=fopen('number.txt','w');fclose(fid);elsef=coun+102;coun=f;endc=0:102;d=coun+c;T=1/8000;N=205;x=0:T:0.5;b=[a,num];set(handles.edit1,'String',b);y=cos(2*pi*f1*x)+cos(2*pi*f2*x);m=abs(fft(y,N));e=m(1:103);stem(d,e,'.');axis([0,coun+50,0,120]);holdoff;fid=fopen('number.txt','a');k=[19,21,23,25,32,35,39];forn=1:7fprintf(fid,'%f',e(k(n)));endfclose(fid);(2)拨号端译码函数forn=1:31;status=feof(fid);ifstatus==1ifn==1set(ftext,'String','无拨号请求');Break；endset(ftext,'String',str);break;else[nums,m]=fscanf(fid,'%f',7);[num]=compare(nums);a=[str,num];str=a;endifn==31set(ftext,'String','输入号码过长，是无效号码');endend上面为检测端的一for循环，作用是对端口输入的信息进行译码。规则为7个一组，并通过compare函数对7个数进行分析，进一步得到7个数所表示的一为号码。for循环中，定义了拨号序列的最大长度为30，否则就将此序列是为无效号码，给予处理，并报告无效信息。其中compare函数如下:(3)检测端译码函数function[num]=compare(nums);[y,k]=max(nums);nums(k)=0;[y,n]=max(nums);ifk<nm=k;k=n;n=m;endm=k;k=m-4;arry=['1''2''3';'4''5''6';'7''8''9';'*''0''#'];ifk<4&k>0ifn<5&n>0num=arry(n,k);endelsenum='(错传号码)';end此函数为检测端的译码函数，其对事现得到的7个数进行分析，其中分析为：定义矩阵A=，在number.Txt中，从开始将数据每7个分为一组，找到较大两个值的序号，以图3.5的第五行为例，7个值分别为：14.24450390.51359510.4162075.2545730.9508101.913778100.420003，可以很容易找到较大两值对应的序号为2和7，并将两个值中较大的减4作为纵坐标，较小的作为横坐标，既有A(2,7-4)=A(2,3)=6。通过比较，与所拨号码一致。3.1.3、结果展示图3-1发送端与检测端的演示图3-2号码2的DFT从上图可以看出，最大两点对应的k为18与34。在表一中可以查到k=18对应的频率为697HZ，k=34对应的频率为1336HZ。在与图2.1拨号电话的频率分配图比较可得，此图为号码2的DFT变换。图3-3号码5的DFT同理，k为20与34，通过与表一与图2.1的比较可得此图为号码5的DFT变换。图3-41511633的DFT序列上图展示了1511633序列的每一个号码的DFT变换。图3-5图3-4所播号码对应的两端接口数据Number.txt为默认的两端通信接口。从上图可以看出，第一行有7个值，而且有第一个，第五个值远远大于其他值，这时，对应的k为18，31，通过上面的查找方法即可得到其对应号码为1。第二行同样有7个值，同理有较大两值对应的k分别为20,34，同理可得到其对应号码为5。另外，通过对表一与图2.1的对比。我们得到一种简单的找到号码的方法如compare函数的描述。3.2.信号滤波演示的程序实现3.2.1.关键部分的概述(1)用窗函数法设计FIR滤波器根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为Δw，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。用窗函数wd(n)将hd(n)截断，并进行加权处理，得到如果要求线性相位特性，则h(n)还必须满足：根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。(2)用双线性变换法设计IIR数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Ｚ平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图3-6所示。图3-6双线性变换的映射关系为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上，可以通过以下的正切变换实现（1）式中,T仍是采样间隔。当Ω1由-π/T经过0变化到π/T时，Ω由-∞经过0变化到+∞，也即映射了整个jΩ轴。将式（1）写成将此关系解析延拓到整个S平面和S1平面，令jΩ=s，jΩ1=s1，则得再将S1平面通过以下标准变换关系映射到Z平面z=es1T从而得到S平面和Z平面的单值映射关系为：(2)(3)式（2）与式（3）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式（1）与式（2）的双线性变换符合映射变换应满足的两点要求。首先,把z=ejω，可得(4)即S平面的虚轴映射到Z平面的单位圆。其次，将s=σ+jΩ代入式（4），得因此由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个jΩ轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式（4）所示，重写如下：上式表明，S平面上Ω与Z平面的ω成非线性的正切关系，如图3-7所示。由图3-7看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。图3-7双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（4）及图3-7所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图3-8所示。图3-8双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。3.2.2.关键部分的设计(1)语音信号的采集：利用windows下的录音机（开始—程序—附件—娱乐—录音机，文件—属性—立即转换—8000KHz，8位，单声道），录制一段自己的话音“信号”，时间控制在1秒左右，然后将音频文件保存“xh.wav”在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。(2)语音信号的频谱分析①首先画出语音信号的时域波形z1=wavread('d:\laoA.wav');plot(z1);图像输出如图3-9图3-9语音信号时域波形②对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性z1=wavread('d:\laoA.wav');y1=z1(1:8192);Y1=fft(y1);n=0:8191;plot(n,Y1);图像输出如图3-10：图3-10语音信号频谱分析图(3)设计数字滤波器和对信号滤波(a)窗函数设计低通滤波器程序设计如下：clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1;Fs=8000;wc=2*pi*fc/Fs;wp=2*pi*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7);N=ceil((As-8)/2.285/wdel);wn=kaiser(N+1,beta);ws=(wp+wc)/2/pi;b=fir1(N,ws,wn);figure(1);freqz(b,1);x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,1000,0,1.0]);title('滤波后信号频谱');subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot(2,2,4);plot(x);title('滤波前信号波形');sound(x,fs,bits);图形分析如下：（b）窗函数设计高通滤波器程序设计如下：clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fp=2800;fc=3000;As=100;Ap=1;Fs=8000;wc=2*pi*fc/Fs;wp=2*pi*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7);N=ceil((As-8)/2.285/wdel);wn=kaiser(N,beta);ws=(wp+wc)/2/pi;b=fir1(N-1,ws,'high',wn);figure(1);freqz(b,1);x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,1000,0,1.0]);title('滤波后信号频谱');subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot(2,2,4);plot(x);title('滤波前信号波形');sound(x,fs,bits);图形分析如下：(c)窗函数设计带通滤波器程序设计如下：clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fp1=1200;fp2=3000;fc1=1000;fc2=3200;As=100;Ap=1;Fs=8000;wp1=2*pi*fp1/Fs;wc1=2*pi*fc1/Fs;wp2=2*pi*fp2/Fs;wc2=2*pi*fc2/Fs;wdel=wp1-wc1;beta=0.112*(As-8.7);N=ceil((As-8)/2.285/wdel);ws=[(wp1+wc1)/2/pi,(wp2+wc2)/2/pi];wn=kaiser(N+1,beta);b=fir1(N,ws,wn);figure(1); freqz(b,1)x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,2000,0,0.0003]);title('滤波后信号频谱')subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot(2,2,4);plot(x);title('滤波前信号波形');sound(x,fs,bits);图形分析如下：(d)双线性变换法设计低通滤波器①选用butter程序设计如下：clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1;Fs=8000;wc=2*fc/Fs;wp=2*fp/Fs;[N,ws]=buttord(wc,wp,Ap,As);[b,a]=butter(N,ws);figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,1000,0,1.0]);title('滤波后信号频谱');subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot(2,2,4);plot(x);title('滤波前信号波形');sound(x,fs,bits);图形分析②选用cheby1程序设计如下：clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1;;Fs=8000;wc=2*fc/Fs;wb=2*fp/Fs;[n,wp]=cheb1ord(wc,wb,Ap,As);[b,a]=cheby1(n,Ap,wp);figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,4000,0,0.03]);title('滤波后信号频谱');subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot(2,2,4);plot(x);title('滤波前信号波形');sound(x,fs,bits);图形分析如下：(e)双线性变换法设计高通滤波器①选用butter程序设计如下：clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fc=2800;fp=3000;As=100;Ap=1;Fs=8000;wc=2*fc/Fs;wp=2*fp/Fs;[N,ws]=buttord(wc,wp,Ap,As);[b,a]=butter(N,ws,'high');figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,1000,0,1.0]);title('滤波后信号频谱');subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot(2,2,4);plot(x);title('滤波前信号波形');sound(x,fs,bits);图形分析如下：②选用cheby1程序设计如下：clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fc=2800;fp=3000;As=100;Ap=1;Fs=8000;wc=2*fc/Fs;wb=2*fp/Fs;[n,wp]=cheb1ord(wc,wb,Ap,As);[b,a]=cheby1(n,Ap,wp,'high');figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,4000,0,0.03]);title('滤波后信号频谱');subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot(2,2,4);plot(x);title('滤波前信号波形');sound(x,fs,bits);图形分析如下：(f)双线性变换法设计带通滤波器①选用butter程序设计如下clear;closeall[z1,fs,bits]=wavread('d:\laoA.wav')y1=z1(1:8192);Y1=fft(y1);fp1=1200;fp2=3000;fc1=1000;fc2=3200;As=100;Ap=1;Fs=8000;wc=[2*fc1/Fs,2*fc2/Fs];wp=[2*fp1/Fs,2*fp2/Fs];[N,ws]=buttord(wc,wp,Ap,As);[b,a]=butter(N,ws,'stop');figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]);title('滤波前信号频谱');subplot(2,2,2);plot(abs(X));axis([0,1000,0,1.0]);title('滤波后信号频谱');subplot(2,2,3);plot(z1);title('滤波前信号波形');subplot