matlab实现DFT(离散傅里叶变换）与优化

任务

《数字图像处理》课程布置：

编写程序(建议Matlab )对某选定图像(自行转换为灰度图)展开傅里叶变换,提取傅里叶变换图像(将频率原点移至图像中心) , 并形成实验报告。

为与matlab中代码一致，便于理解，此处的公式我全部整理为数组从下标1开始。此外，增强特殊性，图像大小为NxM

基本知识

此前尝试用matlab实现傅里叶变换，用基本公式对方波函数进行傅氏变换。但其特点：

函数连续（在matlab进行抽样离散处理，绘制波形）
范围趋近无限（在matlab中把n取大值范围）

可见，对于模拟值用计算机进行傅氏变换需要首先进行离散处理，那么对于本身是离散值的数字图像，则直接利用离散傅里叶变换即可。

[latexpage]
\begin{equation}
F(u,v)=(\frac{1}{M}\cdot\frac{1}{N})\sum_{x=1}^{M}\sum_{y=1}^{N}f(x,y)e^{-i2\pi(\frac{ux}{M}+\frac{vy}{N})}
\end{equation}

这里图像大小为NxM
前面括号内的系数其实可舍去，主要起到归一化的作用，因此形式不必纠结。

首先，将二重变换利用傅氏变换的性质可以分解为两次单重变换：

[latexpage]
\begin{equation}
F(x,v)=(\frac{1}{N})\sum_{y=1}^{N}f(x,y)e^{-i2\pi\frac{vy}{N}}
\end{equation}
\begin{equation}
F(u,v)=(\frac{1}{M})\sum_{x=1}^{M}F(x,v)e^{-i2\pi\frac{ux}{M}}
\end{equation}

与公式相一致：本人代码先按行遍历，对列（y表示）变换；再按列遍历，对行（x表示）变换

过程

读取图像

%读取图像
imgFileName='{{文件路径}}\lena.jpg';
Img=imread(imgFileName);
%图像大小
size(Img)
%图像信息
imgInfo=imfinfo(imgFileName);

测试图片选取经典上图，但大小是压缩的200×200，因此对于初次使用比较友好。
matlab中句末引号可选，不加则会输出左值结果至命令行窗口，如上述代码的size(Img)

转换为灰度图

%% 转换为灰度图
%如果图像是rgb3通道,进行灰度处理
if(imgInfo.('NumberOfSamples')==3)
    grayImg=rgb2gray(Img);
    figure(1)
    imshow(grayImg)
%图像是单一通道,直接使用
elseif(imgInfo.('NumberOfSamples')==1)
    grayImg=(Img);
    figure(1)
    imshow(grayImg)
end

通过在读取的图像信息判断是否为单通道，是则直接使用，否则将RGB图像通过转换函数rgb2gray()转换为灰度图。

三重循环实现

tic
%% 三重循环实现
%代码块
%灰度图大小
[M,N]=size(grayImg);
%保存一维Fourier变换的结果
oneD=zeros(size(grayImg));
%保存二维Fourier变换的结果
twoD=zeros(size(grayImg));
%% 先按行
for row=1:M
    for v=1:N
        for col=1:N
            oneD(row,v)=oneD(row,v) + double(grayImg(row,col)) * exp(-1i*2*pi*(v-1)*col/N);
        end
    end
end
%% 再按列
for col=1:N
    for u=1:M
        for row=1:M
            twoD(u,col) = twoD(u,col) + double(oneD(row,col)) * exp(-1i*2*pi*(u-1)*row/M) ;
        end
    end
end
toc

可以通过tic，toc放在代码前后进行时长测试
注意u和v的范围是从0开始取值

ans =
200 200 3
时间已过 12.674382 秒。

频移函数

function [twoD]= showShiftImg(twoD)
arg=abs(twoD);
[M,N]=size(twoD);
% F=fftshift(arg2D);
% 核心步骤
F=[arg2D(M/2+1:M,N/2+1:N),arg2D(M/2+1:M,1:N/2);
     arg2D(1:M/2,N/2+1:N),arg2D(1:M/2,1:N/2)];
T=log(F+1);
imshow(T,[]);

这里直接显示经过处理后的傅氏变换的图像
核心是将原图四个角点移至中心位置，可以通过matlab中矩阵处理的语法便捷实现
由于非偶数会出现不能被2整除情况，因此需要做取整处理

[latexpage]
图块分为四块，如下所示进行位置调换
\begin{equation}
\begin{bmatrix}
1 &2 \\
3 & 4
\end{bmatrix}\Rightarrow \begin{bmatrix}
4 & 3\\
2 & 1
\end{bmatrix}
\end{equation}

优化

使用上述算法经过讨论，复杂度为O(N^3)，并且在处理更大尺寸的图像（老师提供785×442）时，运行时间极长。

因为了解到matlab对矩阵运算具有高效快速优点，便把思考放在如何将算法通过矩阵运算实现。

思路同样是通过两个一维变换完成，因此只以分析按行遍历，对列(y)展开为例，对思路进行陈述。

观察到公式中的基函数部分其实已经固定，因此直接通过相乘获得某个频率分量（Vn）对应的不同y值得出（某个频率分量）所有计算需要的指数部分，如下所示：

[latexpage]
\begin{equation}
\underbrace{
\left.
Vk*Yn= \begin{bmatrix}
VkY1 & VkY2 & \cdots & VkYN \\
VkY1 & VkY2 & \cdots & VkYN \\
\vdots & \vdots & \ddots & \vdots \\
VkY1 & VkY2 & \cdots & VkN
\end{bmatrix}
\right \}}M
\end{equation}
\begin{equation}
\end{equation}

共N列

下一步将图像信息与对应的分量对应相乘，得到，那么加和则是对每行相加——每一行的x保持不变。最终得到一维转换后的结果：

    %  [
    %     F(x1,v1),F(x1,v2),...,F(x1,vc);
    %     F(x2,v1),F(x2,v2),...,F(x2,vc);
    %     ....
    %     F(xr,v1),F(xr,v2),...,F(xr,vc);
    %   ]

具体代码

tic
%代码块
%灰度图大小
[M,N]=size(grayImg);
%% 先按行
for v=0:N-1
    % 行向量
    Yn=1:N;
    % 列向量
    Vn=(v*ones(1,M)).';
    % 得出指数部分
    ExpY=Vn*Yn;
    % 对应相乘
    corTime=double(grayImg).*exp(-1i*2*pi*ExpY/N);
    % 按行相加
    oneD(:,v+1)=sum(corTime,2);
end
showShiftImg(oneD,1)
%% 再按列
for u=0:M-1
    % 列向量
    Xn=[1:M].';
    % 行向量
    Un=(u*ones(1,N));
    % 得出指数部分
    ExpX=Xn*Un;
    % 对应相乘
    corTime=double(oneD).*exp(-1i*2*pi*ExpX/M);
    % 按列相加
    twoD(u+1,:)=sum(corTime);
end
toc

时间已过 16.837070 秒。

此用时与分钟相比下降很多，但比较系统中的fft()，用时约为0.1s，因此下一步可以思考如何进行快速傅里叶变换的使用。