"实验一二维抽样定理

一 的

通过计算机仿真掌握二维抽样定理(惠特克 香农抽样定理，Whittaker Shanno samplin theorem)，包括带限函数的定义?连续函数的离散化(抽样过程)，以及利用抽样函数重构原函数的过程和还原条件等. 二 理

随时间或空间连续变化的物理量转换为数字信号后，不再是随时间或空间连续变化的物理量，而是一系列离散分布的抽样值阵列.如果一个物理量可以用函数g(x，y)表示，那么该物理量的一系列离散分布的抽样需要满足什么条件才能重构原函数g(x，y)呢 这个答案最早由Whittaker给出，Shannon又将它用于信息论研究，即所谓的惠特克 香农抽样定理.图1 1是连续函数的抽样及重构过程示意图.

图1 1连续函数的抽样及重构

Whittaker和Shannon认为，如果抽样点取得彼此非常靠近，就可以认为这些抽样数据是原函数的表示.对于带限函数只要抽样点之间的间隔不大于某个上限，就可以地重建原函数.

所谓带限函数(band limite functio )，是指这类函数的傅里叶变换只在频率空间的有限区域R上不为零，图1 2给出了一个典型的带限函数在空域及频域的分布，抽样定理适用于带限函数.

图1 2带限函数在空域及频域的分布

考虑二维函数g(x y)在矩形格点上的抽样，抽样函数gs(x y)定义为

gs(x，y)=combxX，yYg(x，y)(1 1

式中comb为梳状函数.抽样函数由δ函数阵列给出，各个δ函数在x方向和y方向上的间隔分别为X和Y.二维函数的抽样过程见图1 3.

图1 3二维函数的抽样过程

设x和y方向上的频率坐标为u和v，则gs(x y)的频谱Gs(u v)可以从函数comb(x/X，y/Y)的变换式与函数g(x，y)的变换式的卷积给出.因为对式(1 1)两边同时作傅里叶变换有

图1 4函数及其频谱

假如函数g(x y)是带限函数，它的频谱Gs只在频率空间(u v)的有限区域R上不为零 抽样函数的频谱不为零的区域可由在频率平面内的每一个(n/X，m/Y)点的周围划出区域R来得到.如果X和Y足够小，则1/X和1/Y的间隔就会足够大，以保障相邻的区域不会重叠.相关函数及其频谱如图1 4所示.

为了确定抽样点之间的较大允许间隔，令2BX和2BY分别表示围住区域R的最小矩形沿u方向和v方向上的宽度，如果抽样点阵的间隔满足

X≤12BX，Y≤12BY(1 6

可保障频谱区域分开而不混频，原函数可恢复，1/(2BX)和1/(2BY)为抽样点阵在u方向和v方向上允许的较大间隔.

用频域中宽度2BX和2BY的位于原点的矩形函数作为滤波函数

H(u，v)=rectu2BXrectv2BY(1 7

让抽样函数gs(x y)的频谱Gs(u v)通过滤波器便能地复原G(u，v)，如图1 5所示.

G(u，v)=Gs(u，v)H(u，v)(1 8

图1 5从抽样函数的频谱恢复原函数的频谱

将G(u，v)=Gs(u，v)H(u，v)对应到空域，有

g(x，y)=gs(x，y)h(x，y)(1 9

式中

gs(x，y)=combxXcombyYg(x，y

=XY∑∞n=-∞∑∞m=-∞g(nX，mY)δ(x-nX)δ(y-mY)(1 10

而且(其中"F-1{}"或"iF{}"表示作逆傅里叶变换

h(x，y)=F-1rectu2BXrectv2BY=4BXBYsinc(2BXx)sinc(2BYy)(1 11

因此

g(x，y)=4BXBYXY∑∞n=-∞∑∞m=-∞g(nX，mY)sinc[2BX(x-nX)]sinc[2BY(y-mY)](1 12

当抽样间隔取较大允许抽样间隔，即

X=12BX，Y=12BY(1 13

最终得到

g(x，y)=∑∞n=-∞∑∞m=-∞g[n/(2BX)，m/(2BY)]sinc{2BX[x-n/(2BX)]}

sinc{2BY[y-m/(2BY)]}(1 14

上式的结果称为Whittaker Shannon抽样定理.它表明，对带限函数在一个间隔合适的矩形阵列上的抽样值，在每一个抽样点上插入一个由sinc函数的乘积构成的插值函数，其权重为相应点上g(x，y)的抽样值，就可以地复原原函数.

图1 6用sinc函数和抽样函数重构原函数

如果带限函数g(x，y 沿u方向和v方向上的带宽分别为2BX和2BY，但抽样时间隔X和Y过大不满足式(1 6)，即X>1/(2BX)，Y>1/(2BY)，称为欠采样(under sampling).此时1/X和1/Y的间隔不够大，频域中函数频谱不为零的区域在每一个(n/X，m/Y)点的周围将会相互重叠，于是，假设仍然用矩形函数rect(u/U)rect(v/V)作为滤波函数，如果取U≥2BX，V≥2BY，将出现混频，反之，如果取U<2BX，V<2BY，将容易出现高频丢失，两者都不能正确复原原函数.图1 6是按抽样定理用sinc函数和抽样函数重构原函数的示意图. 三 容

1 利用Matlab中自带的peaks函数创建一个二维带限函数，通过傅里叶变换观察其频谱，并测量其带宽，理解"带限"的含义;

2 构建二维梳状函数，并显示其空间分布及频谱，观察改变梳状函数的空间间隔——抽样间隔后频谱的变化;

3 利用梳状函数对连续函数抽样，得到该函数的抽样函数，在空域观察抽样函数;

4 观察抽样函数的频谱，并与原连续函数的频谱作比较，体会抽样函数的频谱?梳状函数的频谱，以及连续函数的频谱之间的卷积关系; (5 改变抽样间隔，或调整原连续函数的带宽，观察抽样函数频谱的混叠和分离现象，总结其规律;

6 根据抽样间隔构建二维矩函数滤波器，并对抽样函数的频谱完成滤波和逆傅里叶变换，观察原连续函数的带宽改变，或抽样间隔改变后，利用抽样函数重构原函数的效果，体会欠采样，继而理解抽样定理. 四 参考程序及实验结果

程序流程图

参考程序](CX1_1

1 fxy=cos(peaks(256).2+pi)+1;%构建连续带限函数

2 [rr，cc]=size(fxy);%计算连续函数的大小

3 figure，imshow(fxy，[])%显示连续函数

4 F=fftshift(fft2(fxy));%计算连续函数的频谱

5 figure，plot(abs(F(round(rr/2)+1，:)))，%观察带宽

6 figure，plot(abs(F(:，round(cc/2)+1)))，%观察带宽

7 figure，surfl(abs(F))，shadin interp，colormap(gray);%频谱3D图

8 combxy=zeros(rr，cc);%开始生成comb函数

9 X=4;Y=4;%抽样间隔

10 fo n=1:Y:rr

11 fo m=1:X:cc

12 combxy(n，m)=1;

13 end

14 end

15 figure，imshow(combxy，[]);%显示comb函数

16 C=fftshift(fft2(combxy));%计算comb函数的频谱

17 figure，surfl(abs(C))，shadin interp，colormap(gray);%频谱3D图

18 gxy=zeros(rr，cc);%开始生成抽样函数

19 gxy=fxy.combxy;%生成抽样函数

20 figure，imshow(gxy，[]);%显示抽样函数

21 Gs=fftshift(fft2(gxy));%计算抽样函数的频率

22 figure，surfl(abs(Gs))，shadin interp，colormap(gray);%频谱3D图

23 figure，plot(abs(Gs(:，cc/2+1)))，%观察频谱是否有重叠

24 By=round(rr/2/Y);Bx=round(cc/2/X);%二维矩函数滤波器的宽度

25 H=zeros(rr，cc);%开始生成二维矩函数滤波器

26 H(round(rr/2)+1-By:round(rr/2)+1+By-1，round(cc/2)+1-Bx:round(cc/2)+1+Bx-1)=1;

27 figure，imshow(H，[])%显示二维矩函数滤波器

28 Gsyp=H.Gs;%滤波计算原函数频谱

29 figure，surfl(abs(Gsyp))，shadin interp，colormap(gray);

30 gxyyp=XY.abs(ifft2(Gsyp));%逆傅里叶变换计算原函数

31 figure，imshow(gxyyp，[])%显示还原的原函数

实验结果]

图1 7带限函数为cos(peaks(256)2+π)+1，X=Y=4个像素的实验结果

在上面的实验中，从G(u，v)|v=0和G(u，v)|u=0两条剖线图可以看到:所选带限函数cos[peaks(256)2+π]+1的带宽2BX<64个像素，而2BY≈64个像素.因为图像大小为256256像素，根据快速傅里叶变换的性质，傅里叶变换后较高空间频率为256/2=12 m-1，频域中相邻两个像素点的频率差为128/128= m-1，换言之，所选带限函数的带宽2BX<6 m-1，而2BY≈6 m-1.按照抽样定理，能够重构原函数的条件是抽样间隔至少满足X=1/(2BX)，Y=1/(2BY)，所以至少选择抽样间隔Y=1/(2BY)=256/64=4个像素，为了方便，取X=Y=4个像素.

从Gs(u v)的三维分布图以及Gs(u v)|u=0剖线图可以看到，频域中函数频谱不为零的区域在每一个(n/X，m/Y)点的周围没有相互重叠，所以用频域中宽度为2BX和2BY的位于原点的矩形函数作为滤波函数，可以直接得到原函数的频谱，再作逆傅里叶变换即得到基本没有失真的原带限函数.

实验中傅里叶变换用离散傅里叶变换(快速傅里叶变换FFT)完成，若空域抽样间隔分别为X和Y，则离散信号频谱的幅度应该是原函数频谱的1/(XY)，为保障二维FFT计算幅度正确，必须让用FFT计算得到的结果乘以XY，为此实验程序的第30行乘了XY.

如果加宽带限函数的带宽(方法是将所选带限函数cos[peaks(256)2+π]+1中的2改为4，或更大的数值)，同时保持抽样间隔X=Y=4个像素不变，重复前面的运算，结果如下:

图1 8带限函数为cos[peaks(256)4+π]+1，X=Y=4个像素的实验结果

此时从剖线图 (u v |u=0可以看到带限函数cos[peaks(256)4+π]+1的带宽2BY≈128个像素，按照抽样定理，能够重构原函数的条件是抽样间隔至少满足Y=1/(2BY)，所以至少选择抽样间隔Y=1/(2BY)=256/128=2个像素，但实验中仍然取X=Y=4个像素，不满足抽样定理，应该发生欠采样.

从Gs(u v)的三维分布图以及Gs(u v |u=0剖线图可以看到，频域中函数频谱不为零的区域在每一个(n/X，m/Y)点的周围明显相互重叠，所以用频域中宽度为U=1/X=64和V=1/Y=64个像素的位于原点的矩形函数作为滤波函数，得到原函数的频谱有明显的高频丢失，作逆傅里叶变换复原的原带限函数当然也有明显失真.

如果仍选带限函数为cos[peaks(256)4+π]+1，抽样间隔也仍取X=Y=4个像素，但频域中用宽度为2B =128和2BY=128个像素的位于原点的矩形函数作为滤波函数，可以看到得到的原函数频谱有明显的混频，复原的原带限函数有严重失真，结果如图1 9所示.

图1 9实验结果.带限函数为cos[peaks(256)4+π]+1，X=Y=4个像素，频域中用宽度为2BX=128和2BY=128个像素的矩形函数"