《MATLAB輔助現(xiàn)代工程數(shù)字信號(hào)處理》課件第5章

第5章多采樣率數(shù)字信號(hào)處理5.1信號(hào)的抽取與內(nèi)插

5.2有理因子采樣率轉(zhuǎn)換

5.3采樣率轉(zhuǎn)換的濾波器實(shí)現(xiàn)

5.4數(shù)字濾波器組

5.5小結(jié)

5.1信號(hào)的抽取與內(nèi)插

5.1.1信號(hào)的整數(shù)倍抽取

在時(shí)域中，若原始信號(hào)為x(n)，抽取因子為D，抽取后得到的信號(hào)為y(n)，則整個(gè)信號(hào)的抽取過程可表示為

y(n)＝x(Dn)

(5.1)

信號(hào)的抽取過程可分為兩個(gè)步驟:第一步，將x(n)與一個(gè)周期為D的采樣脈沖序列相乘，即每D個(gè)點(diǎn)中保留一個(gè)點(diǎn)，其他(D－1)個(gè)點(diǎn)為零，以便得到采樣信號(hào)ω(n);第二步，去掉ω(n)中的零點(diǎn)后，便得到一個(gè)低速率的信號(hào)y(n)。中間采樣信號(hào)ω(n)為

(5.2)抽取后的信號(hào)為

y(n)=ω(Dn)

(5.3)

【例5.1】

設(shè)x(n)=2sin(0.089πn)，抽取因子D為3，求抽取輸出信號(hào)y(n)。

MATLAB程序如下:

%MATLABPROGRAM5－1

clc;

n=[0:49];

m=[0:50*3－1];

x=2*sin(0.089*pi*m);

y=x([1:3:length(x)]);

subplot(211);

stem(n，x(1:50));xlabel(′n′);ylabel(′x(n)′);

subplot(212);

stem(n，y);xlabel(′n′);ylabel(′y(n)′);

程序運(yùn)行結(jié)果如圖5.1所示。圖5.1信號(hào)的整數(shù)倍抽取時(shí)域分析為分析下行抽樣信號(hào)的頻譜，先計(jì)算ω(n)的頻譜為(5.4)則抽取信號(hào)的頻譜為(5.5)

抽取信號(hào)的頻譜與原來信號(hào)的頻譜有以下關(guān)系:首先X(ω)作(D－1)次及等間隔平移，其平移間隔為2π/D，然后作疊加平均得到W(ω)，最后頻譜拉伸D倍后即可得到抽取信號(hào)的頻譜。

【例5.2】

利用MATLAB的函數(shù)fir2構(gòu)造一個(gè)帶限輸入序列:f=[0，0.45，0.49，1];幅值mag=[0，1，0，0]。試分析信號(hào)抽取因子D=2的頻域特性。

MATLAB程序如下:

%MATLABPROGRAM5-2

freq=[0，0.45，0.49，1];

mag=[0，1，0，0];

x=fir2(101，freq，mag);

%求取并繪制輸入譜

[Xz，w]=freqz(x，1，512);

subplot(211);

plot(w/pi，abs(Xz));

xlabel(′＼omega/＼pi′);ylabel(′|Xz(w)|′);grid;

y=x([1:2:length(x)]);%產(chǎn)生抽取輸出譜

%求取并繪制抽取輸出譜

[Yz，w]=freqz(y，1，512);

subplot(212);

plot(w/pi，abs(Yz));

xlabel(′＼omega/＼pi′);ylabel(′|Yz(w)|′);grid;

程序運(yùn)行結(jié)果如圖5.2所示。圖5.2整數(shù)倍抽取頻域分析5.1.2信號(hào)的整數(shù)倍內(nèi)插

將采樣率增加到I倍可以通過在信號(hào)的兩個(gè)連續(xù)值之間內(nèi)插(I－1)個(gè)樣本來實(shí)現(xiàn)。內(nèi)插過程可以由許多方式實(shí)現(xiàn)，下面介紹的是一個(gè)保持信號(hào)序列x(n)的譜形狀不變的方式。在時(shí)域中，若已知輸入信號(hào)x(n)，內(nèi)插因子為I，插值后得到的信號(hào)為y(n)，則整個(gè)信號(hào)的內(nèi)插過程可表示為(5.6)其他

內(nèi)插過程是由“填零”方式進(jìn)行的上行采樣的過程。

【例5.3】

設(shè)x(n)=2sin(0.26πn)，內(nèi)插因子I為3，求內(nèi)插輸出信號(hào)y(n)。

MATLAB程序如下:

%MATLABPROGRAM5-3

n=[0:50];

x=2*sin(0.26*pi*n);

y=zeros(1，3*length(x));

y([1:3:length(y)])=x;

subplot(211);

stem(n，x);

xlabel(′n′);ylabel(′x(n)′);

subplot(212);

stem(n，y(1:length(x)));

xlabel(′n′);ylabel(′y(n)′);

程序運(yùn)行結(jié)果如圖5.3所示。圖5.3信號(hào)的整數(shù)倍內(nèi)插時(shí)域分析對(duì)輸出信號(hào)進(jìn)行傅里葉變換，得到信號(hào)y(n)的頻譜為(5.7)

【例5.4】

利用MATLAB的函數(shù)fir2構(gòu)造一個(gè)帶限輸入序列:f=[0，0.46，0.51，1];幅值mag=[0，1，0，0]。試分析信號(hào)內(nèi)插因子I=2的頻域特性。

MATLAB程序如下:

%MATLABPROGRAM5-4

freq=[0，0.46，0.51，1];

mag=[0，1，0，0];

I=2;

x=fir2(99，freq，mag);

%求取并繪制輸入譜

[Xz，w]=freqz(x，1，512，′whole′);

subplot(211);

plot(w/pi，abs(Xz));

xlabel(′＼omega/＼pi′);ylabel(′|Xz(w)|′);grid;

%產(chǎn)生抽取輸出譜

y=zeros(1，I*length(x));

y([1:I:length(y)])=x;

%求取并繪制抽取輸出譜

[Yz，w]=freqz(y，1，512，′whole′);

subplot(212);

plot(w/pi，abs(Yz));

xlabel(′＼omega/＼pi′);ylabel(′|Yz(w)|′);grid;

程序運(yùn)行結(jié)果如圖5.4所示。圖5.4內(nèi)插過程的頻域分析

5.2有理因子采樣率轉(zhuǎn)換

對(duì)于有理因子I/D采樣率轉(zhuǎn)換的情形，可以通過抽取和內(nèi)插串聯(lián)來實(shí)現(xiàn)。先通過內(nèi)插因子I插入信號(hào)，再利用抽取因子D來抽取信號(hào)。內(nèi)插過程的平滑濾波器和抽取過程的反混疊濾波器可以合并為一個(gè)低通濾波器，理想情況下的頻率響應(yīng)為(5.8)其他圖5.5給出了采樣率通過有理因子I/D改變的系統(tǒng)框圖。圖5.5有理因子采樣率轉(zhuǎn)換的系統(tǒng)框圖對(duì)于上述有理因子采樣率的轉(zhuǎn)換，輸出信號(hào)的頻譜可表示為(5.9)其他

5.3采樣率轉(zhuǎn)換的濾波器實(shí)現(xiàn)

5.3.1抽取采樣率轉(zhuǎn)換的濾波器實(shí)現(xiàn)

為消除抽取操作引起的混疊，在執(zhí)行抽樣之前，輸入信號(hào)必須通過低通濾波器H(z)將帶寬限制到|ω|<。因此，實(shí)際的信號(hào)抽取系統(tǒng)如圖5.6所示。圖5.6抽取過程的實(shí)際結(jié)構(gòu)該系統(tǒng)一般稱為抽取器。由圖可知:

(5.10)(5.11)因此(5.12)

若“反混疊”濾波器為理想的濾波器(5.13)其他則有(5.14)

MATLAB信號(hào)處理工具箱提供抽取函數(shù)decimate用于信號(hào)的整數(shù)倍抽取，調(diào)用格式為

y=decimate(x，D)

y=decimate(x，D，n)

y=decimate(x，D，′fir′)

y=decimate(x，D，n，′fir′)

其中，x為輸入信號(hào);D為抽取因子;′fir′為指定的FIR濾波器，缺省時(shí)采用切比雪夫Ⅰ型低通濾波器壓縮頻帶;n為低通濾波器的階數(shù)，當(dāng)采用FIR濾波器時(shí)，缺省時(shí)為30點(diǎn)數(shù)，否則缺省為8階。

【例5.5】

調(diào)頻信號(hào)x(t)=2sin(kπt2)，k=1;0≤t≤T;T=4;采樣率為fs，采樣點(diǎn)數(shù)N=Tfs;fs=4fc。利用MATLAB編程，分析抽取因子分別為2和4時(shí)的情況。

MATLAB程序如下:

%MATLABPROGRAM5-5

k=1;

T=4;

fc=k*T;fs=4*fc;

Ts=1/fs;N=T/Ts;

x=zeros(1，N);

t=[0:N-1];

x=2*sin(k*pi*(t*Ts).^2);%原始輸入信號(hào)x(n)

figure(1);

subplot(221);

stem(t*Ts，x);

D=2;%抽取因子為2

y=decimate(x，D);%抽取輸出信號(hào)y(n)

tnew=[0:N/D－1];

subplot(223);

stem(tnew*D*Ts，y);

X=fft(x);%原始輸入信號(hào)頻譜X(ω)

X=fftshift(X);

subplot(222);

plot((t-N/2)*fs/N，abs(X));

Y=fft(y);%抽取輸出信號(hào)頻譜Y(ω)

Y=fftshift(Y);

subplot(224);

plot((tnew－N/D/2)*fs/N，abs(Y));

figure(2);

D=4;%抽取因子為4

y=decimate(x，D);%抽取輸出信號(hào)y(n)

tnew=[0:N/D－1];

subplot(221);

stem(tnew*D*Ts，y);

Y=fft(y);%抽取輸出信號(hào)頻譜Y(ω)

Y=fftshift(Y);

subplot(222);

plot((tnew-N/D/2)*fs/N，abs(Y));

程序運(yùn)行結(jié)果如圖5.7所示。圖5.7信號(hào)抽取過程的時(shí)頻分析5.3.2內(nèi)插采樣率轉(zhuǎn)換的濾波器實(shí)現(xiàn)

為消除內(nèi)插“填零”過程中引起的“復(fù)制”，一般采用平滑濾波來處理。內(nèi)插過程的實(shí)際結(jié)構(gòu)如圖5.8所示。圖5.8內(nèi)插過程的實(shí)際結(jié)構(gòu)上述系統(tǒng)一般稱為內(nèi)插器。若濾波采用理想低通濾波器(5.15)其他又因(5.16)故有(5.17)

【例5.6】

調(diào)頻信號(hào)x(t)=2cos(kπt2)，k=1;0≤t≤T;T=4;采樣率為fs，采樣點(diǎn)數(shù)N=Tfs;fs=2.5fc。利用MATLAB編程，分析將采樣率提高3倍時(shí)的情況。

MATLAB程序如下:

%MATLABPROGRAM5-6

k=1;

T=4;

fc=k*T;fs=2.5*fc;

Ts=1/fs;N=T/Ts;

x=zeros(1，N);

t=[0:N－1];

x=2*cos(k*pi*(t*Ts).^2);%原始輸入信號(hào)x(n)

subplot(221);

stem(t*Ts，x);

I=3;

y=interp(x，I);%內(nèi)插輸出信號(hào)y(n)

tnew=[0:N*I－1];

subplot(223);

stem(tnew*Ts/I，y);

X=fft(x);%原始輸入信號(hào)頻譜X(ω)

X=fftshift(X);

subplot(222);

plot((t－N/2)*fs/N，abs(X));

Y=fft(y);%內(nèi)插輸出信號(hào)頻譜Y(ω)

Y=fftshift(Y);

subplot(224);

plot((tnew-N*I/2)*fs/N，abs(Y));

程序運(yùn)行結(jié)果如圖5.9所示。圖5.9信號(hào)內(nèi)插過程的時(shí)頻分析5.3.3有理因子采樣率轉(zhuǎn)換的濾波器實(shí)現(xiàn)

利用一個(gè)內(nèi)插器連接一個(gè)抽取采樣器，就可以有效地實(shí)現(xiàn)一個(gè)比率為I/D的采樣率轉(zhuǎn)換器。MATLAB信號(hào)處理工具箱提供了重采樣函數(shù)resample用于有理因子的采樣率轉(zhuǎn)換，調(diào)用格式為

y=resample(x，I，D)

y=resample(x，I，D，n)

y=resample(x，I，D，n，beta)

y=resample(x，I，D，b)

[y，b]=resample(x，I，D)

【例5.7】

調(diào)頻信號(hào)x(t)=3sin(kπt2)，k=1;0≤t≤T;T=4;采樣率為fs，采樣點(diǎn)數(shù)N=Tfs;fs=3fc。利用MATLAB編程，分析I=7，D=3時(shí)有理因子的采樣率轉(zhuǎn)換。

MATLAB程序如下:

%MATLABPROGRAM5-7

k=1;

T=4;

fc=k*T;fs=3*fc;

Ts=1/fs;N=T/Ts;

x=zeros(1，N);

t=[0:N－1];

x=3*sin(k*pi*(t*Ts).^2);%原始輸入信號(hào)x(n)

subplot(221);

stem(t*Ts，x);

I=7;D=3;

y=resample(x，I，D);%輸出信號(hào)y(n)

tnew=[0:N*I/D－1];

subplot(223);

stem(tnew*Ts/I*D，y);

X=fft(x);%原始輸入信號(hào)頻譜X(ω)

X=fftshift(X);

subplot(222);

plot((t－N/2)*fs/N，abs(X));

Y=fft(y);%輸出信號(hào)頻譜Y(ω)

Y=fftshift(Y);

subplot(224);

plot((tnew－N*I/D/2)*fs/N，abs(Y));

程序運(yùn)行結(jié)果如圖5.10所示。圖5.10信號(hào)有理因子轉(zhuǎn)換的時(shí)頻分析

5.4數(shù)字濾波器組

數(shù)字濾波器組是帶有共同輸入或相加輸出的一組數(shù)字帶通濾波器。一般來說，濾波器組可以分為兩類:分析濾波器組和綜合濾波器組。分析濾波器組由一組系統(tǒng)函數(shù)為{Hk(z)}的濾波器按圖5.11(a)排列成并行組構(gòu)成，該濾波器組的頻率響應(yīng)特性將信號(hào)分成相應(yīng)個(gè)數(shù)的子帶。綜合濾波器組由一組系統(tǒng)函數(shù)為{Gk(z)}、相應(yīng)輸入為{yk(n)}的濾波器按圖5.11(b)排列組成，將各濾波器的輸出相加起來構(gòu)成信號(hào){x(n)}。濾波器組經(jīng)常用來實(shí)現(xiàn)譜分析和信號(hào)綜合。圖5.11數(shù)字濾波器組5.4.1均勻?yàn)V波器組

當(dāng)在序列{x(n)}的離散傅里葉變換的計(jì)算中使用濾波器組時(shí)，稱為DFT濾波器組。由N個(gè)濾波器{Hk(z)，k=0，1，2，…，N－1}組成的分析濾波器組稱為均勻DFT分析濾波器組，如果濾波器系統(tǒng)函數(shù){Hk(z)}是從原型濾波器H0(z)導(dǎo)出的，且滿足

(5.18)則濾波器{Hk(z)}的頻率特性可以簡單地通過以2π/N的倍數(shù)均勻地平移原型濾波器的頻域響應(yīng)得到。在時(shí)域中，利用沖激響應(yīng)表征濾波器，可以表示為(5.19)式中，h0(n)是原型濾波器的沖激響應(yīng)，一般來說是FIR濾波器或IIR濾波器。如果H0(z)表示原型濾波器的系統(tǒng)函數(shù)，那么第k個(gè)濾波器的系統(tǒng)函數(shù)為

Hk(z)=H0(ze－j2πk/N)，1<k<N－1

(5.20)

【例5.8】

利用MATLAB設(shè)計(jì)一個(gè)均勻分析濾波器組，濾波器原型為20階低通濾波器，且滿足:歸一化頻率向量f=[0,0.3,0.4,1];幅值向量a=[1,1,0,0];權(quán)重w=[10，1]。

MATLAB程序如下:

%MATLABPROGRAM5-8

clf;

%設(shè)計(jì)原型低通濾波器

b=remez(20，[0，0.2，0.25，1]，[1，1，0，0]，[10，1]);

w=[0:2*pi/255:2*pi];n=[0:20];

fork=[1:4];

c=exp(2*pi*(k-1)*n*i/4);

FB=b.*c;

HB(k，:)=freqz(FB，1，w);

end

%畫出每個(gè)濾波器的幅度響應(yīng)

subplot(221);

plot(w/pi，abs(HB(1，:)));

xlabel(′＼omega/＼pi′);ylabel(′振幅′);

title(′濾波器No.1′);axis([0，2，0，1.1]);

subplot(222);

plot(w/pi，abs(HB(2，:)));

xlabel(′＼omega/＼pi′);ylabel(′振幅′);

title(′濾波器No.2′);axis([0，2，0，1.1]);

subplot(223);

plot(w/pi，abs(HB(3，:)));

xlabel(′＼omega/＼pi′);ylabel(′振幅′);

title(′濾波器No.3′);axis([0，2，0，1.1]);

subplot(224);

plot(w/pi，abs(HB(4，:)));

xlabel(′＼omega/＼pi′);ylabel(′振幅′);

title(′濾波器No.4′);axis([0，2，0，1.1]);

程序運(yùn)行結(jié)果如圖5.12所示。圖5.12均勻分析濾波器組的設(shè)計(jì)下面探討均勻分析濾波器組與綜合濾波器組的具體實(shí)現(xiàn)。對(duì)于均勻分析濾波器組來說，對(duì)每個(gè)帶通濾波器