数字信号处理MATLAB仿真.doc

数字信号处理MATLAB仿真

———————————————————————————————— 作者：
———————————————————————————————— 日期：

13

个人收集整理 勿做商业用途
实验一 数字信号处理的Matlab仿真
一、实验目的
1、掌握连续信号及其MATLAB实现方法；
2、掌握离散信号及其MATLAB实现方法
3、掌握离散信号的基本运算方法，以及MATLAB实现
4、了解离散傅里叶变换的MATLAB实现
5、了解IIR数字滤波器设计
6、了解FIR数字滤波器设计1
二、实验设备
计算机，Matlab软件
三、实验内容
(一）、 连续信号及其MATLAB实现
1、 单位冲击信号

例1.1:t=1/A=50时，单位脉冲序列的MATLAB实现程序如下：
clear all;
t1=—0.5:0。001：0；
A=50；
A1=1/A；
n1=length(t1）；
u1=zeros(1，n1);
t2=0：0。001：A1；
t0=0；
u2=A*stepfun（t2，t0);
t3=A1：0。001：1；
n3=length（t3)；
u3=zeros(1,n3)；
t=[t1 t2 t3］;
u=[u1 u2 u3］；
plot(t，u)
axis(［—0。5 1 0 A+2］）

2、 任意函数

例1.2：用MATLAB画出如下表达式的脉冲序列

clear all;
t=—2：1:3；
N=length（t)；
x=zeros（1，N);
x（1）=0。4；
x(2)=0。8
x（3）=1。2;
x(4）=1.5；
x(5）=1.0；
x(6)=0.7；
stem（t,x);
axis(［-2。2 3。2 0 1.7]）

3、 单位阶跃函数

例1。3：用MATLAB实现单位阶跃函数
clear all;
t=—0。5：0.001：1;
t0=0；
u=stepfun（t,t0）；
plot(t，u）
axis(［-0.5 1 —0.2 1。2］）

4、 斜坡函数

例1.4:用MATLAB实现g(t）=3（t—1)
clear all;
t=0:0。01:3;
B=3；
t0=1；
u=stepfun(t，t0）;
n=length（t）;
for i=1:n
u（i）=B*u（i)＊(t(i）-t0）;
end
plot（t,u)
axis（[—0。2 3。1 —0。2 6.2］）

5、 实指数函数

例1.5:用MATLAB实现
clear all;
t=0:0.001:3；
A=3;
a=0.5；
u=A*exp（a＊t);
plot(t，u)
axis(［-0。2 3.1 -0。2 14］)

6、 正弦函数

例1。6:用MATLAB实现正弦函数f（t)=3cos（10πt+1）
clear all；
t=-0。5：0.001：1；
A=3；
f=5;
fai=1；
u=A＊sin（2*pi*f＊t+fai）;
plot（t,u）
axis(［—0.5 1 -3。2 3.2]）

（二）、离散信号及其MATLAB实现
1、 单位冲激序列

例2.1：用MATLAB产生64点的单位冲激序列
clear all;
N=64;
x=zeros（1,N)；
x(1)=1;
xn=0：N—1；
stem(xn,x)
axis([—1 65 0 1。1])

2、 任意序列

例2.2：用MATLAB画出如下表达式的脉冲序列

clear all;
N=8；
x=zeros(1,N）；
x（1）=8。0;
x（2）=3.4
x（3）=1。8;
x(4)=5.6;
x(5)=2。9;
x（6）=0.7;
xn=0：N—1;
stem（xn，x)
axis(［—1 8 0 8.2])

3、 单位阶跃序列

例2.3:用MATLAB实现单位阶跃函数
clear all；
N=32；
x=ones(1，N);
xn=0:N-1；
stem（xn，x)
axis（[—1 32 0 1.1］)

4、 斜坡序列

例2。4：用MATLAB实现g（n)=3（n-4)点数为32的斜坡序列
clear all;
N=32；
k=4
B=3；
t0=1；
x=［zeros(1,k) ones(1,N-k)]；
for i=1：N
x（i)=B*x（i）＊（i-k）；
end
xn=0：N-1;
stem（xn，x)
axis([—1 32 0 90]）

5、 正弦序列

例2。5:用MATLAB实现幅度A=3，频率f=100，初始相位Φ=1.2,点数为32的正弦信号
clear all；
N=32；
A=3;
f=100；
fai=1.2;
xn=0：N-1;
x=A＊sin（2＊pi＊f*(xn/N)+fai）；
stem(xn,x）
axis([-1 32 —3.2 3.2]）

6、 实指数序列

例2。6：用MATLAB实现,点数为32的实指数序列
clear all；
N=32;
A=3；
a=0。7；
xn=0:N—1；
x=A＊a.^xn；
stem（xn,x）

7、 复指数序列

例2.7:用MATLAB实现幅度A=3,a=0。7，角频率ω=314，点数为32的实指数序列
clear all；
N=32；
A=3；
a=0。7;
w=314；
xn=0：N—1；
x=A*exp((a+j*w)*xn);
stem(xn，x）

8、 随机序列
利用MATLAB产生两种随机信号：
rand(1,N）在区间上产生N点均匀分布的随机序列
randn（1，N）产生均值为0，方差为1的高斯随机序列，即白噪声序列
例2。8:用MATLAB产生点数为32的均匀分布的随机序列与高斯随机序列
clear all;
N=32;
x_rand=rand(1,N）；
x_randn=randn(1,N）;
xn=0：N—1；
figure（1)；
stem（xn,x_rand)
figure(2）;
stem(xn，x_randn）

（三）、离散信号的基本运算
1、 信号的延迟
给定离散信号x(n），若信号y（n）定义为：y(n）=x(n-k)，那么y（n)是信号x（n)在时间轴上右移k个抽样周期得到的新序列。
例3.1：正弦序列y(n)=sin(100n)右移3个抽样周期后所得的序列，MATLAB程序如下：
clear all;
N=32;
w=100;
k=3；
x1=zeros(1,k);
xn=0：N—1;
x2=sin(100*xn);
figure（1)
stem（xn，x2)
x=［x1 x2］；
axis([—1 N -1。1 1.1])
N=N+k;
xn=0：N-1;
figure（2）
stem（xn,x)
axis（［-1 N -1.1 1.1］)

2、 信号相加
若信号，值得注意的是当序列和的长度不相等或者位置不对应时，首先应该使两者的位置对齐，然后通过zeros函数左右补零使其长度相等后再相加
例3。2:用MATLAB实现两序列相加
clear all；
n1=0:3
x1=［2 0.5 0。9 1]；
figure（1)
stem(n1,x1）
axis(［-1 8 0 2。1］ ）
n2=0：7
x2=[ 0 0。1 0。2 0。3 0.4 0.5 0。6 0。7］；
figure(2）
stem(n2，x2）
axis（［-1 8 0 0.8] ）
n=0:7;
x1=［x1 zeros（1,8-length(n1）)］；
x2=［ zeros（1,8—length（n2））,x2];
x=x1+x2;
figure(3)
stem（n，x）
axis(［—1 8 0 2.1]）

3、 信号相乘
信号序列和相乘所得信号的表达式为：
这是样本与样本之间的点乘运算，在MATLAB中可采用“。＊"来实现,但是在信号序列相乘之前,应对其做与相加运算一样的操作。
例3.3：用MATLAB实现上例中两序列相乘
clear all；
n1=0:3
x1=[2 0.5 0.9 1］；
figure（1）
stem(n1,x1)
axis（［-1 8 0 2。1］ ）
n2=0：7
x2=[ 0 0。1 0.2 0.3 0.4 0。5 0.6 0.7]；
figure(2）
stem(n2,x2）
axis([-1 8 0 0。8］ ）
n=0：7;
x1=[x1 zeros（1,8-length(n1))］；
x2=[ zeros(1,8-length(n2)），x2];
x=x1.＊x2;
figure（3)
stem(n，x）
axis（［—1 8 0 0.35])

4、 信号翻转
信号翻转的表达式为：y（n)=x（—n）,在MATLAB中可以用fliplr函数实现此操作
例3.4：用MATLAB实现“信号相加”中的序列翻转
clear all;
n=0:3
x1=［2 0。5 0。9 1]；
x=fliplr(x1);
stem（n，x）
axis([—1 4 0 2。1] ）

5、 信号和
对于N点信号，其和的定义为：
例3.5:用MATLAB实现“信号相加”中的序列和
clear all；
n=0:3
x1=[2 0.5 0。9 1]；
x=sum（x1)

6、 信号积
对于N点信号,其积的定义为：
例3。5：用MATLAB实现“信号相加”中的序列积
clear all；
n=0：3
x1=[2 0。5 0。9 1];
x=prod（x1)

（四）、离散傅里叶变换的MATLAB实现
例４:若是一个N=32的有限序列,利用MATLAB计算它的DFT并画出图形。
N=32;
n=0：N—1;
xn=cos（pi*n/6）;
k=0：N-1；
WN=exp（—j＊2*pi/N)；
nk=n'*k;
WNnk=WN.^nk；
Xk=xn＊WNnk;
figure（1)
stem（n，xn）
figure(2)
stem（k,abs（Xk）)

在MATLAB中，可以直接利用内部函数fft来实现FFT算法，该函数是机器语言，而不是MATLAB指令写成的，执行速度很快。常用格式为：
y=fft（x）
y=fft（x,N)

（五）、IIR数字滤波器设计
1、 基于巴特沃斯法直接设计IIR数字滤波器

例５。1：设计一个10阶的带通巴特沃斯数字滤波器，带通频率为100Hz到200Hz，采样频率为1000Hz，绘出该滤波器的幅频于相频特性，以及其冲击响应图
clear all；
N=10;
Wn=［100 200]/500;
［b,a]=butter（N，Wn，'bandpass'）；
freqz(b,a，128，1000）
figure（2)
[y，t]=impz(b,a，101)；
stem(t，y)

2、 基于切比雪夫法直接设计IIR数字滤波器
例5。2：设计一个切比雪夫Ⅰ型数字低通滤波器，要求：
Ws=200Hz,Wp=100Hz,Rp=3dB,Rs=30dB，Fs=1000Hz

clear all；
Wp=100；
Rp=3；
Ws=200;
Rs=30；
Fs=1000；
［N,Wn]=cheb1ord（Wp/(Fs/2）,Ws/（Fs/2）,Rp,Rs）;
[b，a］=cheby1(N，Rp,Wn);
freqz（b，a,512,1000);

例5。3：设计一个切比雪夫Ⅱ型数字带通滤波器，要求带通范围100—250Hz，带阻上限为300Hz，下限为50Hz,通带内纹波小于3dB，阻带纹波为30 dB，抽样频率为1000 Hz,并利用最小的阶次实现。

clear all;
Wpl=100;
Wph=250；
Wp=[Wpl，Wph];
Rp=3;
Wsl=50;
Wsh=300；
Ws=[Wsl，Wsh］；
Rs=30；
Fs=1000；
[N，Wn]=cheb2ord(Wp/(Fs/2)，Ws/（Fs/2),Rp，Rs）；
[b，a]=cheby2（N,Rp,Wn);
freqz（b,a,512,1000);

（六)、FIR数字滤波器设计
1、、在MATLAB 中产生窗函数十分简单：
（1）矩形窗（Rectangle Window）
调用格式:w=boxcar（n)，根据长度n 产生一个矩形窗w.
（2)三角窗(Triangular Window)
调用格式:w=triang（n） ，根据长度n 产生一个三角窗w.
（3)汉宁窗(Hanning Window）
调用格式：w=hanning(n) ，根据长度n 产生一个汉宁窗w.
(4）海明窗（Hamming Window）
调用格式：w=hamming（n) ,根据长度n 产生一个海明窗w.
（5)布拉克曼窗(Blackman Window）
调用格式:w=blackman(n) ，根据长度n 产生一个布拉克曼窗w。
（6)恺撒窗（Kaiser Window）
调用格式：w=kaiser（n，beta) ，根据长度n 和影响窗函数旁瓣的β参数产生一个恺撒窗w。

2、基于窗函数的FIR 滤波器设计
利用MATLAB 提供的函数firl 来实现
调用格式:firl （n,Wn,’ftype’,Window),n 为阶数、Wn 是截止频率（如果输入是形如［W1 W2］
的矢量时，本函数将设计带通滤波器，其通带为W1<ω〈W2）、ftype 是滤波器的类型(低
通-省略该参数、高通-ftype=high、带阻-ftype=stop）、Window 是窗函数.
例6.1： 设计一个长度为8 的线性相位FIR 滤波器.其理想幅频特性满足

Window=boxcar（8)；
b=fir1(7，0.4，Window）；
freqz(b,1）

例6。2:设计线性相位带通滤波器，其长度N=15，上下边带截止频率分别为W1= 0。3π，w2=0。5
π
Window=blackman（16）；
b=fir1（15，[0。3 0。5]，Window）；
freqz(b，1）

例6.3：MATLAB中的chirp。mat文件中存储信号ｙ的数据，该信号的大部分号能量集中在Fs/4（或二分之一奈奎斯特）以上，试设计一个34阶的FIR高通滤波器,滤除频率低于Fs/4的信号成分,其中滤波器的截止频率为0.48，阻带衰减为30dB，滤波器窗采用切比雪夫窗
clear all；
load chirp
window=chebwin(35，30）；
b=fir1（34，0.48,'high’，window）；
yfit=filter(b,1,y）;
［Py，fy］=pburg(y,10,512,Fs）;
［Pyfit，fyfit]=pburg（yfit,10,512，Fs);
plot（fy,10＊log10（Py）,'.’，fyfit， 10*log10(Pyfit));
grid on
ylabel(‘幅度（dB)’)
xlabel（‘频率(Hz’）
legend(‘滤波前的线性调频信号’， ‘滤波后的线性调频信号'）