计算图像傅立叶变换的过程很简单：首先对每一行做一维FFT,然后对每一列做一维FFT。具体来说,先对第0行的N个点做FFT(实部有值,虚部为0),将FFT输出的实部放回原来第0行的实部,FFT输出的虚部放回第0行的虚部,这样计算完全部行之后,图像的实部和虚部包含的是中间数据,然后用相同的办法进行列方向上的FFT变换,这样N*N的图像经过FFT得到一个N*N的频谱。

I = imread('peppers.png'); % 读入图像
I = I(10+[1:256],222+[1:256],:); % 图像剪切
figure, imshow(I); % 显示原始图像
title('Original Image'); % 设置图像标题
LEN = 31; % 设置长度
THETA = 11; % 设置角度
PSF = fspecial('motion',LEN,THETA); % 生成滤波器
blurred = imfilter(I,PSF,'circular','conv'); % 图像卷积计算
figure, % 新建窗口
imshow(blurred); % 显示处理后的图像
title('Blurred'); % 设置图像标题
wnr1 = deconvwnr(blurred,PSF); % 使用真实PSF进行图像复原
figure, % 新建窗口
imshow(wnr1); % 显示处理后的图像
title('Restored, True PSF'); % 设置图像标题
wnr2 = deconvwnr(blurred,fspecial('motion',2*LEN,THETA));
% 模糊运动长度取值是真实长度的2倍，2*LEN
figure, % 新建窗口
imshow(wnr2); % 显示处理后的图像
title('Restored, "Long" PSF'); % 设置图像标题
wnr3 = deconvwnr(blurred,fspecial('motion',LEN,2*THETA));
% 模糊运动角度取值是真实角度的2倍，2*THETA
figure, % 新建窗口
imshow(wnr3); % 显示处理后的图像
title('Restored, Steep'); % 设置图像标题
noise = 0.1*randn(size(I)); % 生成噪声信号
blurredNoisy = imadd(blurred,im2uint8(noise)); % 加入图像
figure, % 新建窗口
imshow(blurredNoisy); % 显示处理后的图像
title('Blurred & Noisy'); % 设置图像标题
wnr4 = deconvwnr(blurredNoisy,PSF); % 图像恢复
figure, % 新建窗口
imshow(wnr4); % 显示处理后的图像
title('Inverse Filtering of Noisy Data'); % 设置图像标题
NSR = sum(noise(:).^2)/sum(im2double(I(:)).^2); % 计算信噪比
wnr5 = deconvwnr(blurredNoisy,PSF,NSR); % 图像恢复
figure, % 新建窗口
imshow(wnr5); % 显示处理后的图像
title('Restored with NSR'); % 设置图像标题
wnr6 = deconvwnr(blurredNoisy,PSF,NSR/2); % 信噪比降低一半，NSR/2
figure, % 新建窗口
imshow(wnr6); % 显示处理后的图像
title('Restored with NSR/2'); % 设置图像标题
NP = abs(fftn(noise)).^2; % 傅立叶变换
NP