DTMF电话中的拨号信号发生与检测系统仿真精品资料.doc

《DTMF电话中的拨号信号发生与检测系统仿真精品资料.doc》由会员分享，可在线阅读，更多相关《DTMF电话中的拨号信号发生与检测系统仿真精品资料.doc（28页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、DTMF电话中的拨号信号发生与检测系统仿真摘 要双音多频（DTMF）信号是音频电话中的拨号信号，由美国AT&T贝尔公司实验室研制，并用于电话网络中。这种信号制式具有很高的拨号速度，且容易自动检测识别，很快就代替了原有的用脉冲计数方式的拨号制式。这种双音多频信号制式不仅用在电话网络中，还可以用于传输十进制数据的其他通信系统中，用于电子邮件和银行系统中。这些系统中用户可以用电话发送DTMF信号选择语音菜单进行操作。 DTMF信号系统是一个典型的小型信号处理系统，它要用数字方法产生模拟信号并进行传输，其中还用到了D/A变换器；在接收端用A/D变换器将其转换成数字信号，并进行数字信号处理。本文着重介绍

2、了一种快速的频率检测算法GEORTZEL算法的原理，以及该算法在双音多频拨号检测系统中的应用。并根据此算法设计出了相应的滤波器，利用MATLAB在计算机上对双音多频的信令产生及检测接收系统进行了仿真测试。关键词：双音多频，Goertzel算法，数字信号处理，傅里叶变换；目 录摘要 11.绪论31.1 引言31.2 数字信号简介31.3 数字滤波器41.4 频率分析4 1.4.1 傅里叶级数及傅里叶变换 5 1.4.2 Z变换及系统函数5 1.4.3 离散系统傅里叶变换（DFT）62.双音多频（DTMF）通信系统6 2.1 双音多频（DTMF）信号简介6 2.2 双音多频（DTMF）信号的产生7

3、 2.3 双音多频（DTMF）信号的检测83.戈泽尔(Goertzel)算法 103.1 戈泽尔(Goertzel)算法原理 103.2 戈泽尔(Goertzel)算法程序 125.双音多频（DTMF）系统的仿真136.结束语 15参考文献 16 1 绪论1.1 引言双音多频（DTMF）信号，由美国贝尔实验室研制，用于音频电话网络中的拨号信号。一方面这种信号具有非常高的拨号速度，另一方面它便于自动检测识别及电话业务的拓展，所以很快代替了原有的脉冲计数方式的拨号形式。自1876年美国贝尔发明电话以来，电话交换技术发展主要经历了三个阶段：人工交换、机电交换、电子交换。早在1878年就出现了人工供电

4、制交换机，它借助话务员进行话务连接。15年后步进制的交换机问世，这种交换机属于“直接控制”方式，用户通过话机拨号脉冲直接控制步进接续器做升降和旋转动作。从而自动完成用户间的接续。这种交换机虽然实现了自动接续，但存在着速度慢、效率低、杂音大与机械磨损严重等缺点。直到1938年发明了纵横制交换机由直接控制过渡到间接控制方式，随后又出现了电子交换机及程控交换机等。用户的拨号脉冲不再直接控制接线器动作，而先由记发器接收、存储，然后通过标志器驱动接线器，以完成用户间接续。由于标志器记取的是标识信号，所以可以不采用不稳定的脉冲计数方式，逐步发展到数字号码标识，所以快速、稳定灵活的DTMF音频拨号方式取代了

5、脉冲计数方式。1.2 数字信号简介信号是反映消息的物理量，是消息的表现形式。由于非电的物理量可以通过各种传感器较容易地转换成电信号，而电信号又容易传送和控制，所以使其成为应用最广的信号。电信号是指随着时间而变化的电压或电流，因此在数学描述上可将它表示为时间的函数，并可画出其波形。信息可以通过电信号进行传送、交换、存储、提取等。较之传统的方法，数字信号处理系统自身有着无可比拟的优点：（1）便于加密处理（2）便于存储、处理和交换（3）设备便于集成化、微型（4）便于构成综合数字网和综合业务数字网（5）占用信道频带较宽一路模拟电话的频带为4kHz带宽，一路数字电话约占64kHz，这是模拟通信目前仍有生

6、命力的主要原因。1.3 数字滤波器数字滤波器由数字加法器、单位延时和常数乘法器组成的一种算法或装置。数字滤波器的功能是对输入离散信号的数字代码进行运算处理，即把输入序列通过一定的运算变换成输出序列，以达到改变信号频谱的目的。数学公示表示如下： （1.3-1)通常可以用两种方法实现数字滤波器：一种方法是把滤波器所完成的运算编成程序并让计算机执行，也就是采用就算机软件实现；另一种方法是设计专用的数字硬件、专用的数字信号处理器或采用通用的数字信号处理器来实现。存在着两种数字滤波器的的基本类型：无限长单位冲击响应（IIR）滤波器和有限长单位冲击响应（FIR）滤波器。它们的系统函数分别为： (1.3-2

7、) (1.3-3)式（1.3-2）中的H(Z)称为N阶IIR滤波器函数式，（1.3-3）中的H(Z)称为N-1阶FIR滤波器函数。无限长单位冲激响应（IIR）滤波器的单位冲激响应h（n）是无限长的；系统函数H(Z)在有限Z平面（0Z）上有极点存在；结构上存在着输出到输入的反馈，也就是结构上是递归型的。1.4 频率分析 任何给定信号的频率分析，都包含着时域信号向其频率分量的转换。之所以要在频域内对信号进行描述，因为在实现信号处理时，通常都利用以频率响应形式描述的系统。而这些频域信号处理概念时傅立叶变换是必不可少的。同时Z变换对进行数字系统分析和实现也极其重要。1.4.1 傅立叶级数及傅立叶变换任

8、何周期信号X(t)均可以表示为无穷多个谐波调谐相关的正弦波和复指数之和。具有周期T0(秒)周期信号x(t)，其基本的数学表达式为具有系列定义的傅立叶级数： (1.4.1-1)式中Ck为傅立叶级系数，第k个傅立叶系数Ck可以表示为 (1.4.1-2)用周期信号的傅立叶级数通过极限的方法导出的非周期信号的频谱表示式，称为傅立叶变换：傅立叶逆变换： (1.4.1-3)傅立叶正变换： = (1.4.1-4)1.4.2 Z变换及系统函数对于离散时间系统来说，与拉普拉斯对应的变换是是Z变换，Z变换提供了离散时间信号和系统的频域描述，并为数字滤波器的设计和实现，提供了有力的工具。序列x(n)的Z变换定义为

9、： (1.4.2-1)式中z是一个复变量，定义中对n求和是在和之间求和，称为双边Z变换。还有一种称为单边Z变换的定义，如下式： (1.4.2-2)将系统的单位脉冲响应h（n）进行Z变换，得到H（z）,一般称H（z）为系统的系统函数，它表征了系统的频域特性。对N阶差分方程进行Z变换，得到系统函数的一般表示式 (1.4.2-3)令Z=ej即得系统的传输函数H(ej)。1.4.3 离散系统傅立叶变换（DFT）离散傅立叶变换是针对有限长序列或周期序列从存在的，相当把序列的连续傅立叶变换加以抽样，频率的离散化造成时间函数也呈周期，故级数应限制在一个周期内。令，则，我们即可得到离散傅立叶变换对为： (1.

10、4.31) (1.4.32)其中表示有限长序列的抽样点数，或周期序列一个周期的抽样点数。2 双音多频（DTMF）通信系统2.1 双音多频（DTMF）信号过去的电话拨号时靠脉冲计数确定09这10个数字的，拨号速度慢，也不能扩展电话上网其他服务功能。现在均采用双音拨号。每一位号码由两个不同的单音频组成，所有的频率可分成高频带和低频带两组。每一位号码由一个高频信号和一个低频信号叠加形成，DTMF信号有16个编码。一般情况下，声音信号很难造成对DTMF接收器的错误触发。电话中的双音多频信号有两个作用：用拨号信号去控制交换机接通被叫的用户电话机；控制电话机的各种动作，如播放留言、语音信箱等。2.2 双音

11、多频（DTMF）信号的产生DTMF编码器基于两个二阶数字正弦波振荡器，一个用于产生行频，一个用于产生列频。典型的DTMF信号频率范围是7001700Hz，选取8000Hz作为采样频率，DTMF双音频信号由两个二阶数字正弦振荡器产生，一个用来产生行音频信号，另一个产生列音频信号。在电话键盘上的每次按键，会产生式（2.2-1）表示的两个音调之和。规定用8KHz对DTMF信号进行采样，采样后得到时域离散信号为 (2.2-1)两个频率f1,f2唯一确定了被按压的键。具体数字对应的两个频率如表2.2.1所示。表2-1 DTMF拨号频率（第四列常被省略）频低频高1209Hz1336Hz1477Hz1633

12、Hz697Hz123A770Hz456B852Hz789C942Hz*0#D形成上面序列的方法有两种，即计算法和查表法。用计算法求正弦波的序列值容易，但实际中要占用一些计算时间，影响运行速度。查表法是预先将正弦波的各序列值计算出来，寄存在存储器中，运行时只要按顺序和一定的速度取出便可。这种方法要占用一定的存储空间，但是速度快。因为采样频率是8000Hz，因此要求每125ms输出一个样本，得到的序列再送到D/A变换器和平滑滤波器，输出便是连续时间的DTMF信号。DTMF信号通过电话线路送到交换机。2.3双音多频（DTMF）信号的检测在接收端，要对收到的将收到的双音多频信号进行检测，即检测两个正弦

13、波的频率，以判断其对应的十进制数字或者符号。用数字方法进行检测，需要将接收到的模拟音频信号进行A/D变换，恢复为数字信号，然后检测其中的音频频谱来确定所发送的数字。检测方法有两种，一种是用一组滤波器来提取所需频率，判断对应的数字或符号；另一种是用快速傅立叶变换（FFT）算法的DFT对双音多频信号进行频谱分析，有信号的幅度谱，判断信号的两个频率，最后确定对应的数字或符号。以下章节对两种方法进行分析并对Goertzel算法进行详细的推导。由上面可以知道一个DTMF信号是由两个不通频率f1和f2的正弦波组成，它可以用下式表示 (2.3-1)我们可以考虑通过离散傅立叶变换进行信号的频谱分析来检测离散的

14、双音多频码。这是因为只有时域及频域都是离散的情况下，才能适合于在计算机上运算；也就是周期的离散时间信号与离散频率间的变换对。周期性离散时间信号x（n）的离散傅立叶变换周期性离散频率函数X（k）两相邻谱线分量之间的角频率增量与周期之间的关系可表示为， 代表信号的基频。取样频率与取样周期T的关系是 取样角频率 T为时域取样间隔，在一个周期内取样点数为N。在自变量为t及f的情况下，在一个时域中对函数进行取样，两取样点间增量的倒数，必是另一个域中函数的周期。现序列的周期为NT，所以对频谱取样的间距是。以数字频率表示时，则频谱间距是 k0，1，N1令并称之为因子。离散傅立叶变换（DFT）可写成如下形式：

15、， 0 kN1 (2.3-2)DFT的表达式可以直接写成：， k0，1，N1 ( 2.3-3) 一般来说，x（n）和都是复数，X（k）也是复数，每计算一个X（k）值，需要N次复数乘法(x（n）与相乘)以及(N1)次复数加法。而X（k）一共有N个点（k从0取到N1），所以完成整个DFT运算总共需要次复数乘法及N（N1）次复数加法。复数运算实际上是由实数运算来完成的，因此上式可写成： (2.3-4)由此可见，一次复数乘法需用四次实数乘法和两次实数加法；一次复数加法则需两次复数加法。因而每运算一个需要4N次复数乘法和2N2（N-1）=2(2N-1)次实数加法。 所以，直接计算DFT，乘法次数和加法次

16、数都是和成正比的，当N很大时，运算量是很可观的，例如，当N8时，DFT需64次复数乘法，而当N1024时，DFT所需复数乘法为1，048，576次，即一百多万次复数乘法运算，这对实时性很强的信号处理来说，对计算速度的要求太高了。因而需要改进对DFT的计算方法，以大大减少运算次数。仔细观察DFT的运算量可以看出，利用系数的以下固有特性，可以减少DFT的运算量。(1)的共轭对称性： (2)的周期性： (3) 的可约性：由此得出，这样，（1）利用这些特性，使DFT运算中有些项可以合并；（2）利用的周期性和对称性，可以将长序列的DFT分解为短序列的DFT。而前面已经说到，DFT的运算量使与成正比的，所

17、以N越小越有利，因而小点数的DFT比大点数的DFT的运算量要小。基于这样的思路，快速傅立叶变换（FFT）发展起来了。如果按时间抽取的FFT算法，系统的运算量是，运算量明显比DFT运算量减少。但它有一个缺点，需要全部抽样数据到达时，才能进行计算；不适合双音多频信号频谱的计算，满足不了双音拨号系统实时性的要求。因此开发出了适合于DMFT的特殊算法Goertzel算法。3 戈泽尔(Goertzel)算法3.1 戈泽尔(Goertzel)算法原理戈泽尔算法利用DFT中的旋转因子的周期性，将DFT的运算转换成一种线性滤波运算。下面推导戈泽尔算法的计算公式和实现结构。假设长度为N的序列x（n）的N点DFT

18、用X（k）表示，因为 （3.1-1）因此 k=0,1,12,N-1 (3.1-2)按照上式定义序列 （3.1-3）这个方程可以解释为有限长序列x(n)，与序列的卷积，因此表示为 （3.1-4）令 （3.1-5）则 （3.1-6）由上式，将看成是序列x（n）通过单位脉冲响应为的滤波器的输出，根据式（3.1-2）和式（3.1-3），得到 (3.1-7)那么，x（n）的DFT的第k点就是序列x(n)通过滤波器输出的第n=N点样值。那么N点DFT就是这N个滤波器分别对序列x(n)的响应序列的第N点输出。对式（3.1-5）进行Z变换，得到滤波器的系统函数 （3.1-8）该滤波器是一个一阶纯极点滤波器，极

19、点为，极点频率为。该滤波器的结构图如图3-1所示，戈泽尔算法的原理方框图如图3-2所示。图3-1单极点滤波器图3-2戈泽尔算法的滤波器组在图3-1中存在一次复乘算法，为了避免复乘算法，将一阶纯极点滤波器变为二阶滤波器，推导如下 （3.1-9）按照图3-2所示。再按照该结构图，可以用两个差分方程表示该二阶滤波器，即 （3.1-10）= （3.1-11）对于音调检测，只有的幅值平方是需要的，通过实现这种需求，可以使算法作进一步简化，由（3.1-11）式，的平方幅值可以计算为： （3.1-12）这样输入信号是实序列，用式（3.1-10）计算中间变量，用式（3.1-12）计算输出信号的幅度，避免了复数

20、的运算。由此得到=。因为有8种音频要检测，所以需要8个式（3.1-8）表示的滤波器，或者8个式（3.1-9）表示的滤波器。8个滤波器的中心频率分别对应8种音频。3.2 Goertzel算法程序下面是按（3.1-12）式编写的Goertzel算法子程序。它根据输入序列x（n）和指定的DFT样本的序号k计算待求的DFT样本X。其核心语句只有四句，具体如下：function X=gfft(x,k)N=length(x);x1=x,0;d1=2*cos(2*pi*k/N)v=filter(1,1,-d1,1,x1);w=exp(-i*2*pi*k/N);x=v(N+1)-w*v(N);5 双音多频（D

21、TMF）系统的仿真DTMF信号的产生与识别仿真实验在MATLAB环境下进行，编写仿真程序，运行程序，送入8位电话号码，程序自动产生每一位号码数字相应的DTMF信号，并送出双频声音，再用DFT进行谱分析，显示每一位号码数字的DTMF信号的DFT幅度谱，安照幅度谱的最大值确定对应的频率，再安照频率确定每一位对应的号码数字，最后输出16位电话号码，程序名为exp8。下面用MATLAB程序对该系统进行模拟，程序分四段：第一段（第27行）设置参数，并读入8位电话号码；第二段（第920行）根据输入的8位电话号码产生时域离散DTMF信号，并连续发出8位号码对应的双音频声音；第三段（第2225行）对时域离散D

22、TMF信号进行频率检测，画出幅度谱；第四段（第2633行）根据幅度谱的两个峰值，分别查找并确定输入的8位电话号码。程序如下：% DTMF双频拨号信号的生成程序：%clear all;clc;tm=1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68; % DTMF信号代表的16个数N=205;K=18,20,22,24,31,34,38,42;%16个基频对应的16个k值f1=697,770,852,941; % 行频率向量f2=1209,1336,1477,1633; % 列频率向量TN=input(键入16位电话号码= ); % 输入16位数字TNr=0; %接收端电

23、话号码初值为零for l=1:8; %分别对每位号码数字处理：产生信号，发声，检测 d=fix(TN/10(8-l); %计算出第l位号码数字 TN=TN-d*10(8-l); for p=1:4; for q=1:4; if tm(p,q)=abs(d); break,end % 检测与l位号码相符的列号q endif tm(p,q)=abs(d); break,end % 检测与l位号码相符的行号pEnd n=0:1023; % 为了发声，加长序列 x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 构成双音频信号 sound(x,8

24、000); % 发出声音 pause(0.1) %相邻信号响声之间加0.1秒停顿 % 接收检测端的程序 X=goertzel(x(1:N),K+1); % 用Goertzel算法计算16点DFT样本 val = abs(X); % 列出八点DFT的摸 subplot(4,2,l); stem(K,val,.);grid;xlabel(k);ylabel(|X(k)|) % 画出16点DFT(k)的幅度 axis(10 50 0 120) limit = 80; %基频检测门限为80 for s=5:8; if val(s) limit, break, end % 查找列号 end for r=

25、1:4; if val(r) limit, break, end % 查找行号end18 TNr=TNr+tm(r,s-4)*10(16-l); %将16位号码表示成一个16位数，以便显示enddisp(接收端检测到的号码为：) % 显示接收到的16位电话号码disp(TNr)编译程序，弹出操作界面，点击操作界面的数字按键，程序开始运行，数据经连接函数与接收端程序连接输入键值，经Goertzel算法的处理将频谱叠加，在频谱图显示窗口显示输入频率（axes2）和叠加频率(axes1),同时接收端程序检根据合成频率检测出键值，在editor窗口显示检测到输入的数字。6 结束语本文采用GEORTZE

26、L算法，实现了对双音多频信号的检测，在MATLAB环境下，对双音多频信号的解码检测过程进行了模拟仿真。本文只是完成了双音多频拨号系统的软件仿真，如何实现高性价比的硬件系统是下一步需要研究和完成的工作。在做MATLAB界面，做好以后同时会产一个M文件，只要把代码编写在对应的控件中即可以方便的实现仿真。试验结果表明基于MATLAB 的DTMF 编解码具有速度快、精度高、稳定性好以及便于灵活应用等特点。而利用MATLAB语言，可以较好地、快速地进行算法的仿真。参考文献1 高西全、丁玉美、阔永红.数字信号处理原理、实现及应用.电子工业出版社,20062 程佩青.数字信号处理教程.清华大学出版社,200

27、1年, 第二版3 Joyce van de Vegte,侯正信,王国安译.数字信号处理基础.电子科工业出版社,2003年4 赵霞.双音多频信号产生及解码研究.微电子学,2001年,第6期5 陈怀琛.数字信号处理教程MATLAB释疑与实现.电子工业出版社,2004附录资料：不需要的可以自行删除各类滤波器的MATLAB程序一、 理想低通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(

28、Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 二、理想高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(doub

29、le(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 三、 Butterworth低通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;fo

30、r i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 四、 Butterwo

31、rth高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r(i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(u

32、int8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 五、 高斯低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double

33、(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 六、 高斯高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(

34、IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facel

35、ighting,phong); 七、 梯形低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=zeros(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft

36、2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 八、 梯形高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=ones(size(IA);Hd

37、(r=D0 & r(i,j)=D1 Hd(i,j)=(D0-r(i,j)/(D0-D1); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 九、 用其他方法编写的理想低通、理想高通、Butterwort

38、h低通、同态滤波程序1、 理想低通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d=d0 h=1; else h=0; end y(i,j)=g(i,j)*h; endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 2、 理想高通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);n=2;d0=10;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2