如何用matlab进行滤波器设计和仿真.docx
《如何用matlab进行滤波器设计和仿真.docx》由会员分享，可在线阅读，更多相关《如何用matlab进行滤波器设计和仿真.docx（21页珍藏版）》请在冰豆网上搜索。
如何用matlab进行滤波器设计和仿真
华中科技大学
本科课程考试答题本
考生姓名
考生学号
班级
联系电话
Email：
考试科目
考试日期
一、报告摘要：
本报告首先论述如何用matlab进行滤波器设计和仿真，并且将参数转换成DSP可用的滤波器参数。
然后，将这组参数在CCS上进行纯软件仿真。
随之是本次实验过程中代码拼接的问题与解决方法。
软硬件流程图及核心代码的陈述。
将实际滤波的结果进行展示。
与理想滤波器的比较之间的差异与改良。
最后是对于本课程的感想
二、滤波器基础知识以及基于matlab的参数设计：
设计一个能够将1000Hz一下保留，1500Hz以上的频率滤出的滤波器，利用理想的fir滤波器截止频率可以设定为（1000+1500）/2=1250Hz，利用哈明窗设计本滤波器，阶数为21阶，
function[a]=ccc（）
clearall;
clf;
N=1024;%设置采样点数
fs=6000;%设置采样频率
dt=1/fs;
fork=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）;
end
wp=pi*1000*2/6000;%归一化的通带频率
ws=pi*1500*2/6000;%归一化的阻带频率
wc=（wp+ws）/2;%理想的滤波器的截止频率
M=21;%设置阶数
b=fir1（M-1,wc/pi,hamming（M））%;利用fir1函数得到滤波系数b
[H,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）。
滤波系数：
A0
A1
A2
A3
A4
A5
A6
A7
A8
A9
A10
0.0013
-0.0026
-0.0058
0.0032
0.0211
0.0089
-0.0470
-0.0608
0.0726
0.3007
0.4169
A11
A12
A13
A14
A15
A16
A17
A18
A19
A20
0.3007
0.0726
0.0608
-0.0470
0.0089
0.0211
0.0032
-0.0058
-0.0026
0.0013
转换后得到的滤波系数如下：
a0002AHa1-0054Ha2-00BEHa30068H
a402B4Ha50124Ha6-0605Ha7-07C9H
a8094CHa9267DHa10355CHa11267DH;
a12094CHa13-07C9Ha14-0605Ha150124H
a1602B4Ha170068Ha18-00BEHa19-0054H
a20002AH
三、算法验证：
将matlab得到的的滤波器系数转换成十六进制数，并设计一个21阶的fir滤波器在ccs2.0中进行仿真。
STM#K_CIR,BK;设置循环区
STM#1,AR0;设置便宜的步长
STM#inputdata,ORIGIN
STM#bufferdatax,INPUT
STM#filterdata,OUTPUT
RPT#K_A-2
MVDD*ORIGIN+,*INPUT+0%；放前20个数据，并将INPUT指向下一位，同时将ORIGIN指向下一位，
STM#K_DATA_SIZE-K_A-1,BRC
RPTBfilter_end-1
MVDD*ORIGIN+,*INPUT；在第一次滤波时将第21个数据放入，并且在之后的每次滤波中均放进一个数据，用于本次滤波
RPT#K_A-2
MAR*INPUT-0%；将INPUT左移动20个，第一次滤波时正好将INPUT指向x
（1），就能保证滤波的进行，
MPY*INPUT+0%,#a20,B
LDB,A；将该数存进累加器A，
MPY*INPUT+0%,#a19,B
ADDB,A
………该段为a18至a1与对应的值相乘，与a19的操作完全相同，为简化文档，省去。
MPY*INPUT+0%,#a0,B
ADDB,A
STHA,*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）程序有浮点数运算，必须加上SSBXFRCT。
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）-eresets表示从resets标号开始执行程序。
中断向量表的设置如下：
.sect"vectors"
.refjump
.refstarts
.defresets
k_stack_size.set200
k_stack.usect"stack_section",k_stack_size
system_stack.setk_stack+k_stack_size
resets:
bdstarts
stm#system_stack,sp
nmi:
rete
nop
nop
nop
sint17.space4*16
sint18.space4*16
sint19.space4*16
sint20.space4*16
sint21.space4*16
sint22.space4*16
sint23.space4*16
sint24.space4*16
sint25.space4*16
sint26.space4*16
sint27.space4*16
sint28.space4*16
sint29.space4*16
sint30.space4*16
int0:
rete
nop
nop
nop
int1:
rete
nop
nop
nop
int2:
rete
nop
nop
nop
tint0:
bjump
nop
nop
将中断向量
将表保存，后缀名为.asm。
2）数据空间的分配：
buferdatax.usect"filterdata",k_bufsize*2;循环缓冲区,
out_data.usect"filterdata",1;AD输出信号,即输入计算机的采样信号
d_cr0_send.usect"filiterdata",1;保存AD初始化状态字1
d_cr1_send.usect"filterdata",1;保存AD初始化状态字2
d_temp.usect"filterdata",1;用于保存临时读入的数据
3）设置中断向量表的首地址。
这个地址的高九位放在PMST寄存器中，而低七位由产生外部中断时，随之一起产生的中断类型码乘以4得出。
starts:
;;;设置中断向量表首地址高九位;;;;
ldmpmst,a
and#k_temp,a
or#k_iptr,a
stla,pmst
4）打开外部中断。
外部中断有两级控制，一个是DSP内核外的中断屏蔽寄存器，他可以对外部可屏蔽中断按位进行屏蔽。
另一级是DSP内的状态寄存器的中断屏蔽位。
STM#0008h,IMR
nop
nop
RSBXINTM;软件将CPU内部外部中断使能位打开
nop
SSBXFRCT;对小数整数位置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初始化
portr0FFFFh,*ar2;使AD的CS不使能
nop
nop
nop
portw*ar1+,07FFh;写入AD的寄存器CR0
rpt#8
nop
portw*ar1,07FFh;写入AD的寄存器CR1
portr0FFFFh,*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
bwait
一旦产生中断，在允许中断的情况下，就会进入中断服务程序。
首先要打开AD的读使能。
然后从AD读一个采样数据。
之后，为了防止溢出，要先清除掉高六位，用这个数据进行滤波。
rsbxxf;AD读使能,为接下来的读数做准备,同时开始下一次转换
rpt#5
nop
portr07FFFh,*ORIGIN
rpt#5
nop
LD*ORIGIN,16,A
SFTAA,6
STHA,*ORIGIN
7）滤波
MVDD*ORIGIN,*INPUT
RPT#K_A-2
MAR*INPUT-0%
MPY*INPUT+0%,#a20,B
LDB,A
MPY*INPUT+0%,#a19,B
ADDB,A
MPY*INPUT+0%,#a18,B
ADDB,A
MPY*INPUT+0%,#a17,B
ADDB,A
MPY*INPUT+0%,#a16,B
ADDB,A
MPY*INPUT+0%,#a15,B
ADDB,A
MPY*INPUT+0%,#a14,B
ADDB,A
MPY*INPUT+0%,#a13,B
ADDB,A
MPY*INPUT+0%,#a12,B
ADDB,A
MPY*INPUT+0%,#a11,B
ADDB,A
MPY*INPUT+0%,#a10,B
ADDB,A
MPY*INPUT+0%,#a9,B
ADDB,A
MPY*INPUT+0%,#a8,B
ADDB,A
MPY*INPUT+0%,#a7,B
ADDB,A
MPY*INPUT+0%,#a6,B
ADDB,A
MPY*INPUT+0%,#a5,B
ADDB,A
MPY*INPUT+0%,#a4,B
ADDB,A
MPY*INPUT+0%,#a3,B
ADDB,A
MPY*INPUT+0%,#a2,B
ADDB,A
MPY*INPUT+0%,#a1,B
ADDB,A
MPY*INPUT+0%,#a0,B
ADDB,A
SFTAA,-4;sfta算术左移
STHA,*OUTPUT;本段只能一次完成一个数的滤波
ssbxxf;关闭读使能
nop
portr0FFFFh,*ar2;AD片选不使能
nop
nop
portw*OUTPUT,0BFFFh;将滤波数据通过DA进行输出
nop
nop
rete
.end
由于滤波断已经在之前已经做过详细的解释了，这里不再赘述，不同的是之前的滤波程序先放了20个x（n），而在这里先放了20个零，并且在滤波后得到的数据进行了简单的处理，右移四位，防止D/A输出错误。
8）滤波之后，要关闭读使能，并且关闭AD片选使能。
之后，将滤波得到的数
据在DA输出。
然后中断返回，完成服务程序。
ssbxxf;关闭读使能
nop
portr0FFFFh,*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.5kHz信号不通过，截止频率在1.25KHz，但是由于matlab的数据和ccs的数据的差别，滤波特性上略有差别，在1.5kHz的衰减并不是特别理想，理论上这里会衰减到-20dB，也就是0.1倍。
关于后面的出现的随频率间歇的波形，可以从matlab仿真得到的幅频特性曲线看出来，幅值并不是随频率的递增而递减，而是会出现上升然后在下降，但总体是呈下降的趋势。
可以说明，设计的参数基本上能够实现低通的滤波的要求。
改进设想：
根据窗函数法的原理，可以通过增加滤波的阶数，使滤波的特性更加接近理想化。
八、关于本课程的感想
1）通过本实验了解DSP芯片的基本结构，与其他处理器的差别，使用专用的硬件乘法器，采用哈佛结构，指令系统也相对来说比较特别，相比8086系列，有较大的差别。
2）通过本实验，学习了ccs软件，能够利用ccs进行简单的DSP处理器仿真，能后完成简单的功能。
3）关于软件调试，实验中难免遇到错误，而有些错误只是设计时候小问题，例如循环次数等等，这些问题会造成较大的错误，通过调试，可以通过内存分配，地址存放的数值进行判断，相比盲目的寻找错误，这样更加高效可靠。
当然调试过程也不能心浮气躁，否则给自己也会带来很多的麻烦。
4）实验过程中，要了解整个程序的流程，任何一个疏忽都可能导致难以想象的错误，也可能使原本简单的过程变得十分复杂。
评分
题号
得分
题号
得分
第1部分（10分）：
报告摘要，论述报告的组成等信息。
第2部分（15分）：
滤波器系数的获得，包括系数结果、获得系数的方法和工具、以及对系数的高级程序验证方法。
第3部分（15分）：
滤波器系数的DSP纯软件算法验证结果，注意CCS结果展示以及核心代码论述。
第4部分（15分）：
实时滤波器系统代码合成过程的实际记录，论述你在合成代码过程中遇到的所有问题和解决方法。
第5部分（15分）：
实时滤波器系统的软件流程和硬件组成图，核心代码陈述。
第6部分（10分）：
实时滤波器系统的结果展示。
第7部分（10分）：
实时滤波器结果与其理论滤波结果不同的分析，以及改进设想。
第8部分（10分）：
对本课程的感想。
总分:
评卷人:
升级会员 