发布时间 : 星期二 文章数字信号处理上机指导 - 1dl更新完毕开始阅读a620bd4b852458fb770b56fc
第5章 常用数字信号处理函数
§5-1 波形发生器
1．基本信号序列的表示
离散时间信号是时间上不连续的“序列”，一般直接用x(n)表示。一个长度为N的序列x(n)在MATLAB中可以用一个N维行向量或列向量表示。常用的基本信号序列产生如下：
1n?01） 单位抽样序列 ?(n)???0 n?0
?这一序列可用MATLAB中的zeros函数实现
x=[1,zeros(1,N-1)];
1n?02） 单位阶跃序列 u(n)???0 n?0
?这一序列可用MATLAB中的ones函数实现
x=ones(1,N); 3） 实指数序列 x(n)?an ?n,a?R
MATLAB实现：
n = 0 ; N-1;
x = a.^ n ;
4） 复指数序列 x(n)?e(c?j?)n ?n
MATLAB实现：
n=0 : N-1;
x = exp((c+j*w)*n) ; 5） 随机序列
rand(1, N)产生在[0,1]上均匀分布的长度为N的随机序列。
randn(1, N)产生长度为N且均值为0方差为1的高斯随机序列，即白噪声序列。 2．基本波形的产生
1）三角波或锯齿波发生函数：sawtooth()
语法格式：sawtooth(t ,width)。产生以2π为周期幅值范围在[-1 , +1]之间的三角波或锯齿波。参数t为时间向量；width是[0，1]之间的数，它决定于函数在一个周期内上升部分和下降部分的比例。width=0.5产生三角波，width=1产生锯齿波，此时函数可简写为sawtooth(t)。 2）方波发生函数square()
语法格式：square(t)。产生以2π为周期幅值范围在[-1，+1]之间的方波，参数t为时间向量。
3）sinc发生函数：sinc()
语法格式：sinc(t)
?sin(pi*t)/(pi*t);t?0 sinc(t)??
t?01?例 信号产生举例 function fun1 clear all
t=0:0.0001:0.1; figure(1);
- 12 -
x1=sawtooth(2*pi*50*t); %在[0,0.1]时间内产生5个周期的锯齿波 subplot(2,2,1) plot(t,x1)
x2=sawtooth(2*pi*50*t,0.5); %在[0,0.1]时间内产生5个周期的三角波 subplot(2,2,2) plot(t,x2)
x3=square(2*pi*50*t); %在[0,0.1]时间内产生5个周期的方波 subplot(2,2,3) plot(t,x3)
axis([0,0.1,-1.2,1.2]) t=-4:0.1:4;
x4=sinc(t); %产生抽样函数 subplot(2,2,4) plot(t,x4)
运行结果如下上图。 §5-2 序列的操作
在数字信号处理中基本的序列运算及其MATLAB实现如下表5-1所示：
表5-1一些常用的序列运算及其matlab实现
信号加 信号乘 数乘 折叠 抽样和 抽样积 信号能量 序列运算 x(n)=x1(n)+x2(n) x(n)=x1(n)×x2(n) y(n)=k×x(n) y(n)=x(-n) n2其MATLAB实现 x=x1+x2 x=x1.*x2 y=k*x y=fliplr(x) y=sum(x(n1:n2)) y=prod(x(n1:n2)) Ex=sum(abs(x).^2)) y?y??x(n) n?n1n2?n?n1x(n) 2?Ex??n???x(n) 2信号功率 §5-3 常用数字信号处理函数
Px?1NN?1?n?0x(n) Px=sum(abs(x).^2))/n 对于线性离散系统，其表达方式有多种： 传递函数法 H(z)?b0?b1z?1?1???bMz???aNz?M?Na0?a1z
零极点增益法 H(z)?K(z?y1)(z?y2)?(z?yM)(z?x1)(z?x2)?(z?xN)带余式的部分分式展开法
rNr1r2 H(z)??????k1?k2z?1???kM?N?1z?(M?N)
?1?1?11?p1z1?p2z1?pNzNMk差分方程法 y(n)??ak?1y(n?k)?m?0?bmx(n?m)
数字信号处理即研究信号通过系统后的变化，其中经常用到的函数有：
- 13 -
1． x=roots(a)
欲求多项式y(z)=a0+a1z –1+a2z –2+? 的根，可用命令x=roots(a)来实现，其中a=[a0 , a1 , a2 , ?]是它的系数向量，x便是y(z)的根向量。
例 求多项式y(z)=1-1.6z –1+0.65z –2-0.05z –3的根。编程运行如下： >> a=[1, -1.6 , 0.65 , -0.05] ; x=roots(a)
x =
1.0000 0.5000
0.1000
从而原多项式可表达为y(z)=(1-z –1)(1-0.5z –1)(1-0.1z –1)。 2． b=poly(x)
欲求出y(z)=(z-x1)(z-x2)?的多项式表达，可用命令b=poly(x)来实现，其中x=[x1 , x2 , ?]是它的根向量，b是欲求多项式的系数向量。 例 某滤波器的零极点增益法表达式为H(z)?(z?0.1)(z?0.5)(z?1)，求其传递函数表达
(z?0.3)(z?2)(z?5)式。MATLAB程序及运行结果如下：
function [a,b]=fun2(x,y,k) b=k*poly(y); a=poly(x);
>> y=[0.1,0.5,1]; x=[0.3,2,5]; k=1; [ b, a]=fun2(y,x,k) b =
1.0000 -1.6000 0.6500 -0.0500 a =
1.0000 -7.3000 12.1000 -3.0000
?1?2?3从而滤波器的传递函数表达式为H(z)?1?1.6z?0.65z?0.05z 。
?1?2?31?7.3z?12.1z?3.0z3．[r , p , k] = residuez(b , a)
由传递函数表达式求出带余式的部分分式展开式时可用命令[r, p , k] = residuez(b , a) 来实现，其中b、a分别是原系统传递函数表达式中的分子、分母系数向量，而r、p、k分别是该系统的带余式的部分分式展开式中的系数向量。 例 对 H(z)??4?8z?1作部分分式展开。程序和显示输出如下： 1?6z?1?8z?2>> b=[-4 , 8] ; a=[1 , 6 , 8] ; [r , p , k]=residuez(b , a) r = -12 8 p = -4 -2 k =
[ ]
这就是说，系统的部分分式展开形式为：
H(z)??128。 ??1?11?4z1?2z?4．y (n) = conv (x , h)
求两个序列的线性卷积 y(n)?x(n)*h(n)?k????x(k)h(n?k) 可直接采用命令y (n) =
conv (x , h )来实现。其中特别要注意的是，函数conv默认序列从n = 0开始。若{x (n) : nx1 ≤
n ≤ nx2} , {h(n) : nh1 ≤ n ≤ nh2}，则卷积结果{y (n) : ny1 ≤ n ≤ ny2}，其中ny1= nx1+ nh1 , ny2= nx2+ ny2。
例 已知x (n) = {3 , 4 , 5} , h (n) = {2 , 6 , 7 , 8}，求它们的线性卷积。
>> x=[3,4,5]; h=[2,6,7,8]; n1=-1:1; n2=-1:2; y=conv(x,h) y =
- 14 -
6 26 55 82 67 40 >>nb3=n1(1)+n2(1);
>>ne3=n1(length(x))+n2(length(h)); >>n3=nb3:ne3 n3 =
-2 -1 0 1 2 3
即序列x(n)与h(n)的线性卷积为y(n)={6 , 26 , 55 , 82 , 67 , 40}。 5．X=fft(x , N) 与x=ifft(X , N)
N?1采用FFT算法计算一维序列x(n)的N点离散Fourier变换(DFT) X(k)??x(n)Wn?0nkN可以
直接用命令X=fft (x , N)来实现，其逆运算用命令x = ifft(X , N)来实现。计算中若序列x（或X）的长度大于N时，截断x（或X）；否则补零。
例 模拟信号x(t)=2sin(4πt)+5cos(8πt)，求其64点的DFT的幅值谱及相位谱。 MATLAB实现：function y=fun3(N)
n=0:N-1; t=0.01*n; q=n*2*pi/N; x=fun30(t); y=fft(x,N); subplot(3,1,1);
plot(t,x); title('原信号'); grid subplot(3,1,2);
plot(q,abs(y)); title('幅值'); grid subplot(3,1,3);
plot(q,angle(y));title('相位'); grid Function x=fun30(t)
x=2*sin(4*pi*t)+5*cos(8*pi*t); >>y=fun3(64);
6．[H , W] = freqz (b , a , n , Fs)
函数freqz基于FFT算法计算以传递函数形式表达的数字滤波器的频率响应。命令freqz (b , a , N , Fs)能够自动绘制以b、a为系数的数字滤波器的N点幅频和相频曲线(Fs是抽样频率，此项可缺省)，命令[H , w] = freqz (b , a , n)计算并返回数字滤波器的N点频率响应H (abs ( H ) )是幅频响应，angle (H)是相频响应），相应的N个分布在(0~π)范围内的频率点记录在w中。
例 绘制系统H(z)?1的幅频和相频特性曲线。
1?0.9z?1>> b=1; a=[1,-0.9]; freqz(b , a , 512)
- 15 -