很多人想找找不到的2维DFT矩阵算法和手动Fspecial退化函数实现

今天起床脑仁疼，趁着脑袋不清醒把这个月的工作总结一下haha，提前祝csdn和各位小老板双11快乐哈。

基于课业要求或者想了解matlab内库的童鞋想自己写DFT的和退化模型的可以看下。

咱们只谈干货，二维DFT和退化模型的基础知识什么的就不在这说了，可以去看别的博主。

---

如果尝试去写的话遇到的问题可能有：①使用四个for循环的2维DFT运行时间比较长，我用实验室的i5试了下大概100分钟的样子。考虑到matlab的老本行，还是用矩阵写。②矩阵的话"市面上"只有一维的，对于二维的其实不难猜到它的形式其实就是一维的推广。对一个MXN的矩阵进行运算得到一个MXN的可能会猜到就是MMXMNXNN，我写了两张A2的草稿纸证明了下是完全正确的。③退化函数Fspecial可能有人用的贼6，但是其内部机理还是有点说头的。好像如果只是按照给的H函数去做不会做出斜向的退化效果，这个时候就要湿考辽。

---

%%
%二维DFT
clc;
clear;
Q=imread('im.bmp');
I1=rgb2gray(Q);                                  %读取后转为灰度图
I=double(I1);                                    
[M,N] = size(I);
%构造两个矩阵分别为MxM，NxN，在I前后分别相乘
m=[1:M];         
n=[1:N]                                      
subplot(221);imshow(I1);title('时域图像');
WM=exp(-j*2*pi/M);
mm=m'*m;
WMmm=WM.^mm;
WN=exp(-j*2*pi/N);
nn=n'*n;
WNnn=WN.^nn;                                     %是nn里的元素分别乘方，指数为WN
F=WMmm*I*WNnn;                                   %F里的元素均为复数
FRE=abs(F);                                    %取模并缩小范围
FRE=log(FRE+1);
%min=min(FRE(:));  
%max=max(FRE(:));
FRE_=im2double(FRE);%转换为unit8便于显示
FRE_=FRE_/25;
y=fftshift(FRE_);
subplot(222);imshow(y);title('频域DFT图像');
%imwrite(y,'dft队员.bmp');

使用冈萨雷斯书里经典的配图，效果如下：

可能有一些亮度上的缺陷，自行调整就好了。

clc,clear;
[I,map]=imread('im.bmp'); % 原始图像X数据类型为uint8
f=double(rgb2gray(I));
[M,N]=size(f);
% imshow(X,[]);title('ImageOrg.bmp');
H = zeros(M,N);
T=1;a=0.1;b=0.1;
for i = 1:M
   for k = 1:N
       ii=i-M/2+0.001;
       kk=k-N/2+0.001;
       H(i,k) = (T/(pi*(ii*a+kk*b))) * sin(pi*(ii*a+kk*b))*exp(-j*pi*(ii*a+kk*b));
   end
end
F=fft2(f);
F=fftshift(F);
G=F.*H;
g=ifftshift(G);
gg=ifft2(g);
ggg=abs(gg)/200;
%gg=abs(g)/3;
imshow(ggg);title('result.bmp');

先解释下突破点，就是H函数其实原点在中心，而不是简单地把它带进去算，得到矩阵再相乘。然后有个小问题是会有ii=kk=0，这样分母就会为0，得到NAN，因此在后面加个很小的数（比如0.001）就行了。

退化没有用自己的DFT，变换部分直接用的库函数，效果如下：

用fspecial函数的效果如下，自己写的这个就算是成功了罢：