二维抽样定理

fxy=cos(peaks(256).*2+pi)+1;    % 构建连续带限函数
[rr,cc]=size(fxy);              % 计算连续函数大小
figure(1);imshow(fxy,[]);title("二维连续带限函数")         % 显示连续函数
F=fftshift(fft2(fxy));          % 计算连续函数的频谱
% figure(2);plot(abs(F(round(rr/2)+1,:)));    % x方向带宽
% figure(3);plot(abs(F(round(cc/2)+1,:)));    % y方向带宽
figure(2);surfl(abs(F)),shading interp,colormap("gray");title("带限函数频谱图")    % 带限函数频谱3D图
combxy=zeros(rr,cc);            % 生成comb函数
X=4;Y=4;                        % 抽样间隔
n=1:X:rr;
m=1:Y:cc;
combxy(n,m)=1;
figure(3);imshow(combxy,[]);title('comb函数');            % 显示comb函数
C=fftshift(fft2(combxy));       % 计算comb函数的频谱
figure(4);surfl(abs(C)),shading interp,colormap("gray");title("comb函数频谱图")  % comb函数频谱3D图
gxy=fxy.*combxy;                % 生成抽样函数
figure(5);imshow(gxy,[]);title("抽样函数");               % 显示抽样函数
Gs=fftshift(fft2(gxy));         % 抽样函数频谱
figure(6);surfl(abs(Gs)),shading interp,colormap("gray");title("抽样函数频谱图")  % comb函数频谱3D图
figure(7);plot(abs(Gs(:,cc/2+1)));  % 观察频谱是否有重叠
By=round(rr/2/Y);Bx=round(cc/2/X);  % 二维矩阵函数滤波器的宽度
H=zeros(rr,cc);
H(round(rr/2)+1-By:round(rr/2)+By,round(rr/2)+1-Bx:round(rr/2)+Bx)=1;
figure(8);imshow(H,[]);             % 显示二维矩阵函数滤波器
Gsyp=H.*Gs;                     % 滤波计算原函数频谱
figure(9);surfl(abs(Gsyp));shading interp,colormap("gray");title("滤波后原函数频谱");
gxyyp=X*Y.*abs(ifft2(Gsyp));    % 逆傅里叶变换计算原函数
figure(10);imshow(gxyyp,[]);title('还原的原函数');% 显示还原的原函数