基于MATLAB的PMSM直接转矩控制的建模方法的研究精品资料.doc
目 录1绪论12永磁同步电动机的原理23永磁同步电机的直接转矩控制33.1永磁同步电机的数学模型33.2永磁同步电机直接转矩控制系统的实现43.2.1转矩增量与定子电压空间矢量关系模型43.2.2定子磁链控制43.3逆变器开关时间控制模型63.4永磁同步电机直接转矩控制的系统74系统仿真模型的组建94.1仿真系统94.2其他模型的建立125仿真结果及其分析146结束语16参考文献16致谢181 绪 论直接转矩控制（DTC）是在空间矢量调速理论的基础上发展起来的一种新型交流电动机调速策略，其基本思想是根据交流电动机的转矩要求，直接选择合适的定子电压空间矢量，实现交流电动机电磁转矩的快速响应。由于直接在定子两相静止坐标系统下分析交流电动机的数学模型，将定子磁链与电磁转矩作为被控制量，根据给定转矩与实际转矩以及给定定子磁链与实际定子磁链的偏差来直接选择电压矢量,从而避免了矢量控制中许多复杂的矢量变换计算。所以直接转矩控制策略具有控制方式简单、转矩响应快、便于实现全数字化等优点。直接转矩控制在异步电动机调速系统中的应用已经比较成熟，但在永磁同步电动机（PMSM）伺服控制系统中的应用研究相对滞后。由于永磁同步电动机具有体积小、重量轻、运行可靠、功率密度高等诸多优点，将DTC控制策略应用于永磁同步电机控制中，以提高电机的快速转矩响应，成为研究者关注的课题究的热点课题。由于电机转矩和磁链的计算对控制系统性能影响较大，为了获得满意的转矩计算，仿真研究是最有效的工具和手段。本文中利用MATLAB软件下的Simulink仿真工具对PMSM DTC系统进行仿真；同时还详细地介绍了DTC系统中各控制计算单元的模型的建立，并分析控制系统的性能。2 永磁同步电动机的原理永磁同步电动机的转子旋转速度与定子绕组所产生的旋转磁场的速度是一样的。正由于这样，同步电动机的电流在相位上是超前于电压的，即同步电动机是一个容性负载。同步电动机在结构上大致有两种：永磁同步电动机的转子磁钢的几何形状不同，使得转子磁场在空间的分布可分为正弦波和梯形波两种。因此，当转子旋转时，在定子上产生的反电动势波形也有两种：一种为正弦波；另一种为梯形波。这样就造成两种同步电动机在原理、模型及控制方法上有所不同，为了区别由它们组成的永磁同步电动机交流调速系统，习惯上又把正弦波永磁同步电动机组成的调速系统称为正弦型永磁同步电动机（PMSM）调速系统；而由梯形波（方波）永磁同步电动机组成的调速系统，在原理和控制方法上与直流电动机系统类似，故称这种系统为无刷直流电动机（BLDCM）调速系统。永磁同步电动机转子磁路结构不同，则电动机的运行特性、控制系统等也不同。根据永磁体在转子上的位置的不同，永磁同步电动机主要可分为，表面式和内置式。在表面式永磁同步电动机中，永磁体通常呈瓦片形，并位于转子铁心的外表面上，这种电机的重要特点是直、交轴的主电感相等；而内置式永磁同步电机的永磁体位于转子内部，永磁体外表面与定子铁心内圆之间有铁磁物质制成的极靴，可以保护永磁体。这种永磁电机的重要特点是直、交轴的主电感不相等。因此这两种电机的性能有所不同。3 永磁同步电机的直接转矩控制3.1 永磁同步电机的数学模型直接转矩控制是对定转子磁链间的夹角也就是对转矩角进行控制。先推导转矩角与电机电磁转矩的数学关系。其中要用到与转子同步旋转的d-q坐标系，d轴正方向为转子磁链方向；与定子同步旋转的x-y坐标系，x轴正方向为定子磁链方向；两相-静止坐标系，轴正方向与电机a轴重合。如下图3-1所示。 图3-1 PMSM在不同坐标系下的矢量图在建立PMSM数学模型前，先对PMSM作如下假设：忽略电动机铁心的饱和，不计电动机中的涡流和磁滞损耗，转子无阻尼绕组，则永磁同步电机在d-q坐标系下的基本方程如下： （3-1） (3-2) (3-3) (3-4) (3-5) (3-6)式中：、为定子磁链d、q轴分量；、为定子绕组d、q轴等效电感；、为定子电流d、q轴分量；、为定子电压d、q轴分量；为转子磁链；为定子绕组电阻；P为微分算子；r为转子机械角速度；为电磁转矩；为电机极对数；为负载转矩；J为电机转动惯量；B为粘滞系数。经过进一步推导，得转矩控制的基本方程式： (3-7)对于隐极式永磁同步电机，lq=ld，故(37)式可写为： (3-8) 对于永磁同步电机，转子磁链不变，当定子磁链不变时，电机的电磁转矩只与转矩有关。由于定子磁链的状态比转子磁链的状态容易改变，因此可以通过控制定子磁链的状态来改变转矩角，从而实现对转矩的控制。3.2 永磁同步电机直接转矩控制系统的实现3.2.1 转矩增量与定子电压空间矢量关系模型与轴的夹角为，则在-轴系中，定子的合成空间磁链矢量以角速度 旋转，设逆时针方向为正方向。为了获得尽可能小的转矩脉动，定子的合成磁链空间矢量最好采圆形轨迹，参考磁链的幅值为常量，如图3-2所示。图3-2 转矩增量与定子电压空间矢量关系3.2.2 定子磁链控制在永磁同步电机中施加零电压矢量时，定子磁链停止转动，转矩几乎维持恒定。利用这一特点，可以采用具有零电压矢量的新型开关表实现直接转矩控制，把零电压矢量看成保持当前转矩的作用，利用零电压矢量来减小转矩的波动，以及减少逆变器开关次数和转矩的脉动。开关表中输出量与输入量的关系见表3-1，其中v(xxx)中x=1表示逆变器相应相的上臂导通，0表示下臂导通。的分布见图3-3。区域的划分是以-坐标系的轴为基准确定的，因此定子磁链的位置可用其-轴分量、估算的定子旋转磁链来确定。v1v6为逆变器输出电压的空间电压矢量,当施加电压矢量与夹角小于/2时，将使磁链幅值增加；当大于/2时,磁链幅值减小当电压矢量超前于时，转矩增加；落后于时，转矩减小。表3-1 逆变器开关表12345601v1（001）v2（001）v3（001）v4（001）v5（001）v6（001）0v0（000）v7（111）v0（000）v7（111）v0（000）v7（111）-1v1（001）v5（101）v4（110）v6（001）v2（010）v3（011）11v1（001）v2（010）v3（011）v1（001）v5（101）v4（110）0v7（111）v0（000）v7（111）v0（000）v7（111）v0（000）-1v5（101）v4（110）v6（001）v2（010）v3（011）v1（001）图3-3 空间电压矢量和分区用、分别表示电机磁链和转矩的给定值和实际值的误差状态，当给定值比实际值大时状态为1，否则状态为0，则由、的状态以及磁链所处分区的位置，便可按表1选择开关电压矢量，并为了方便在MATLAB/simulink中实现查表的位置我们假设一个变量的到： （3-9）3.3 逆变器开关时间控制模型如图2所示，设定当前采样周期的定子磁链矢量为，下一采样周期的定子磁链矢量为，所以给定转矩和估计转矩之间存在误差。将式(3-8)的定子磁链增量分解到-轴系的两坐标轴上, （3-10）再利用Clarke逆变换，将式(3-10)映射到定子三相坐标系A、B、C轴系上，把定子磁链增量分配到三相绕组，如： （3-11）定子合成磁链矢量的增量在各相绕组上的分量惟一地由参考电压在该绕组上的作用时间决定，参考电压在各相绕组上的作用时间是定子合成磁链矢量空间位置角的函数，如： （3-12）其中为定子合成磁链矢量相对于A相绕组（轴）的空间角位移。一个周期内定子参考电压矢量的作用时间根据式（3-8）和（3-9）求得，即 （1-13）上式中，虽然ds和都是空间矢量，但是二者方向相同，所以直接采用其幅值进行计算。是参考电压空间矢量的幅值为直流母线电压。3.4 永磁同步电机直接转矩控制的系统永磁同步电机直接转矩控制的系统结构如图3-4所示，是一个包含速度和转矩的双闭环系统。主要包括电机转速pi调节器、定子磁链和转矩观测器、转矩与磁链滞环比较控制器、最佳开关逻辑选择表、逆变器等。以和分别作为转矩和定子磁链的给定值，、分别为利用-坐标系的电压和电流估算得的转矩和磁链实际值，在调速系统中可作为转速调节器的输出，与实际转矩比较后经转矩滞环控制器输出转矩增减信号。当需要增大电磁转矩时，转矩控制器输出1，逆变器输出电压所形成的空间电压矢量使向前转动，由于电机的电磁时间常数小于机电时间常数，使定子磁链转速快于转子磁链转速，其结果是增大，增大了实际电磁转矩。转矩控制器输出0时，将减小实际转矩。进入稳态后，转矩给定值与实际值之差落在转矩控制器的滞环内，转速平均值也稳定为同步转速。为给定的定子磁链，与实际定子磁链相比较后，经磁链滞环控制器输出磁链增减信号，定子磁链保持在一个规定的范围内。这两个控制器的输出共同确定开关表的输出值，控制逆变器的pwm波输出。在控制过程中，在每一个采样周期开始，非零导通时间所对应的逆变器开关同时导通，向PMSM三相绕组供电，然后根据每相绕组的通电时间控制对应功率开关截止。一个周期内每个桥臂只有一个功率开关动作一次，具有恒定的开关频率，实际是一种恒频脉宽调制方式。该控制算法彻底避免了传统DTC控制以及SVM控制策略中磁链工作扇区的判别以及矢量合成等复杂的数据和逻辑计算，特别适合于数字PWM控制器件，可以直接用于DSP、FPGA等控制平台。图3-4 永磁同步电机直接转矩控制的系统结构4 系统仿真模型的组建4.1 仿真系统本文主要是利用Matlab软件进行仿真，Matlab的含义是矩阵实验室（MatrixLaboratory）。它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个方便的、界面友好的用户环境。在这个环境下，对所要求求解的问题，用户只需要简单地列出数学表达式，其结果便以数值或图形方式显示出来。Matlab的推出得到了各个领域专家学者的广泛关注，其强大的扩展功能更为各个工程领域提供了分析和设计的基础。Matlab包括被称作工具箱（Toolbox）的各类应用问题的求解工具。随着Matlab版本的不断升级，其所含的工具箱的功能也越来越丰富，因此应用范围也越来越广泛，成为涉及数值分析的各类设计不可或缺的工具。Simulink是基于Matlab的框图设计环境，可以用来对各种动态系统进行建模、分析和仿真，它的建模范围广泛，可以针对任何能用数学来描述的系统进行建模，例如航空航天动力学系统、卫星控制制导系统、通信系统、船舶及汽车等，其中包括了连续、离散，条件执行，事件驱动，单速率、多速率和混杂系统等。Simulink提供了利用鼠标拖放的方法来建立系统框图模型的图形界面，而且还提供了丰富的功能块以及不同的专业模块集合，利用Simulink几乎可以做到不书写一行代码即完成整个动态系统的建模工作。除此之外，Simulink还支持Stateflow，用来仿真事件驱动过程。Simulink是从底层开发的一个完整的仿真环境和图形界面，是模块化了的编程工具，它把Matlab的许多功能都设计成一个个直观的功能模块，把需要的功能模块用连线连起来就可以实现需要的仿真功能了。也可以根据自己的需要设计自己的功能模块，Simulink功能强大，界面友好是一种很不错的仿真工具。Simulink仿真具有以下的特点：（1）交互建模Simulink提供了大量的功能块，方便用户快速地建立动态系统模型，建模时只需要使用鼠标拖放库中的功能块，并将它们连接起来；用户可以通过将块组成子系统建立多级模型；对块和连接的数目没有限制。（2）交互仿真Simulink框图提供了交互性很强的非线性仿真环境。用户可以通过下拉菜单执行仿真，或者用命令行进行批处理。仿真结果可以在运行的同时通过示波器或者图形窗口显示。（3）能够扩充和定制Simulink的开放式结构允许用户扩充仿真环境的功能。（4）与Matlab和工具箱集成由于Simulink可以直接利用Matlab的数学、图形和编程功能，用户可以直接在Simulink下完成诸如数据分析、过程自动化、优化参数等工作。工具箱提供的高级设计和分析能力可以通过Simulink的屏蔽手段在仿真过程中执行。（5）专用模型库Simulink的模型库可以通过专用元件集进一步扩展。在MATLAB/simulink中建立了如图4-1的系统仿真模型。其中包括它包括3/2变换、磁链估算和转矩估算等子系统。进行磁链估算时，磁链初值不宜为0，否则仿真会出错。因此在磁链估算子系统中，要给积分模块赋一个初值。图4-1 基于MATLAB/simulink的PMSM DTC系统的仿真模型（6）区段判断的实现定子磁链矢量所在的区段我们可以根据磁链在坐标上的分量进行判定，由的正负确定定子磁链矢量的象限，再由和决定定子磁链矢量的具体位置，其实现模块如图4-1所示。其中的MATLAB函数模块是用来调用MATLAB中求反正切的函数，开关模块是一个2选1的输出，其输出再经过图4-2到图4-4子系统便可以得到区段结果。表4-1为磁链位置所对应的区段值。图4-2 角度计算图4-3 磁链区段的确定表4-1磁链区段和角度的关系角度区段角度区段-/2, -/66-/6, /61/6, /22/2,5 /635/6 ,7/64-/2,-/65图4-4 转矩和磁链误差信号4.2 其他模型的建立按照以上相类似的方法，我们对逆变器及其驱动信号、坐标的变换u和u的获取、电机磁链的估算和转矩的估算等等。建立相应的Simulink模型其子模块如图4-5到图4-8。图4-5 逆变器及其驱动信号图4-6 静止abc坐标系到坐标系的变换图4-7定子磁链估算模型图4-8 转矩估算模型5 仿真结果及其分析永磁同步电机的主要参数为：电阻rs=1.652，电感ld=lq=8.5mH，永磁磁链，极对数p=4，转动惯量j=0.0008kg·m2。系统给定转速n=500r/min，零负载启动。在t=0.5s突加4n·m的负载。仿真所得的转矩、转速、电流波形分别如图5-1到图5-3。在逆变器和PMSM子模块间，接入电压测量装置以观测的相间电压，因为当Simulink模块与PSB模块相连时，要求接入一个电气测量模块，否则仿真会出现错误。图5-1 转矩波形图5-2 转速波形图5-3 电流波形可以看出电机在通电以后，迅速到达最大转矩，然后很快回到稳定值，在0.5s时，负载转矩突变为4n·m，输出转矩波形跟随给定，实现了对电机的瞬时力矩的控制。转速迅速达到给定值500rad/s，在0.5s负载转矩发生变化时，转速有微小的震荡后回到给定值；从电流波形可看出：三相电流按正弦规律变化，说明电机的运行是正常的，在0.5s定子电流发生变化，可以看出此控制系统响应快，稳定性好，抗干扰性能高，仿真波形与理论分析情况一致。6 结束语与传统的建模方法相比采用MATLAB的PMSM直接转矩控制的建模方法可以在simulink环境下建立系统的仿真模型并进行仿真实现，仿真结果表明：波形符合理论分析，系统运行平稳，具有良好得静、动态性能。这为进一步分析和设计永磁同步电机直接转矩控制系统提供了有效的手段和工具。参考文献1 田淳，胡育文永磁同步电机直接转矩控制理论及控制方案的研究电工技术学报，2002（2）：8112 王沫然Simulink 4建模及动态仿真北京：电子工业出版社，20023 姜平CNC系统中指数升降速控制曲线的算法微特电机，2003，31（5）：16174 范治田新型数字化可编程频率合成器DDS现代电子技术2003，（12）：26305 周俊峰，陈涛基于FPGA的直接数字频率合成器的设计和实现国外电子元器件，2003，（1）：466 唐任远现代永磁电机理论与设计北京：机械工业出版社，19977 王成元矢量控制交流伺服取得电动机北京：机械工业出版社，19958 万文斌，徐衍亮，唐任远。永磁同步电动机的高性能电流控制器J.中国电机工程学报，2000，（12）：24-27.9王江，王家军，许镇林。基于逆变器死区特性的永磁同步电动机系统的自适应变结构控制J.中国电机工程学报，2000，21（8）：37-41.10 周杨忠，胡育文，田蕉。永磁同步电机控制系统中变比例系数转矩调节器设计究J。中国电机工程学报，2004，24（9）：204-208.11孙笑辉，张增科，韩曾晋。基于直接转矩控制的感应电动机转矩脉动最小化方法研究J。中国电机工程学报，2002，22（8）：110-112.12 杨家强，黄进。异步电动机直接转矩控制转矩脉动最小化方法研究J。电工技术学报，2004，19（9）：23-29.13 刘军，楚小刚，白华煜基于参考磁链电压空间矢量调制策略的永磁同步电机直接转矩控制研究J。电工技术学报，2005，20（6）：11-15.14 刘军，刘丁，吴浦升，等。基于模糊控制调节电压矢量作用时间策略的永磁同步电机直接转矩控制仿真研究J。中国电机工程学报，2004，24（10）：148-152.15 贾洪平，贺益康。永磁同步电机滑模变结构直接转矩控制J。电工技术学报，2006，21（1）：1-616 Takahashi I, Noguchi T.A new quick response and higlr effreciency control strategy of an motorJ. IEEE Transactions on Industrial Applications, 1986,22(5):821-827.17 Kin JH,CHoi JW. Novel rotor Flux Observer Characteristic Function in Complex Vector Space for Field Oriented Induction motor Drives J.IEEE. Trans Ind. Applicat,20092,38,(5):1334 - 134318 Hay lock JA, M edrow BC .Enhanced Current Control of High- Speed PM Machine Drives Through the Use of Flux Controllers J, IEEE, Trans Ind A ppicat,1999,35(5):1030-103819 Peter Vas. Vector control of ac machinesM. Oxford: Clarendon Press,1990.20 Zhong L, Rahman M F.Analysis of direct torque control in Permanent magnet synchronous motor drivesJ.IEEE Transactions on Power Electronics, 1997,12(3):528-535.21 Rahman M F, Zhong L .Comparison of torque responses of the Interior permanent magnet synchronous motor under PWM current and direct torque controlsC.IECON99 Proceedings,San Jose, CA,1999.22 Tiitinen P, Surandra M, The next generation motor control method,DTC direct torque controlC. Proceedings of the 1996 International Conference on Power Electronics, Drives and Energy Systems for Industrial Growth. IEEE. Part. 1. New York, USA, 1996:37-43致 谢本论文最终得以顺利完成，非常感谢我的指导教师宋教授。从论文选题直到论文的最终完成，她都给予我尽心尽力的指导。宋教授严谨的治学态度深深的影响着我，对我今后的学习、工作、生活必将产生影响。借此机会特向宋教授表示最诚挚的感谢。感谢平顶山学院学院的所有领导和老师。你们严谨的学风、渊博的知识、诲人不倦的品格一直感染和激励着我不断上进，使我大学的时光充实而有意义。“海纳百川，取则行远”，在这所美丽的校园里，不断成长，在这里我所学到的，必将使我受益终生。在本论文的写作中，我也参照了大量的著作和文章，许多学者的科研成果及写作思路给我很大启发，在此向这些学者们表示由衷的感谢。本论文虽然几经修改，但由于才疏学浅，本文疏漏之处在所难免，还望各位老师批评指正。附录资料：不需要的可以自行删除各类滤波器的MATLAB程序一、 理想低通滤波器IA=imread('lena.bmp');f1,f2=freqspace(size(IA),'meshgrid');Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r>0.2)=0;Y=fft2(double(IA);Y=fftshift(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(r<0.2)=0;Y=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'); 三、 Butterworth低通滤波器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=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'); 四、 Butterworth高通滤波器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(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);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(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(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','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=zeros(size(IA);Hd(r<D0)=1;for i=1:size(IA,1) for j=1:size(IA,2) if r(i,j)>=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(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(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=ones(size(IA);Hd(r<D1)=0;for i=1:size(IA,1) for j=1:size(IA,2) if r(i,j)>=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'); 九、 用其他方法编写的理想低通、理想高通、Butterworth低通、同态滤波程序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); if d<=d0 h=0; else h=1; 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); 3、 Butterworth低通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=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(
