《MATLAB辅助现代工程数字信号处理》课件第5章

第5章多采样率数字信号处理5.1信号的抽取与内插

5.2有理因子采样率转换

5.3采样率转换的滤波器实现

5.4数字滤波器组

5.5小结

5.1信号的抽取与内插

5.1.1信号的整数倍抽取

在时域中，若原始信号为x(n)，抽取因子为D，抽取后得到的信号为y(n)，则整个信号的抽取过程可表示为

y(n)＝x(Dn)

(5.1)

信号的抽取过程可分为两个步骤:第一步，将x(n)与一个周期为D的采样脉冲序列相乘，即每D个点中保留一个点，其他(D－1)个点为零，以便得到采样信号ω(n);第二步，去掉ω(n)中的零点后，便得到一个低速率的信号y(n)。中间采样信号ω(n)为

(5.2)抽取后的信号为

y(n)=ω(Dn)

(5.3)

【例5.1】

设x(n)=2sin(0.089πn)，抽取因子D为3，求抽取输出信号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)′);

程序运行结果如图5.1所示。图5.1信号的整数倍抽取时域分析为分析下行抽样信号的频谱，先计算ω(n)的频谱为(5.4)则抽取信号的频谱为(5.5)

抽取信号的频谱与原来信号的频谱有以下关系:首先X(ω)作(D－1)次及等间隔平移，其平移间隔为2π/D，然后作叠加平均得到W(ω)，最后频谱拉伸D倍后即可得到抽取信号的频谱。

【例5.2】

利用MATLAB的函数fir2构造一个带限输入序列:f=[0，0.45，0.49，1];幅值mag=[0，1，0，0]。试分析信号抽取因子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)]);%产生抽取输出谱

%求取并绘制抽取输出谱

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

subplot(212);

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

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

程序运行结果如图5.2所示。图5.2整数倍抽取频域分析5.1.2信号的整数倍内插

将采样率增加到I倍可以通过在信号的两个连续值之间内插(I－1)个样本来实现。内插过程可以由许多方式实现，下面介绍的是一个保持信号序列x(n)的谱形状不变的方式。在时域中，若已知输入信号x(n)，内插因子为I，插值后得到的信号为y(n)，则整个信号的内插过程可表示为(5.6)其他

内插过程是由“填零”方式进行的上行采样的过程。

【例5.3】

设x(n)=2sin(0.26πn)，内插因子I为3，求内插输出信号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)′);

程序运行结果如图5.3所示。图5.3信号的整数倍内插时域分析对输出信号进行傅里叶变换，得到信号y(n)的频谱为(5.7)

【例5.4】

利用MATLAB的函数fir2构造一个带限输入序列:f=[0，0.46，0.51，1];幅值mag=[0，1，0，0]。试分析信号内插因子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;

%产生抽取输出谱

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;

程序运行结果如图5.4所示。图5.4内插过程的频域分析

5.2有理因子采样率转换

对于有理因子I/D采样率转换的情形，可以通过抽取和内插串联来实现。先通过内插因子I插入信号，再利用抽取因子D来抽取信号。内插过程的平滑滤波器和抽取过程的反混叠滤波器可以合并为一个低通滤波器，理想情况下的频率响应为(5.8)其他图5.5给出了采样率通过有理因子I/D改变的系统框图。图5.5有理因子采样率转换的系统框图对于上述有理因子采样率的转换，输出信号的频谱可表示为(5.9)其他

5.3采样率转换的滤波器实现

5.3.1抽取采样率转换的滤波器实现

为消除抽取操作引起的混叠，在执行抽样之前，输入信号必须通过低通滤波器H(z)将带宽限制到|ω|<。因此，实际的信号抽取系统如图5.6所示。图5.6抽取过程的实际结构该系统一般称为抽取器。由图可知:

(5.10)(5.11)因此(5.12)

若“反混叠”滤波器为理想的滤波器(5.13)其他则有(5.14)

MATLAB信号处理工具箱提供抽取函数decimate用于信号的整数倍抽取，调用格式为

y=decimate(x，D)

y=decimate(x，D，n)

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

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

其中，x为输入信号;D为抽取因子;′fir′为指定的FIR滤波器，缺省时采用切比雪夫Ⅰ型低通滤波器压缩频带;n为低通滤波器的阶数，当采用FIR滤波器时，缺省时为30点数，否则缺省为8阶。

【例5.5】

调频信号x(t)=2sin(kπt2)，k=1;0≤t≤T;T=4;采样率为fs，采样点数N=Tfs;fs=4fc。利用MATLAB编程，分析抽取因子分别为2和4时的情况。

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);%原始输入信号x(n)

figure(1);

subplot(221);

stem(t*Ts，x);

D=2;%抽取因子为2

y=decimate(x，D);%抽取输出信号y(n)

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

subplot(223);

stem(tnew*D*Ts，y);

X=fft(x);%原始输入信号频谱X(ω)

X=fftshift(X);

subplot(222);

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

Y=fft(y);%抽取输出信号频谱Y(ω)

Y=fftshift(Y);

subplot(224);

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

figure(2);

D=4;%抽取因子为4

y=decimate(x，D);%抽取输出信号y(n)

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

subplot(221);

stem(tnew*D*Ts，y);

Y=fft(y);%抽取输出信号频谱Y(ω)

Y=fftshift(Y);

subplot(222);

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

程序运行结果如图5.7所示。图5.7信号抽取过程的时频分析5.3.2内插采样率转换的滤波器实现

为消除内插“填零”过程中引起的“复制”，一般采用平滑滤波来处理。内插过程的实际结构如图5.8所示。图5.8内插过程的实际结构上述系统一般称为内插器。若滤波采用理想低通滤波器(5.15)其他又因(5.16)故有(5.17)

【例5.6】

调频信号x(t)=2cos(kπt2)，k=1;0≤t≤T;T=4;采样率为fs，采样点数N=Tfs;fs=2.5fc。利用MATLAB编程，分析将采样率提高3倍时的情况。

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);%原始输入信号x(n)

subplot(221);

stem(t*Ts，x);

I=3;

y=interp(x，I);%内插输出信号y(n)

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

subplot(223);

stem(tnew*Ts/I，y);

X=fft(x);%原始输入信号频谱X(ω)

X=fftshift(X);

subplot(222);

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

Y=fft(y);%内插输出信号频谱Y(ω)

Y=fftshift(Y);

subplot(224);

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

程序运行结果如图5.9所示。图5.9信号内插过程的时频分析5.3.3有理因子采样率转换的滤波器实现

利用一个内插器连接一个抽取采样器，就可以有效地实现一个比率为I/D的采样率转换器。MATLAB信号处理工具箱提供了重采样函数resample用于有理因子的采样率转换，调用格式为

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】

调频信号x(t)=3sin(kπt2)，k=1;0≤t≤T;T=4;采样率为fs，采样点数N=Tfs;fs=3fc。利用MATLAB编程，分析I=7，D=3时有理因子的采样率转换。

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);%原始输入信号x(n)

subplot(221);

stem(t*Ts，x);

I=7;D=3;

y=resample(x，I，D);%输出信号y(n)

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

subplot(223);

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

X=fft(x);%原始输入信号频谱X(ω)

X=fftshift(X);

subplot(222);

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

Y=fft(y);%输出信号频谱Y(ω)

Y=fftshift(Y);

subplot(224);

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

程序运行结果如图5.10所示。图5.10信号有理因子转换的时频分析

5.4数字滤波器组

数字滤波器组是带有共同输入或相加输出的一组数字带通滤波器。一般来说，滤波器组可以分为两类:分析滤波器组和综合滤波器组。分析滤波器组由一组系统函数为{Hk(z)}的滤波器按图5.11(a)排列成并行组构成，该滤波器组的频率响应特性将信号分成相应个数的子带。综合滤波器组由一组系统函数为{Gk(z)}、相应输入为{yk(n)}的滤波器按图5.11(b)排列组成，将各滤波器的输出相加起来构成信号{x(n)}。滤波器组经常用来实现谱分析和信号综合。图5.11数字滤波器组5.4.1均匀滤波器组

当在序列{x(n)}的离散傅里叶变换的计算中使用滤波器组时，称为DFT滤波器组。由N个滤波器{Hk(z)，k=0，1，2，…，N－1}组成的分析滤波器组称为均匀DFT分析滤波器组，如果滤波器系统函数{Hk(z)}是从原型滤波器H0(z)导出的，且满足

(5.18)则滤波器{Hk(z)}的频率特性可以简单地通过以2π/N的倍数均匀地平移原型滤波器的频域响应得到。在时域中，利用冲激响应表征滤波器，可以表示为(5.19)式中，h0(n)是原型滤波器的冲激响应，一般来说是FIR滤波器或IIR滤波器。如果H0(z)表示原型滤波器的系统函数，那么第k个滤波器的系统函数为

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

(5.20)

【例5.8】

利用MATLAB设计一个均匀分析滤波器组，滤波器原型为20阶低通滤波器，且满足:归一化频率向量f=[0,0.3,0.4,1];幅值向量a=[1,1,0,0];权重w=[10，1]。

MATLAB程序如下:

%MATLABPROGRAM5-8

clf;

%设计原型低通滤波器

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

%画出每个滤波器的幅度响应

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]);

程序运行结果如图5.12所示。图5.12均匀分析滤波器组的设计下面探讨均匀分析滤波器组与综合滤波器组的具体实现。对于均匀分析滤波器组来说，对每个带通滤波器