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居然不降反升,所以要想达到更有效的应用,需要继续优化。

参考文章

4 comments

Leave a comment

Your email address will not be published. Required fields are marked *