matlab实现FFT（快速傅里叶变换）

任务

接上篇，尝试使用快速傅里叶变换实现空域到频域转换

过程

判定2的幂次

由算法可知，对于长度为2的幂次的序列可以进行变换，但对于大多数情况非2的幂次应该如何解决？答案是补零至2的幂次。

判定

简单的判定方法，利用2的幂次（M）在二进制中的独特性，首位为1，余位为0；且减1后（M-1）原首位为0，余位为1，则有M&M-1=0

matlab中的按位与函数为bitand()，与&不同

 bitshift(num,n); num表示需要移位的数,十进制; n为正左移,负右移。

求大于一个数字的二进制的最高位

思路，通过移位，当非2次幂数的第i位移出使原数为0，则最高位为i。可以将2次幂数-1转换为相同情况。

函数实现

首先强制类型转换，转换为固定的uint8类型。

matlab没有三目运算符，则用if…else语句实现。将特殊情况0单独列出，

%% 判定int类型数M是否为2的n次幂，
% 若是，state=1，返回n
% 若否，state=0，返回最临近的n次幂
function [state,n]=is2power(M)
% 强制类型转换
M=uint8(M);
if M==0
    % M为0的特殊情况
    state=0;
    n=0;
elseif bitand(M,M-1)==0
    state=1;
    n=highestOneSite(M-1);
else
    state=0;
    n=highestOneSite(M);
end
end

调用临近最高位的求取函数highestOneSite(M)

%% 返回大于等于某数M的2的n次幂的n
function [n]=highestOneSite(M)
% 强制类型转换
M=uint8(M);
for i=0:8
    if bitshift(M,-i)==0
        n=i;
        break
    end
end
end

位的对换

算法实现需要将输入数据排列为满足算法需要的次序，可以推导发现，满足的映射规则为第x个数为x的二进制数颠倒后的结果。

a是0~255的整数

b=dec2bin(a,8);%b是8位’0′,’1’字符串

h=b(1:4); %高四位

l=b(5:8); %低四位

MATLAB中怎么将一个常数（0到255）转为8位二进制数–dukinkin 

但在matlab中，由于数组下标从1开始，需要进行-1的操作进行映射

即F(abcd）《-》f(dcba)

bottom-up实现

足足思索了几个小时，终于捋清x点变换的逻辑。

[latexpage]
\begin{equation}
F(u)=\sum_{x=1}^{N}f(x)W_{N}^{ux}
\end{equation}

指数部分函数W表示

其中W可视为函数：

[latexpage]
\begin{equation}
W_{N}^{ux}=\exp [-j 2 \pi ux/ N]
\end{equation}

matlab代码为：

%% 指数部分
function W=W(N,ux)
W=exp(-1i*2*pi*ux/N);
end

首先从一维开始理解，由于需要变换的矩阵长度已经过判定2的幂次并补零至最邻近的2的幂次（N=2^n)，则可通过不断奇偶二分，从原始的N点变换，分解为两个N/2点变换，再到N/4点变换、……、8点变换、4点变换、2点变换。最后的1点变换，因为公式（1）中u为0，则指数项为0，即1点变换视为原值。

上面的思路确定了如何从需要变换的函数得到分解因子，本次实现通过从因子逐层自底向上推出。重点在于清楚公式中的分子的表示，以及在存储矩阵中的位置。

function T=fftL1(Xn)
%% 得出共有多少点
[r,~]=size(Xn);
% 找到最近的2次幂N=2^n
[~,n]=is2power(r);
N=2^n;
% 拓宽至2次幂
Xn_2=zeros(N,1);
Xn_2(1:r,1)=Xn;
% 运行存储空间
T=zeros(N,1);
%1点变换(对换位置)
for i=1:N
    T(i)=Xn_2(reverseBit(i-1,n)+1);
end
%从2^point点变换依次开始
for point=1:n
    %组数
    groups=N/2^point;
    %组员数
    members=2^point;
    for group=1:groups
        %对子数
        pairs=2^(point-1);
        for pair=1:pairs
            site=(group-1)*members+1;
            M=pairs;
            u=pair-1;
            even=T(site+u);
            odd=T(site+u+M)*W(2*M,u);
            T(site+u)=1/2*(even+odd);
            T(site+u+M)=1/2*(even-odd);
        end
    end
end
end

二维实现

利用matlab实现DFT(离散傅里叶变换）与优化 的读取图像和灰度图转换，得到了灰度图像grayImg

tic
%代码块
%灰度图大小
[M,N]=size(grayImg);
[state,n]=is2power(M);
R=2^n
[state,n]=is2power(N)
C=2^n
%扩展图
newImg=zeros(R,C);
newImg(1:M,1:N)=grayImg;
%保存一维Fourier变换的结果
oneD=zeros(R,C);
%保存二维Fourier变换的结果
twoD=zeros(R,C);
%% 算法开始
for row =1:R
    oneD(row,:)=(fftL1((newImg(row,:)).')).';
end
for col =1:C
    twoD(:,col)=fftL1(oneD(:,col));
end
toc

再利用之前的频移函数showShiftImg(twoD)得到频移到中心点的图像。

matlab技巧——如何利用其他文件夹下的函数

addpath D:\Matlab
即可使用路径D:\Matlab下的函数

结果展示

ans =
442 785 3
时间已过 17.532276 秒。

现在运行给定图像需要的时间实际为17.5s，相比于优化后的DFT的12s居然不降反升，所以要想达到更有效的应用，需要继续优化。

参考文章

• 《图像工程 （上）<图像处理和分析>》 章毓晋
• matlab中调用其他文件夹中的方法