如何用matlab进行滤波器设计和仿真.docx
《如何用matlab进行滤波器设计和仿真.docx》由会员分享,可在线阅读,更多相关《如何用matlab进行滤波器设计和仿真.docx(21页珍藏版)》请在冰豆网上搜索。
![如何用matlab进行滤波器设计和仿真.docx](https://file1.bdocx.com/fileroot1/2023-2/24/62b18064-ef12-49d5-98a1-da7cfe19e6a4/62b18064-ef12-49d5-98a1-da7cfe19e6a41.gif)
如何用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分):
对本课程的感想。
总分:
评卷人: