数字信号处理实验-用MATLAB设计循环卷积系统仿真

用MATLAB 设计循环卷积系统仿真

专 业： 电子信息工程 学 号： 姓 名：

2014年12月

一、实验目的

（1） 熟悉使用MATLAB 软件。

（2） 学会调用MATLAB 信号处理工具的设计函数。 （3） 对循环卷积有更深的认识和理解。

二、实验原理和步骤

1、卷积的定义： 任意信号

都可以根据不同需要进行不同的分解。如信号

可以分解为直

流分量和交流分量，也可以分解为奇分量和偶分量，或分解为实部分量和虚部分

量。如果信号费解为冲击信号，那么信号分解为一系列不同强度，不同时延的冲击信号的叠加，这个过程称为卷积积分。

一般而言，如果有两个函数f 1(t )和f 2(t )，则它们的积分

y (t )=⎰f 1(τ)f 2(t -τ)d τ称为f 1(t )与f 2(t )的卷积积分，简称卷积，表达式为：

-∞+∞

y (t )=f 1(t )*f 2(t )，即:

y (t )=f 1(t )*f 2(t )=⎰f 1(τ)f 2(t -τ)d τ

-∞

+∞

2、线性卷积的运算：

卷积运算是线性时不变系统分析的重要工具，很多滤波器的设计中都要用到卷积运算。给出线性卷积运算的定义，设有离散信号x(n)和y(n)，其线性卷积为：

线性卷积有四步运算：①卷积运算时，y(n)要先反折得到y(-n);②m>0表示y(-n)序列右移，m

值。线性卷积运算简

式中 “∗”表示线性卷积运算符。

由线性卷积的定义

，等式右边是乘积求和形式，，

因而考虑能否用矩阵相乘的形式来表示线性卷积。假设序列x(n) 长度为4点，y(n)

长度为3点，x(n) 除区间之外皆为零，y(n) 除区间之外皆为零，用矩阵的形式来表达线性卷积Z:

0 0 0

0 0

Z=

0 0

0 0 0

x(n),y(n)序列长度不同，则将短序列补0使两者相同。

3、循环卷积的运算

有限长序列的循环移位是指y((m-n))

，也就是先让序列y(n)以N 为周

期进行周期延拓，再进行反折，然后朝右移位，只朝一个方向移位的原因是：对周期序列向右移动一个位置，也就相当于向左移动了N －1个位置, 最后取（0，N －1）的N 个值就得到了循环移位后的N 个序列值。

设有序列x(n)和y(n)，其N 点循环卷积为：

由于循环移位的关系最后得到的循环卷积的长度就是N 点，m 取[0，1，2，„，N-1]。

循环卷积的简洁表示为：

式中表示循环卷积运算符。

例如N=4的循环卷积如下：

其中，N ≥length(y(n))。

值得说明的是，当N ≥length(y(n))+ length(x(n))-1时，圆周卷积的值等于线性卷积。

三、MATLAB 设计循环卷积

1、循环卷积的分析：

两个序列的循环卷积可以分三个步骤完成：

(1)初始化：确定循环点数N ，测量输入2个序列的长度，长度小于N 的在后面补0。

(2)循环右移函数：将序列x(n)循环右移，一共移N 次（N 为循环卷积的循环次数），最后将每次循环成的新序列组成一个矩阵V 。

(3)相乘：将x(n)移位后组成的矩阵V 与第二个序列h(n)对应相乘，即得循环卷积结果。

2、根据循环卷积分析设计流程图：

循环卷积流程图如图1所示: （1）、主流程图

(2)

图4.1循环卷积流程图

（3）、据循环卷积流程图设计matlab 源代码：

function y=myconv(x1,x2)

x1=input('x1='); x2=input('x2='); N=input('N=');

x1=[x1,zeros(1,N-length(x1))]; x2=[x2,zeros(1,N-length(x2))]; V=circlel(x2) Z=x1*V;

stem(Z');xlabel('n');ylabel('Z');grid on; title('循环卷积结果Z') 运行程序，输入序列x1，x2 x1=[-1 2 3 -5] x2=[6 7 -10 4 12] 循环卷积结果：

[10 -55 42 -33 -69 86] 运行图形如图2所示

图4.2循环卷积运行结果

四、实验总结

总结本次数字信号处理实训，我受益匪浅。 首先就是方案的确定。由于这个学期我学习了数字信号处理这门课程，课程中我了解到要实现两信号的卷积，可以通过定义来实现，也可以通过DFT 来计算线性卷积。对于有限长序列，存在两种形式的卷积：线性卷积与圆周卷积。由于圆周卷积可以采用DFT 的快速算法——快速傅里叶变换进行运算，运算速度上有很大的优越性。

其中，设计线性卷积有4个步骤，反折、移位 、相乘、求和，而循环卷积则是通过循环移位后得到的矩阵与序列相乘。根据课上老师讲的求法，画出了思路的流程图，然后根据流程图写出程序，事半功倍。

将结果与直接调用matlab 自带的函数比较，结果显示，自己设计的程序是正确的。

通过这次实验，我对卷积和循环的定义、原理、以及实现方法都有了深入的认识。同时也对MATLAB 软件产生了更加浓厚的兴趣。

本次实训，检验了自己的能力，加强了逻辑思维的能力，不过我也发现了自身存在的一些问题，比如在MATLAB 软件的应用上还有一些功能不懂如何运用的地方，但是在老师和同学的帮助下，我认真学习，并且懂得了许多以前不懂的matlab 的运用。还有很多matlab 的强大功能，希望能在日后好好学习，取得更好的成绩，也希望日后老师能不厌其烦的指导我，给予我更大的支持。

五、参考文献

[1]丁玉美，高西全。《数字信号处理》，第三版，西安电子科技大学出版社。

[2]楼顺天，李博函. 基于MATLAB 的系统分析与设计一信号处理[M].西安:西安电子科技大学出版社，1998.81--88.

用FFT 对信号作频谱分析

专 业： 电子信息工程 学 号： 1208040230 姓 名： 许 先 举

2014年12月

一、 实验目的

学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便正确用FFT 。

二、 实验原理

用FFT 对信号作频谱的分析是学习信号进行谱分析的重要问题是频谱分辨率D 和分析误差。频谱分辨率直接和FFT 的变化区间N 有关，因为FFT 能够实现的频率分辨率是2π/N≤D 。可以根据此式选择FFT 的变换区间N 。误差主要来自于用时，离散谱的包络才能逼近与连续谱，因此N 要适当选择大一些。

周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT ，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量的选择信号的观察时间长一点。

对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的时间长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

三、 实验步骤及内容

（1）、对以下序列进行谱分析：

x 1(n)=R4(n)

x2(n)=

x 3(n)=

n+1 0≤n ≤3 8-n 4≤n ≤7

0 其它n 4-n 0≤n ≤3

n-3 4≤n ≤7

0 其它

n

选择FFT 的变换区间N 为8和16两种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论、分析与比较。 （2）、对以下周期序列进行谱分析：

x 4(n) = cos(π/4)*n ； x5(n)= cos[(π/4)*n]+ cos[(π/8)*n] 选择FFT 的变换区间N 为8和16两种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论、分析与比较。

（3）、对模拟周期信号进行频谱分析：

x 6(n)= cos(8πt)+ cos(16πt)+ cos(20πt)

选择采样频率Fs=64Hz，FFT 的变换区间N 为16、32、64三种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论和分析。

四、思考题

（1）、对于周期序列，如果周期不知道，如何用FFT 进行谱分析？

答：周期信号的周期预先不知道时, 可先截取M 点进行DFT, 再将截取长度扩大1倍截取, 比较结果, 如果二者的差别满足分析误差要求, 则可以近似表示该信号的频谱, 如果不满足误差要求就继续将截取长度加倍, 重复比较, 直到结果满足要求。 （2）、如何选择FFT 的变换区间（包括周期信号与非周期信号）

答：一、对于非周期信号：有频谱分辨率F ，而频谱分辨率直接和FFT 的变换区间有关，因为FFT 能够实现的频率分辨率是2π/N...因此有最小的N>2π/F。就可以根据此式选择FFT 的变换区间。

二、对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT ，得到的离散谱才能代表周期信号的频谱。

五、实验结果如下：

x 1(n)=R4(n)

x2(n)=

x 3(n)=

n+1 0≤n ≤3 8-n 4≤n ≤7

0 其它

n 4-n 0≤n ≤3

n-3 4≤n ≤7

0 其它

n

FFT 的变换区间N 为8和16两种情况进行频谱分析

实验结果图形与理论分析相符。 （2）对以下周期序列进行谱分析：

x4(n)=cos[(π/4)*n]

x5(n)= cos[(π/4)*n]+ cos[(π/8)*n]

选择FFT 的变换区间N 为8和16两种情况进行频谱分析，分别打印出幅频

特性曲线，并进行讨论、分析与比较。 实验结果如下：

（3）对模拟周期信号进行频谱分析：

x6(n)= cos(8πt)+ cos(16πt)+ cos(20πt)

选择采样频率Fs=64Hz，FFT 的变换区间N 为16、32、64三种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论、分析与比较。

实验结果如下：

六、实验中代码：

x1n=[ones(1,4)]; %产生R4(n)序列向量 X1k8=fft(x1n,8); %计算x1n 的8点DFT X1k16=fft(x1n,16); %计算x1n 的16点DFT %以下绘制幅频特性曲线 N=8;

f=2/N*(0:N-1); figure(1);

subplot(1,2,1);stem(f,abs(X1k8),'.'); %绘制8点DFT 的幅频特性图 title('(1a) 8点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度'); N=16;

f=2/N*(0:N-1);

subplot(1,2,2);stem(f,abs(X1k16),'.'); %绘制8点DFT 的幅频特性图 title('(1a) 16点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度'); %x2n 和 x3n

M=8;xa=1:(M/2); xb=(M/2):-1:1;

x2n=[xa,xb]; %产生长度为8的三角波序列x2(n) x3n=[xb,xa];

X2k8=fft(x2n,8); X2k16=fft(x2n,16); X3k8=fft(x3n,8); X3k16=fft(x3n,16); figure(2); N=8;

f=2/N*(0:N-1);

subplot(2,2,1);stem(f,abs(X2k8),'.'); %绘制8点DFT 的幅频特性图 title('(2a) 8点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度'); subplot(2,2,3);stem(f,abs(X3k8),'.'); %绘制8点DFT 的幅频特性图 title('(3a) 8点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度'); N=16;

f=2/N*(0:N-1);

subplot(2,2,2);stem(f,abs(X2k16),'.'); %绘制8点DFT 的幅频特性图 title('(2a) 16点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度'); subplot(2,2,4);stem(f,abs(X3k16),'.'); %绘制8点DFT 的幅频特性图 title('(3a) 16点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度'); %x4n 和 x5n N=8;n=0:N-1;

x4n=cos(pi*n/4);

x5n=cos(pi*n/4)+cos(pi*n/8); X4k8=fft(x4n,8); X4k16=fft(x4n,16); X5k8=fft(x5n,8); X5k16=fft(x5n,16); figure(3); N=8;

f=2/N*(0:N-1);

subplot(2,2,1);stem(f,abs(X4k8),'.'); %绘制8点DFT 的幅频特性图 title('(4a) 8点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度'); subplot(2,2,3);stem(f,abs(X5k8),'.'); %绘制8点DFT 的幅频特性图 title('(5a) 8点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度'); N=16;

f=2/N*(0:N-1);

subplot(2,2,2);stem(f,abs(X4k16),'.'); %绘制8点DFT 的幅频特性图 title('(4a) 16点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度'); subplot(2,2,4);stem(f,abs(X5k16),'.'); %绘制8点DFT 的幅频特性图 title('(5a) 16点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度'); %x8n

Fs=64; T=1/Fs;

N=16;n=0:N-1; %对于N=16的情况 nT = n*T;

x8n=cos(8*pi*nT)+cos(16*pi*nT)+cos(20*pi*nT) X8k16=fft(x8n,16); N=16;

f=2/N*(0:N-1); figure(4);

subplot(2,2,1);stem(f,abs(X8k16),'.'); %绘制8点DFT 的幅频特性图 title('(8a) 16点DFT[x_8(n)]');xlabel('ω/π');ylabel('幅度'); N=32;n=0:N-1; %对于N=16的情况 nT = n*T;

x8n=cos(8*pi*nT)+cos(16*pi*nT)+cos(20*pi*nT) X8k32=fft(x8n,32); N=32;

f=2/N*(0:N-1);

subplot(2,2,2);stem(f,abs(X8k32),'.'); %绘制8点DFT 的幅频特性图 title('(8a) 32点DFT[x_8(n)]');xlabel('ω/π');ylabel('幅度'); N=64;n=0:N-1; %对于N=16的情况

nT = n*T;

x8n=cos(8*pi*nT)+cos(16*pi*nT)+cos(20*pi*nT) X8k64=fft(x8n,64); N=64;

f=2/N*(0:N-1);

subplot(2,2,3);stem(f,abs(X8k64),'.'); %绘制8点DFT 的幅频特性图 title('(8a) 64点DFT[x_8(n)]');xlabel('ω/π');ylabel('幅度');

七、实验体会

通过实验，我知道了用FFT 对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D 和分析误差。频谱分辨率直接和FFT 的变换区间N 有关，因为FFT 能够实现的频率分辨率是2л／N ≤D 。可以根据此式选择FFT 的变换区间N 。误差主要来自于用FFT 作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N 较大时，离散谱的包络才能逼近于连续谱，因此N 要适当选择大一些。

周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT ，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。

对模拟信号进行频谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的普分析进行。

八、参考文献

丁玉美，高西全。《数字信号处理》，第三版，西安电子科技大学出版社。

IIR 数字滤波器设计及软件实现

专 业： 电子信息工程 学 号： 1208040230 姓 名： 许 先 举

2014年12月

一、实验目的

（1）熟悉用双线性变换法设计IIR 数字滤波器的原理与方法；

（2）学会调用MATLAB 信号处理工具箱中滤波器设计函数（或滤波器设计分析工具fdatool ）设计各种IIR 数字滤波器，学会根据滤波需求确定滤波器指标参数。

（3）掌握IIR 数字滤波器的MATLAB 实现方法。

（3）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。

二、实验原理

设计IIR 数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标； ②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。MATLAB 信号处理工具箱中的各种IIR 数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter 、cheby1 、cheby2 和ellip 可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函数直接设计IIR 数字滤波器。

本实验的数字滤波器的MATLAB 实现是指调用MATLAB 信号处理工具箱函数filter 对给定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n）。

三、 实验内容及步骤

（1）调用信号产生函数mstg 产生由三路抑制载波调幅信号相加构成的复合信号st ，该函数还会自动绘图显示st 的时域波形和幅频特性曲线，如图10.4.1所示。由图可见，三路信号时域混叠无法在时域分离。但频域是分离的，所以可以通过滤波的方法在频域分离，这就是本实验的目的。

图1 三路调幅信号st 的时域波形和幅频特性曲线

（2）要求将st 中三路调幅信号分离，通过观察st 的幅频特性曲线，分别确定可以分离st 中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0.1dB, 阻带最小衰减为60dB 。

提示：抑制载波单频调幅信号的数学表示式为

1

s (t ) =cos(2πf 0t )cos(2πf c t ) =[cos(2π(f c -f 0) t ) +cos(2π(f c +f 0) t )]

2其中，cos(2πf c t ) 称为载波，f c 为载波频率，cos(2πf 0t ) 称为单频调制信号，f 0为调制正弦波信号频率，且满足f c >f 0。由上式可见，所谓抑制载波单频调幅信号，就是2个正弦信号相乘，它有2个频率成分：和频f c +f 0和差频f c -f 0，这2个频率成分关于载波频率f c 对称。所以，1路抑制载波单频调幅信号的频谱图是关于载波频率f c 对称的2根谱线，其中没有载频成分，故取名为抑制载波单频调幅信号。容易看出，图10.4.1中三路调幅信号的载波频率分别为250Hz 、500Hz 、1000Hz 。如果调制信号m(t)具有带限连续频谱，无直流成分，则

s (t ) =m (t ) cos(2πc f t 就是一般的抑制载波调幅信号。其频谱图是关于载波频率) f c 对称的2个边带（上下边带），在专业课通信原理中称为双边带抑制载波 (DSB-SC) 调幅信号, 简称双边带 (DSB) 信号。如果调制信号m(t)有直流成分，则s (t ) =m (t )cos(2πf c t ) 就是一般的双边带调幅信号。其频谱图是关于载波频率f c 对称的2个边带（上下边带），并包含载频成分。

（3）编程序调用MATLAB 滤波器设计函数ellipord 和ellip 分别设计这三个椭圆滤波器，并绘图显示其幅频响应特性曲线。

（4）调用滤波器实现函数filter ，用三个滤波器分别对信号产生函数mstg 产生的信号st 进行滤波，分离出st 中的三路不同载波频率的调幅信号y 1(n)、y 2(n)和y 3(n)， 并绘图显示y1(n)、y2(n)和y3(n)的时域波形，观察分离效果。

四、信号产生函数mstg 清单

function st=mstg

%产生信号序列向量st, 并显示st 的时域波形和频谱

%st=mstg 返回三路调幅信号相加形成的混合信号，长度N=1600 N=1600 %N为信号st 的长度。

Fs=10000;T=1/Fs;Tp=N*T; %采样频率Fs=10kHz，Tp 为采样时间 t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;

fc1=Fs/10; %第1路调幅信号的载波频率fc1=1000Hz, fm1=fc1/10; %第1路调幅信号的调制信号频率fm1=100Hz

fc2=Fs/20; %第2路调幅信号的载波频率fc2=500Hz fm2=fc2/10; %第2路调幅信号的调制信号频率fm2=50Hz fc3=Fs/40; %第3路调幅信号的载波频率fc3=250Hz, fm3=fc3/10; %第3路调幅信号的调制信号频率fm3=25Hz xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号 xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号 xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号 st=xt1+xt2+xt3; %三路调幅信号相加 fxt=fft(st,N); %计算信号st 的频谱

%====以下为绘图部分，绘制st 的时域波形和幅频特性曲线==================== subplot(3,1,1)

plot(t,st);grid;xlabel('t/s');ylabel('s(t)');

axis([0,Tp/8,min(st),max(st)]);title('(a) s(t)的波形' ) subplot(3,1,2)

stem(f,abs(fxt)/max(abs(fxt)),'.' );grid;title('(b) s(t)的频谱' ) axis([0,Fs/5,0,1.2]);

xlabel('f/Hz');ylabel(' 幅度' )

五、实验程序框图如图2所示

图2 实验程序框图

六、滤波器参数及实验程序清单

（1）、滤波器参数选取

观察图1可知，三路调幅信号的载波频率分别为250Hz 、500Hz 、1000Hz 。带宽（也可以由信号产生函数mstg 清单看出）分别为50Hz 、100Hz 、200Hz 。所以，分离混合信号st 中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的指标参数选取如下：

对载波频率为250Hz 的条幅信号，可以用低通滤波器分离，其指标为

带截止频率

f p =280Hz ，通带最大衰减αp =0.1dB dB ； f s =450Hz ，阻带最小衰减αs =60dB dB ，

阻带截止频率

对载波频率为500Hz 的条幅信号，可以用带通滤波器分离，其指标为 带截止频率

f pl =440Hz ，f pu =560Hz ，通带最大衰减αp =0.1dB dB ；

f sl =275

Hz ，

阻带截止频率

f su =900Hz ，Hz ，阻带最小衰减

αs =60dB dB ，

对载波频率为1000Hz 的条幅信号，可以用高通滤波器分离，其指标为 带截止频率

f p =890Hz ，通带最大衰减αp =0.1dB dB ； f s =550Hz ，阻带最小衰减αs =60dB dB ，

阻带截止频率

说明：（1）为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是

尽量使滤波器过渡带宽尽可能宽。

（2）与信号产生函数mstg 相同，采样频率Fs=10kHz。 （3）为了滤波器阶数最低，选用椭圆滤波器。

按照图2 所示的程序框图编写的实验程序为exp4.m 。 （2）、实验程序清单 %实验4程序exp4.m

% IIR数字滤波器设计及软件实现 clear all ;close all

Fs=10000;T=1/Fs; %采样频率

%调用信号产生函数mstg 产生由三路抑制载波调幅信号相加构成的复合信号st st=mstg;

%低通滤波器设计与实现========================================= fp=280;fs=450;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord 计算椭圆DF 阶数N 和通带截止

频率wp

[B,A]=ellip(N,rp,rs,wp); %调用ellip 计算椭圆带通DF 系统函数系数向量B 和A

y1t=filter(B,A,st); %滤波器软件实现

% 低通滤波器设计与实现绘图部分

figure(2);subplot(3,1,1);

myplot(B,A); %调用绘图函数myplot 绘制损耗函数曲线

yt='y_1(t)';

subplot(3,1,2);tplot(y1t,T,yt); %调用绘图函数tplot 绘制滤波器输出波形 %带通滤波器设计与实现

====================================================

fpl=440;fpu=560;fsl=275;fsu=900;

wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;

[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord 计算椭圆DF 阶数N 和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp); %调用ellip 计算椭圆带通DF 系统函数系数向量B 和A y2t=filter(B,A,st); %滤波器软件实现

figure(3);subplot(3,1,1);

myplot(B,A); %调用绘图函数myplot 绘制损耗函数曲线

yt='y_2(t)';

subplot(3,1,2);tplot(y2t,T,yt); %调用绘图函数tplot 绘制滤波器输出波形 %高通滤波器设计与实现================================================ fp=890;fs=600;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord 计算椭圆DF 阶数N 和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp,'high' ); %调用ellip 计算椭圆带通DF 系统函数系数向量B 和A

y3t=filter(B,A,st); %滤波器软件实现

figure(4);subplot(3,1,1);

myplot(B,A); %调用绘图函数myplot 绘制损耗函数曲线

yt='y_3(t)';

subplot(3,1,2);tplot(y3t,T,yt); %调用绘图函数tplot 绘制滤波器输出波形

调用的子函数:

(1)myplot:计算时域离散系统损耗函数并绘制曲线图。函数清单如下： function myplot(B,A)

[H,W]=freqz(B,A,1000);

m=abs(H);

plot(W/pi,20*log10(m/max(m)));grid on ;

xlabel('\omega/\pi');ylabel(' 幅度(dB)')

axis([0,1,-80,5]);title(' 损耗函数曲线' );

(2)tplot：时域序列连续曲线绘图函数，将采样序列绘图。函数清单如下： function tplot(xn,T,yn)

n=0:length(xn)-1;t=n*T;

plot(t,xn);

xlabel('t/s');ylabel(yn)

axis([0,t(end),min(xn),1.2*max(xn)]);

七、实验程序运行结果

实验程序exp4.m 运行结果如图104.2所示。由图可见，三个分离滤波器指标参数选取正确，算耗函数曲线达到所给指标。分离出的三路信号y1(n)，y2(n)和y3(n)的波形是抑制载波的单频调幅波。

(a) 低通滤波器损耗函数及其分离出的调幅信号y 1

(t)

(b) 带通滤波器损耗函数及其分离出的调幅信号y 2(t)

(c)高通滤波器损耗函数及其分离出的调幅信号y 3(t)

图3 实验程序exp4.m 运行结果

八、参考文献

丁玉美，高西全。《数字信号处理》，第三版，西安电子科技大学出版社。

FIR 数字滤波器设计与软件实现

专 业： 电子信息工程

学 号： 1208040230

姓 名： 许 先 举

2014年12月

一、实验目的

（1）掌握用窗函数法设计FIR 数字滤波器的原理和方法。

（2）掌握用等波纹最佳逼近法设计FIR 数字滤波器的原理和方法。

（3）掌握FIR 滤波器的快速卷积实现原理。

（4）学会调用MATLAB 函数设计与实现FIR 滤波器。

二、实验内容及步骤

（1）认真复习第七章中用窗函数法和等波纹最佳逼近法设计FIR 数字滤波器的原理；

（2）调用信号产生函数xtg 产生具有加性噪声的信号xt ，并自动显示xt 及其频谱，如图1所示；

图1 具有加性噪声的信号x(t)及其频谱如图

（3）请设计低通滤波器，从高频噪声中提取xt 中的单频调幅信号，要求信号幅频失真小于0.1dB ，将噪声频谱衰减60dB 。先观察xt 的频谱，确定滤波器指标参数。

（4）根据滤波器指标选择合适的窗函数，计算窗函数的长度N ，调用MATLAB 函数fir1设计一个FIR 低通滤波器。并编写程序，调用MATLAB 快速卷积函数fftfilt 实现对xt 的滤波。绘图显示滤波器的频响特性曲线、滤波器输出信号的幅频特性图和时域波形图。

（4）重复（3），滤波器指标不变，但改用等波纹最佳逼近法，调用MATLAB 函数remezord 和remez 设计FIR 数字滤波器。并比较两种设计方法设计的滤波器阶数。

提示：（1）、MATLAB 函数fir1的功能及其调用格式请查阅教材；

（2）、采样频率Fs=1000Hz，采样周期T=1/Fs；

（3）、根据图1(b)和实验要求，可选择滤波器指标参数：通带截止频率fp=120Hz，阻带截至频率fs=150Hz，换算成数字频率，通带截止频率

ωp =2πf p T=0.24π，通带最大衰为0.1dB ，阻带截至频率ωs =2πf s T=0.3π，阻带最小衰为60dB 。

（4）、实验程序框图如图2所示。

图2 实验程序框图

三、 滤波器参数及实验程序清单

1、滤波器参数选取

根据实验指导的提示③选择滤波器指标参数：

通带截止频率fp=120Hz，阻带截至频率fs=150Hz。代入采样频率Fs=1000Hz，换算成数字频率，通带截止频率ωp =2πf p T=0.24π，通带最大衰为0.1dB ，阻带截至频率ωs =2πf s T=0.3π，阻带最小衰为60dB 。所以选取blackman 窗函数。与信号产生函数xtg 相同，采样频率Fs=1000Hz。

按照图2 所示的程序框图编写的实验程序为exp2.m 。

2、实验程序清单

% FIR数字滤波器设计及软件实现

clear all;close all;

%==调用xtg 产生信号xt, xt长度N=1000,并显示xt 及其频谱,========= N=1000;xt=xtg(N);

fp=120; fs=150;Rp=0.2;As=60;Fs=1000; % 输入给定指标

% (1) 用窗函数法设计滤波器

wc=(fp+fs)/Fs; %理想低通滤波器截止频率(关于pi 归一化）

B=2*pi*(fs-fp)/Fs; %过渡带宽度指标

Nb=ceil(11*pi/B); %blackman窗的长度N

hn=fir1(Nb-1,wc,blackman(Nb));

Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性

ywt=fftfilt(hn,xt,N); %调用函数fftfilt 对xt 滤波

%以下为用窗函数法设计法的绘图部分（滤波器损耗函数，滤波器输出信号波形)

f=[0:1023]*Fs/1024;

figure(2)

subplot(2,1,1)

plot(f,20*log10(Hw/max(Hw)));grid;title('(a) 低通滤波器幅频特性')

axis([0,Fs/2,-120,20]);

xlabel('f/Hz');ylabel('幅度')

t=[0:N-1]/Fs;Tp=N/Fs;

subplot(2,1,2)

plot(t,ywt);grid;

axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_w(t)');

title('(b) 滤除噪声后的信号波形')

% (2) 用等波纹最佳逼近法设计滤波器

fb=[fp,fs];m=[1,0]; % 确定remezord 函数所需参数f,m,dev

dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];

[Ne,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez 函数所需参数

hn=remez(Ne,fo,mo,W); % 调用remez 函数进行设计

Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性

yet=fftfilt(hn,xt,N); % 调用函数fftfilt 对xt 滤波

%以下为用等波纹设计法的绘图部分（滤波器损耗函数，滤波器输出信号波形)

figure(3);subplot(2,1,1)

f=[0:1023]*Fs/1024;

plot(f,20*log10(Hw/max(Hw)));grid;title('(c) 低通滤波器幅频特性')

axis([0,Fs/2,-80,10]);

xlabel('f/Hz');ylabel('幅度')

subplot(2,1,2);plot(t,yet);grid;

axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_e(t)');

title('(d) 滤除噪声后的信号波形')

四、 实验程序运行结果

用窗函数法设计滤波器，滤波器长度 Nb=184。滤波器损耗函数和滤波器输出yw(nT)分别如图3(a)和（b ）所示。

用等波纹最佳逼近法设计滤波器，滤波器长度 Ne=83。滤波器损耗函数和滤波器输出ye(nT)分别如图3(c)和（d ）所示。

两种方法设计的滤波器都能有效地从噪声中提取信号，但等波纹最佳逼近法设计的滤波器阶数低得多，当然滤波实现的运算量以及时延也小得多，从图3(b)和（d ）可以直观地看出时延差别。

图3

五、参考文献

丁玉美，高西全。《数字信号处理》，第三版，西安电子科技大学出版社。