第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 -