Fourier_01_图像傅里叶变换的公式推导（一）SP——指正与补充
posted on 14 Feb 2018 under category 2018_02
图像傅里叶变换的公式推导SP——指正与补充
1 Point_out_Q1 fftshift(fft(fftshift(function)))的形式
昨天在下更新一波以后，LWL大佬指出，在fftshift变换的问题上面，国内和国外使用方式有些不一致。当时我的第一反应是在周期序列的背景下，shift的运算规整方式不同不影响操作和计算，为了确定fftshift(fft(function))和fftshift(fft(fftshift(function)))两种写法究竟会带来什么样的区别，在下设计了一段代码来进行验证：
clear all
clc
%%%---LWL大佬的指正与补充---%%%
%%%---读取归一化lena图---%%%
I = imread('lena.jpg');
I = rgb2gray(I);
I = im2double(I);
figure;
%%%---正变换---%%%
%%%---fftshift(fft(function))形式---%%%
I_fft = fftshift(fft2((I)));
I_fft_d = (log(abs(real(I_fft)+1)));
subplot(221);imshow(I_fft_d, [ ])
title('I-fft-d')
%%%---fftshift(fft(fftshift（function)))形式---%%%
I_fft2 = fftshift(fft2(fftshift(I)));
I_fft_d2 = (log(abs(real(I_fft2)+1)));
subplot(222);imshow(I_fft_d2, [ ])
title('I-fft-d2')
%%%---逆变换---%%%
I_ifft = (ifft2(ifftshift(I_fft)));
subplot(223);imshow(I_ifft, [ ])
title('I-ifft')
I_ifft2 = ifftshift(ifft2(ifftshift(I_fft2)));
subplot(224);imshow(I_ifft2, [ ])
title('I-ifft2')
%%%---校验---%%%
sub_1 = I - I_ifft;
sub_2 = I - I_ifft2;%% 两个方式变换回去的结果都是和原信息一致
sub_3 = I_fft - I_fft2;%% 相角相反
sub_4 = I_ifft - I_ifft2;%% 结果一致
在这段代码中，I_fft系列变量采用fftshift(fft(function))的方式，I_fft2系列变量采用fftshift(fft(fftshift(function)))的方式，正反变换之后使用sub_1、sub_2以及sub_4验证，原图和两种变化得到的图结果都是一样的（计算误差范围内）。运行结果如下：
而对于频谱的部分，sub_3的结果并不是足够小，反而还挺大的。不着急，找这类问题最好的解决方式还是回到数据上来分析。
既然是有较大值的点，那么不妨从最大值部分来找。在命令行输入：
[a, b] = max(max(sub_3));
[c, d] = max(max(sub_3’));
找到最大值点的位置(d, b)，然后回到I_fft和I_fft2的数据上，观察到如下图：
发现两者在区域内的数字，虚部和实部数值上相等但是符号相反。初步得到一个推断（如代码注释）：这两种变换初始相角相差π。（相角这个部分准备在物理意义里面说的，提前在这里提及应该可以建立一点点小印象）
最终得到的结论是：这两种变换方式不会改变的实质，在操作上和计算上确实没有问题，在关心频谱实部的情况下，两者得到的频域显示结果是一样的。
2 Point_out_Q1 离散傅里叶变换矩阵
昨天的文章结尾处提到计算速度的问题，WY大佬留言提到“把循环改成矩阵，内存换速度”的建议。对此在下下去查证了一下离散傅里叶变换矩阵的相关资料，老实说，现在还没消化透彻，借助MATLAB已有的函数勉强演示一下离散傅里叶变换矩阵的使用，大家如果感兴趣，请参照离散傅里叶变换的矩阵表示及其运算量这个课件，以及MATLAB的help文档中dftmtx函数的使用。
简单的描述一下构造离散傅里叶变换矩阵的过程。在昨天提到的二维离散傅里叶变换的式子中，可以观察到exp(-2pi*(function))的式子是有规律的，对于这种规律可以总结成如下公式：
对于式子WN我们可以构建出矩阵形式和空频的转换关系如下：
基本的公式原理得到以后，要设计实验验证离散傅里叶变换矩阵的正反变换效果，以及掌握基本的函数操作方法。
在下对于这个问题，使用dftmtx制造了一个离散傅里叶变换矩阵H_ft，以及逆变换矩阵H_ift进行变换然后和fft的变换方式对比结果，代码如下：
clear all
clc
%%%---WY大佬的指正与补充---%%%
%%%---设计冲激函数---%%%
n = 256;
x = 1:1:n;
% J = 0.5*sin(x/n*2*pi)+0.5;%% 开始准备用正弦函数，但是频谱特征一般
J = (x == 128);%% 单位冲激函数
figure;
subplot(331);plot(J);
title('冲激函数J')
%%%---通过离散傅里叶变换矩阵进行变换---%%%
H_ft = dftmtx(n);
J_ft = fftshift(H_ft*J');
J_ft_d = log(abs(real(J_ft))+1);%% _d标志表示display，用于显示的结果，不计算
subplot(332);plot(1:1:256, J_ft_d)
title('J-ft-d')
%%%---逆变换---%%%
H_ift = conj(dftmtx(n))/n;
J_ift = H_ift*ifftshift(J_ft);
subplot(333);plot(1:1:256, abs(J_ift)')
title('abs(J-ift)')
%%%---fft进行傅里叶变换---%%%
J_fft = fftshift(fft(J));
J_fft_d = log(abs(real(J_fft))+1);
subplot(335);plot(1:1:256, J_fft_d)
title('J-fft-d')
%%%---逆变换---%%%
J_ifft = ifft2(fftshift(J_fft));
subplot(336);plot(1:1:256, abs(J_ifft))
title('abs(J-ifft)')
%%%---校验参数---%%%
sub_5 = J - J_ift';
sub_6 = J - J_ifft;
sub_7 = J_ift' - J_ifft;
sub_8 = J_fft - J_ft';
subplot(337);plot3(1:1:256, real(J_fft), imag(J_fft));title('J-fft')
subplot(338);plot3(1:1:256, real(J_ft'), imag(J_ft'));title('J-ft')
subplot(339);plot3(1:1:256, real(sub_8), imag(sub_8));title('J-fft - J-ft')
结果如图：
这个结果图片和sub的计算结果表明dftmtx和fft运算结果的一致。需要提醒大家的是注意矩阵运算的行和列，有几个地方需要进行转置运算，否则会报出“Index exceeds matrix dimensions.”的报错。
这里放出结果图主要是想展示一下两种方式得到的频谱图的差别。在做sub_8运算的时候，发现在这次的变换中也出现了和今天SP验证中一样的问题，就是实部相等，虚部不一致（如图J-fft – J-ft）。
因为这一次是一维函数，用x轴和y轴来表示real部分和imag部分，可以表示出一维情况下，两种不同变换方式的复数频谱结果，观察到漂亮的螺旋结构（好吧，在下的审美只能这样了），两个螺旋看上去方向相反，差值结果具有变化规律。这样子，前面的相角差值的推断得到了进一步的验证，并且在这个变换的对比中，正好也出现了“相角相差π”的情况。
回到我们要解决的主要问题：离散傅里叶变换矩阵的变换效果，在这种操作过程下和快速傅里叶变换得到的结果一致，这种“内存换时间”的做法，非常值得使用。当然这种方法正是快速傅里叶变换的核心思想，蝶形的运算模式剩下了大量的计算量。并且，因为这个矩阵和快速傅里叶变换极强的关联性，dftmtx(n)得到的结果其实就是调用fft(eyes(n))（在MATLAB的源码里面就是这样调用的）。在下今天只吃透了这么多，具体的数学过程请看本章开始的参考资料。
非常感谢LWL和WY两位大佬的指正和补充！也希望其他走过路过的盆友们多多指教！明天（不会）更新更基础的卷积运算部分，下一个chapter见~