如何用matlab进行滤波器设计和仿真.doc

华 中 科 技 大 学本科课程考试答题本 考生姓名 考生学号 班 级 联系电话 Email： 考试科目 考试日期 一、报告摘要： 本报告首先论述如何用matlab进行滤波器设计和仿真，并且将参数转换成DSP可用的滤波器参数。然后，将这组参数在CCS上进行纯软件仿真。随之是本次实验过程中代码拼接的问题与解决方法。软硬件流程图及核心代码的陈述。将实际滤波的结果进行展示。与理想滤波器的比较之间的差异与改良。最后是对于本课程的感想二、滤波器基础知识以及基于matlab的参数设计：设计一个能够将1000Hz一下保留，1500Hz以上的频率滤出的滤波器，利用理想的fir滤波器截止频率可以设定为（1000+1500）/2=1250Hz，利用哈明窗设计本滤波器，阶数为21阶，functiona=ccc()clear all;clf;N=1024;%设置采样点数fs=6000;%设置采样频率dt=1/fs;for k=1:N; f1=1500; f2=1000;y(k)=sin(2*pi*f1*k*dt)+sin(2*pi*f2*k*dt)+sin(2*pi*(f1+200)*k*dt);%+sin(2*pi*1250*k*dt);endwp=pi*1000*2/6000;%归一化的通带频率ws=pi*1500*2/6000;%归一化的阻带频率wc=(wp+ws)/2; %理想的滤波器的截止频率M=21; %设置阶数b=fir1(M-1,wc/pi,hamming(M)%;利用fir1函数得到滤波系数bH,F=freqz(b,1,N); %figure(1);plot(F/pi,20*log10(abs(H);grid;xlabel(f(Hz);ylabel(幅度);%画出幅频特性yy1=conv(b,y); %做卷积即为滤波y=fft(y,N); %做fft运算得到频域特性pyy=y.*conj(y); %乘以共轭，得到的是y的幅值的平方f=(0:(N/2-1);figure(2);plot(f*fs/N,pyy(1:N/2); y=fft(yy1,N);pyy=y.*conj(y);f=(0:(N/2-1);figure(3);plot(f*fs/N,pyy(1:N/2);执行以上函数后得到以下图像：可以看到在（2*1000）/6000=0.33处的衰减基本不变，而在（1000+1500）/6000=0.42处衰减值已经大于3dB了，而在（2*1500）/6000=0.5衰减值已经接近20dB了。显然本滤波器与理想的滤波器有些差异，当然随着阶数的提高会更加接近理想的滤波器。滤波前后的变化，可以清楚地看到1000Hz处的幅值并没有什么衰减，而在1500Hz以上的分量基本全部被滤除（在1500Hz只衰减不到20dB）。滤波系数：A0A1A2A3A4A5A6A7A8A9A100.0013 -0.0026-0.00580.00320.02110.0089-0.0470-0.06080.07260.30070.4169A11A12A13A14A15A16A17A18A19A200.30070.07260.0608-0.04700.00890.02110.0032-0.0058-0.00260.0013转换后得到的滤波系数如下：a0 002AH a1 -0054H a2 -00BEH a3 0068H a4 02B4H a5 0124H a6 -0605H a7 -07C9H a8 094CH a9 267DH a10 355CH a11 267DH ;a12 094CH a13 -07C9H a14 -0605H a15 0124H a16 02B4H a17 0068H a18 -00BEH a19 -0054H a20 002AH三、算法验证：将matlab得到的的滤波器系数转换成十六进制数，并设计一个21阶的fir滤波器在ccs2.0中进行仿真。STM#K_CIR,BK;设置循环区STM#1,AR0 ;设置便宜的步长STM#inputdata,ORIGINSTM#bufferdatax,INPUTSTM#filterdata,OUTPUTRPT #K_A-2 MVDD *ORIGIN+, *INPUT+0%；放前20个数据，并将INPUT指向下一位，同时将ORIGIN指向下一位， STM #K_DATA_SIZE-K_A-1,BRC RPTB filter_end-1 MVDD *ORIGIN+, *INPUT；在第一次滤波时将第21个数据放入，并且在之后的每次滤波中均放进一个数据，用于本次滤波 RPT #K_A-2 MAR *INPUT-0% ；将INPUT左移动20个，第一次滤波时正好将INPUT指向x(1)，就能保证滤波的进行， MPY *INPUT+0%, #a20,B LD B, A ；将该数存进累加器A， MPY *INPUT+0%, #a19,BADD B, A该段为a18至a1与对应的值相乘，与a19的操作完全相同，为简化文档，省去。 MPY *INPUT+0%, #a0,B ADD B, A STH A, *OUTPUT+；将高十六位放进滤波后数据的存储空间，并指向下个空间接收下次滤波数据本过程核心在于循环缓冲区和循环块设置，因为y(n)=a0*x(n)+a1*x(n-1)+a20*x(n-20)，本实验现放20个数据，也就是x(n)的前20个，此过程中INPUT都将在存数阶数后指向下个存储空间，随后放入第21个数，不移动，在滤波前将INPUT左移20个空间，就能够充分保证上述式子的正确执行，所以本实验中将y(21)作为第一个滤波数据，这样也自然会导致只能滤出 K_DATA_SIZE-K_A个数据，这也正是将循环块大小定为K_DATA_SIZE-K_A-1的原因。这些正是滤波程序的核心部分，加上其他的配置，就可以完成滤波了。*值得注意的是，本段只是将滤波器进行简单的验证以判断通过matlab得到的fir滤波器系数能否正确的进行滤波。与本实验的真正的代码还是有些许差别。经过ccs2.0验证得到的图5如下所示，可以看出能够满足低通的要求但是原来在0.41处的截止频率变为了现在的0.333，也就是有一定的误差。但是已经达到了滤波的能力。波形如下： 四、代码合成 由于在本实验开始之前就已经学习过ccs滤波器的仿真，其中包括有IIR和FIR滤波器的仿真，同时还学习了A/D、D/A转换和中断逻辑的联合使用。这些都在实验之前便已经完成了。 在本实验中，要完成的目标是能够将信源的信号通过DSP处理进行滤波，所以要有A/D和D/A的转换，而且数据通过采样获得，必然有定时中断。整个具体流程将在第五部分详述。 代码拼接，核心在于将滤波程序加在A/D的后面，将滤波得到的数据进行D/A转换就可以了。因为在之前已经做好了通过定时器实现中断控制将数据A/D、D/A转换的部分，所以此次试验不用修改终端服务程序的部分。这样就完成了粗略的合成，具体的还有很多修改。1） 因为本过程是对采样的一个数立即进行A/D，卷积，D/A这些操作，只用各定义一个用来存放A/D和D/A的存储空间即可。代码中将in_data定义为存储A/D后的数据存储空间，out_data为D/A后的数据存储空间。所以不用计数，也不用在WAIT中设置条件跳转。2） 将滤波段程序加在A/D的后面，并将存放20个0的初始化操作放在第一次进定时器中断之前，第21个数据在第1次定时中断A/D后放入，之后的每个数据都是在A/D放入。3） 依然是对PMST寄存器的高9位操作最后放入PMST中。4） 开定时器中断，置数0008H，放入IMR寄存器，并允许中断，并服务程序中的代码对应。5） A/D，D/A的初始化、A/D、D/A的代码，这部分不用修改，尤其要意端口地址。6） 程序有浮点数运算，必须加上SSBX FRCT。7） 内存分配，为了方便将涉及滤波的数据放在filter_vars的数据段中。在CMD文件也要加如对应的定义。8） 完成整个程序，具体实现将在第五部分详述。关于本实验遇到的问题:1） 在刚开始的时候由于对整体的过程没有足够清楚地了解，认为只用在A/D之后加上滤波断就可以完成，想通过原有的结构，放进256个采样点后进行滤波，以失败告终，因为这样不仅逻辑上比较复杂，而且效率低，采样控制更加复杂，可行性不高，最终弄清原理后，个分配一个空间即可，在代码中将in_data定义为存储A/D后的数据存储空间，out_data为D/A后的数据存储空间。2） 溢出错误，在刚开始设计时忽略了本实验系统的特性，A/D是10位数据线，而DSP处理器内部的寄存器是16位的，D/A是12位数据线，三者之间传递数据如果加处理可能出错，A/D的数据传入后，高六位并没有做出要求，在滤波，D/A这些过程也会产生对应的问题。本实验过程将A/D后的数据左移6位，并将滤波后的数据右移4位。3） 数据溢出错误，与上面的溢出错误有本质差别，上面的错误在于数据位数不同造成的，此处的错误是因为数据为负值时，被D/A后一然后会作为一个正值输出，此时在示波器上原本为负值位置上会产生一个巨大的正值波形。解决方法是加上直流电平，防止输出数据为负值。五、软硬件流程软件流程：1） 定时器： 定时器状态寄存器：寄存器0-15位总共16位。第三位是时钟中断屏蔽位。要将此位屏蔽关闭。 定时器硬件构成： PRD是一个4位寄存器，TDDR是一个16为寄存器。两个寄存器都用于存放分频比。每次向计时器状态寄存器写入开始计时指令后，TDDR就会被装入新的值，并且将此值装入到减法计数器PSC中。每来一个CPU时钟时（只要计时器没有停止工作），就把PSC自减1.当PSC减到0时，就会想TIM高位产生借位信号，使之自减1，同时重新装载TDDR中的数值。而TIM的初值是由PRD设置的，当TIM被减到0时，就会通过TINT产生定时中断2）AD：XF AD的地址为7FFFFH，也就是只要地址线上出现这个地址，CS就会变成低电平，AD被选中。REFP是AD参考电压，也就是AD的量程。由于RD接XF，而XF由CPU中状态寄存器的一位是相同的，所以要想读使能，就必须用软件将XF置成低电平。这样，再发出读信号，就可以从数据线上读出AD数据。而读完后要软件关闭RD。核心代码叙述：1） -e resets表示从resets标号开始执行程序。中断向量表的设置如下：.sect vectors .ref jump .ref starts .def resetsk_stack_size .set 200k_stack .usect stack_section,k_stack_sizesystem_stack .set k_stack+k_stack_sizeresets: bd starts stm #system_stack,spnmi: rete nop nop nopsint17 .space 4*16sint18 .space 4*16sint19 .space 4*16sint20 .space 4*16sint21 .space 4*16sint22 .space 4*16sint23 .space 4*16sint24 .space 4*16sint25 .space 4*16sint26 .space 4*16sint27 .space 4*16sint28 .space 4*16sint29 .space 4*16sint30 .space 4*16 int0: rete nop nop nopint1: rete nop nop nopint2: rete nop nop nop tint0: b jump nop nop将中断向量将表保存，后缀名为.asm。2）数据空间的分配：buferdatax .usect filterdata,k_bufsize*2 ;循环缓冲区,out_data .usect filterdata,1 ;AD输出信号,即输入计算机的采样信号d_cr0_send .usect filiterdata,1;保存AD初始化状态字1d_cr1_send .usect filterdata,1;保存AD初始化状态字2d_temp .usect filterdata,1 ;用于保存临时读入的数据3）设置中断向量表的首地址。这个地址的高九位放在PMST寄存器中，而低七位由产生外部中断时，随之一起产生的中断类型码乘以4得出。starts: ;设置中断向量表首地址高九位; ldm pmst,a and #k_temp,a or #k_iptr,a stl a,pmst4）打开外部中断。外部中断有两级控制，一个是DSP内核外的中断屏蔽寄存器，他可以对外部可屏蔽中断按位进行屏蔽。另一级是DSP内的状态寄存器的中断屏蔽位。 STM #0008h,IMR nop nop RSBX INTM ;软件将CPU内部外部中断使能位打开 nopSSBX FRCT ;对小数整数位置1 nop nop 5）对AD进行初始化，以及各个地址寄存器值的初始化STM #d_temp,ar2 ;将临时单元地址放在ar2 ;将AD初始化状态字放入相应的内存单元,为AD初始化做准备STM #d_cr0_send,ar1 ST #k_cr0_send,*ar1+ ST #k_cr1_send,*ar1 STM #d_cr0_send,ar1 STM#out_data,ORIGIN ;设置滤波器循环缓冲区 stm #k_cir,bk;循环缓冲区大小 stm #1,ar0 ;循环缓冲区步长 stm #buferdatax,INPUT;将循环缓冲区首地址放在INPUT寄存器内 RPT #k_a-1; STM #0, *INPUT+0% ;输入信号的初始状态全为0 ;对AD初始化 portr 0FFFFh,*ar2 ;使AD的CS不使能 nop nop nop portw *ar1+,07FFh ;写入AD的寄存器CR0 rpt #8 nop portw *ar1,07FFh ;写入AD的寄存器CR1 portr 0FFFFh,*ar2 ;使AD的CS不使能 rpt #10 nop nop在这之后，AD已经开始工作。6）对定时器进行初始化：stm #0010h,tcr ;停止定时器的工作 stm #0683h,prd ;设置PRD分频比. stm #0825h,tcr ;设置TDDR分频比为5,始计时.总分频为10000,原始时钟频率为60MHz,所以分得6000Hz nop nop TCR的最低十六进制位和PRD整个寄存器的乘积就是分频比。 最后一条指令将两个分频比的值装入到分频寄存器中，并且开始计时。 上述步骤之后，就要进循环，等待定时中断的到来。每进行一个计时周期，就会产生一个定时中断。 wait: nop nop b wait一旦产生中断，在允许中断的情况下，就会进入中断服务程序。 首先要打开AD的读使能。然后从AD读一个采样数据。之后，为了防止溢出，要先清除掉高六位，用这个数据进行滤波。rsbx xf ;AD读使能,为接下来的读数做准备,同时开始下一次转换rpt #5 nop portr 07FFFh,*ORIGIN rpt #5 nop LD *ORIGIN,16,A SFTA A,6 STH A,*ORIGIN 7）滤波MVDD *ORIGIN, *INPUTRPT #K_A-2MAR *INPUT-0%MPY *INPUT+0%, #a20,B LD B, AMPY *INPUT+0%, #a19,BADD B, AMPY *INPUT+0%, #a18,BADD B, AMPY *INPUT+0%, #a17,BADD B, AMPY *INPUT+0%, #a16,BADD B, AMPY *INPUT+0%, #a15,BADD B, AMPY *INPUT+0%, #a14,BADD B, AMPY *INPUT+0%, #a13,BADD B, AMPY *INPUT+0%, #a12,BADD B, AMPY *INPUT+0%, #a11,BADD B, AMPY *INPUT+0%, #a10,BADD B, AMPY *INPUT+0%, #a9,BADD B, AMPY *INPUT+0%, #a8,BADD B, AMPY *INPUT+0%, #a7,BADD B, AMPY *INPUT+0%, #a6,BADD B, AMPY *INPUT+0%, #a5,BADD B, AMPY *INPUT+0%, #a4,BADD B, AMPY *INPUT+0%, #a3,BADD B, AMPY *INPUT+0%, #a2,BADD B, AMPY *INPUT+0%, #a1,BADD B, AMPY *INPUT+0%, #a0,BADD B, ASFTA A,-4 ;sfta算术左移STH A, *OUTPUT ;本段只能一次完成一个数的滤波ssbx xf ;关闭读使能nopportr 0FFFFh,*ar2 ;AD片选不使能nopnopportw *OUTPUT,0BFFFh ;将滤波数据通过DA进行输出nopnoprete.end由于滤波断已经在之前已经做过详细的解释了，这里不再赘述，不同的是之前的滤波程序先放了20个x(n)，而在这里先放了20个零，并且在滤波后得到的数据进行了简单的处理，右移四位，防止D/A输出错误。8）滤波之后，要关闭读使能，并且关闭AD片选使能。之后，将滤波得到的数据在DA输出。然后中断返回，完成服务程序。 ssbx xf ;关闭读使能 nop portr 0FFFFh,*ar2 ;AD片选不使能 nop nop portw *OUTPUT,0BFFFh ;将滤波数据通过DA进行输出 nop nop Rete六、结果展示:打开工程，载入工程，运行程序，将信号发生器接在A/D输入上，将示波器的CH1和CH2分别接在A/D输入和D/A输出上。将信号的频率设置在1kHz，观察示波器。观察到示波器的CH1和CH2上展现的一个完整的正弦波形。表明1KHz时能够通过。调节频率扫描旋钮，减小频率，可以看到在减小频率到某个值得时候就会出现混乱，没有负的波形，反而出现一个很大的正值，这就是前面提到的溢出，此时调节一个合适的直流电平就可以解决这个问题。调节直流电平后，再次调节频率至1kHz，可以发现这时波形正常，继续增大频率，可以看到CH2的峰峰值慢慢的减小，至1.25kHz时，波形减小到大约为CH1的0.75倍。继续增大频率，可以看到CH2的峰峰值继续在减小，到1.5KHz时，峰峰值已经到了CH1的不到0.15倍。继续增大频率，可以看到CH1基本成为一条直线。但是随着频率的继续增大，每隔一段频率可以看到，CH1的上面依然会有非常小的波形，只是这个值远远不到0.1倍。七、结果分析和改进：可以清楚的看到，滤波器确实完成了低通的特性要求，保证了1kHz信号的通过，1