实验二因果性数字系统的时域实现.doc

贵州大学实验报告
学院： 专业： 班级：
姓名

学号

实验组

实验时间
2010、10
指导教师
刘姐
成绩

实验项目名称
实验二：因果性数字系统的时域实现
实验目的
通过本实验的学习使学生熟悉和掌握数字信号处理系统中因果性数字系统的时域实现。

实验要求
必修
实验原理
根据数字滤波器直接形式I的结构，先实现FIR滤波器部分（nonrec.m），再实现纯递归的IIR滤波（rec.m），级联后可实现上面差分方程的运算。
实验仪器
计算机一台；
Matlab软件
实验步骤
1、编写nonrec.m函数文件，实现FIR滤波。这里给定，，求。
2、编写rec.m函数文件，实现纯递归IIR滤波。这里给定，，，，求单位取样响应。
3、用nonrec.m和rec.m函数编制dfilter.m函数文件，构成完整的一般IIR滤波器程序，并完成下列信号滤波

这里给定系统函数
，计算。
4、用help查看内部函数filter.m，了解调用格式，重做3，并和你编写的dfilter.m的结果进行比较。
实验内容
1、编程实现FIR滤波

2、编写时实现纯递归IIR滤波；
差分方程：

系统函数：

3、调用库函数filter.m实现IIR滤波
实验数据
1. 编制nonrec.m函数文件,实现FIR滤波y(n)=h(n)*x(n).这里给定h(n)=R8(n),x(n)=nR16(n),求y(n).
nonrec.m函数文件:
function y=nonrec(x,h)
x=[x,zeros(1,length(h)-1)]; %补零
w=zeros(1,length(h));
for I=1:length(x)
for j=length(h):-1:2
w(j)=w(j-1); %得到每一延时单元输出变量
end
w(1)=x(i);
y(i)=w*h’; %Bcofi与wi对应相乘
end

主程序文件:
x=0:15;
h=ones(1,8);
y=nonrec(x,h);
n=0:22;
stem(n,y);

分析：线性卷积y(n)=x(n)*h(n)的长度为16+8-1=23，可利用y(n)=∑h(m)x(n-m)直接计算得
n(n+1)/2, n≤7
y(n)= 4(2n-7), 8≤n≤15
(n+8)(23-n)/2, 16≤n≤22
即 y=[ 0 1 3 6 10 15 21 28 36 44 52 60 68 76 84 92 84 75
65 54 42 29 15] ,与曲线相符。

2. 编制rec.m函数文件,实现纯递归IIR滤波 y(n)=x(n)+∑aky(n-k). 这里给定a1=2rcosw0, a2=-r2, r=0.95, w0=π/8, 求单位抽样响应h(n).

rec.m函数文件:
function y=rec(x,a,n)
x=[x,zeros(1,n-length(x))]; %补零到所需长度
sum=0;
w=zeros(1,length(a));
for i=1:n
y(i)=sum+x(i); % 递归
for j=length(a):-1:2
w(j)=w(j-1); %延时
end
w(1)=y(i);
sum=w*a';
end

主程序文件:
x=[1];
a=[2*0.95*cos(pi/8),-0.95^2];
h=rec(x,a,75); %取h(n)的长度为75点
n=0:74;
stem(n,h);

分析计算:
由题意, a1=2*0.95*cos(π/8), a2=-0.952, 所以,
得到系统函数 H(z)=1/[1-1.9cos(π/8)z-1+0.952z-2],
做逆Z变换得 h(n)=0.95ncos(πn/8)+ctg(π/8)*0.95nsin(πn/8),
利用MATLAB直接画h(n), 即使用下列语句
n=0:74;
h=0.95.^n.*cos(pi.*n./8)+cot(pi/8).*(0.95.^n).*sin(pi.*n./8);
stem(n,h);

3. 用nonrec.m和rec.m函数编制dfilter.m 函数文件, 构成完整的一般IIR滤波程序, 并完
成下列信号滤波:
x(n)=cos(2πn/5)R64(n)
这里给定系统函数
H(z)=(1-2z-1+z-2)/(1-0.5z-1+0.5z-2)
求y(n).

dfilter.m函数文件:
function y=dfilter(x,b,a,n)
y1=nonrec(x,b);
y=rec(y1,a,n);

主程序文件:
n=0:63;
x=cos(2*pi/5*n);
b=[1,-2,1]; %由H(z)得到系数a,b
a=[0.5,-0.5];
y=dfilter(x,b,a,64); %取y(n)的长度为64点
stem(n,y);

4. 用help查看内部函数filter.m, 了解调用格式, 重做3, 并与你编写的dfilter.m结果进行比较 .

用help可以看到内部函数为Y = FILTER(B,A,X), 且有 The filter is a "Direct Form II Transposed" implementation of the standard difference equation:
a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
-a(2)*y(n-1) - ... - a(na+1)*y(n-na)
因此, 调用内部函数filter时, 要对原系数a做适当变化: a1(1)=1, a1(i)=-a(i), i>1.
n=0:63;
x=cos(2*pi/5*n);
b=[1,-2,1];
a=[0.5,-0.5];
y=dfilter(x,b,a,64); %调用自己编的dfilter函数
a1=[1,-0.5,0.5]; %a变为a1, 用于调用内部函数filter
y1=filter(b,a1,x);
subplot(2,1,1);
stem(n,y);
subplot(2,1,2);
stem(n,y1);

5. 重叠保留法做线性卷积
x=0:134; %x为长序列
y=[];
h=[ones(1,8),zeros(1,8)]; %h为单位抽样响应
H=fft(h);
for i=0:14
if (i==0) xk=[zeros(1,7),0:8]; %x序列的前面补零M-1点(M-1=7)
else xk=x(9*i-6:9*i+9); %截取xk(n)中的N点(N=16)
end
XK=fft(xk);
YK=H.*XK;
yk=ifft(YK); %利用FFT和IFFT做循环卷积
y=[y,yk(8:16)]; %截取yk(n)中的未混叠部分并衔接
end
stem(x,y);
利用MATLAB内部函数filter对此结果进行验证, 输入语句filter(h,[1],x),
则可得到如下曲线. 与重叠保留法所得结果相同.
实验总结
1、内部函数filter.m的调用格式是什么？与编写的dfilter.m调用格式是否一致？差异在何处？
答: 内部函数filter.m的调用格式为Y = FILTER(B,A,X), 与我编写的dfilter.m调用格式略有差异. 我编写的dfilter的格式为dfilter(x,b,a,n), 多一参数n, 含义为输出序列的长度, 可以在调用时选择由x(n)所激励的输出y(n)的长度. 而filter输出的y(n)长度与x(n)一样.
2、在编写dfilter.m文件时，要调用nonrec.m和rec.m文件，调用时先调用哪个比较合理？为什么？
答: 应先调用nonrec.m, 因为nonrec输出的结果要作为rec的一个参数输入, 所以必须先调用nonrec.m.
3、实验中如何合理地选取单位取样响应的点数？实验步骤3中给定的是具有什么性质的滤波器？
答: 选取的点数应能体现出h(n)曲线的变化趋势和特点, 正好反映出h(n)由较大值衰减到零的特性, 因此, 我选择75点, 得到上面二.2中的曲线.
指导教师意见

签名： 年 月 日