傅丽叶变换

       因为这些都是从《数字图像处理基础 . 阮秋琦》（注：这本书确实写的不错，虽然没有具体的算法实现，并且有些小错误，但是通俗易懂）里搬过来的，但又是进行图像傅丽叶变换之前需要了解的一些知识，所以不敢写原创，只能算是转载吧！如果想直接了解“图像傅丽叶变换”的算法实现请看下一章《傅丽叶变换（二）》

 数字图像处理的方法主要分为两大类：

   一个是空间域处理法(或称空域法)，

   一个是频域法(或称变换域法)。

   在频域法处理中最为关键的预处理便是变换处理。

目前，在图像处理技术中正交变换被广泛地运用于图像特征提取、图像增强、图像复原、图像识别以及图像编码等处理中。

傅里叶变换是大家所熟知的正交变换。在一维信号处理中得到了广泛应用。把这种处理方法推广到图像处理中是很自然的事。下面从傅丽叶的基本概念讲起，再讲图像处理中的傅丽叶变换算法实现。

一维傅里叶变换

定义

傅里叶变换在数学中的定义是严格的。设f(x)为x的函数，如果满足下面的狄里赫莱条件：

    （１）具有有限个间断点；

    （２）具有有限个极值点；

（３）绝对可积。

则有下列二式成立

式中x是时域变量，u为频率变量。

如令， 则有

通常把以上公式称为傅里叶变换对。

函数f(x)的傅里叶变换一般是一个复数，它可以由式（5）表示

或写成指数形式

把|F(x)|叫做f(x)的傅里叶谱，而叫相位谱。

傅丽叶变换的几个概念：

    （１）只要满足狄里赫莱条件，连续函数就可以进行傅里叶变换，实际上这个条件在工程运用中总是可以满足的。

    （２）连续非周期函数的傅里叶谱是连续的非周期函数，连续的周期函数的傅里叶谱是离散的非周期函数。

二维傅丽叶变换

定义

傅里叶变换可推广到二维函数。如果二维函数f(x,y)满足狄里赫莱条件，那么将有下面二维付里哀变换对存在：

与一维傅里叶变换类似，二维傅里叶变换的幅度谱和相位谱如下式

式中：F(u,v)是幅度谱；是相位谱；E(u,v)是能量谱。

傅里叶变换的性质 

傅里叶变换有许多重要性质。这些性质为实际运算处理提供了极大的便利。这里，仅就二维傅里叶变换为例列出其主要的几个性质。 

具有可分性

这个性质说明一个二维傅里叶变换可用二次一维傅里叶变换来实现。 

线性 

傅里叶变换是线性算子，即 

 

共轭对称性 

 如果F(u, v)是f(x, y)的傅里叶变换，F*(-u, -v)是f(-x, -y) 傅里叶变换的共轭函数，那么 

F(u, v) = F*(-u, -v)

旋转性

 如果空间域函数旋转的角度为 ，那么在变换域中此函数的傅里叶变换也旋转同样的角度，即 

比例变换特性 

如果F(u, v) 是f(x, y)的傅里叶变换。a和b分别为两个标量，那么 

帕斯维尔（Parseval）定理

这个性质也可称为能量保持定理。如果F(u, v) 是f(x, y)的傅里叶变换，那么有下式成立 

 这个性质说明变换前后并不损失能量

相关定理 

    如果，f(x),g(x)为两个一维时域函数；f(x,y)和g(x,y)为两个二维空域函数，那么，定义下二式为相关函数

由以上定义可引出傅里叶变换的一个重要性质。这就是相关定理，即

式中F(u, v) 是f(x, y)的傅里叶变换，G(u, v)是g(x, y) 的傅里叶变换，G*(u, v)是G(u, v)的共轭，g*(x, y)是g(x, y)的共轭。

卷积定理 

如果f(x)和g(x)是一维时域函数，f(x,y)和g(x,y)是二维空域函数，那么，定义以下二式为卷积函数，即

由此，可得到傅里叶变换的卷积定理如下

式中 F(u,v) 和 G(u,v)  分别是f(x)和g(x)的傅里叶变换。

---

离散傅里叶变换

离散傅里叶变换使得数学方法与计算机技术建立了联系，这就为傅里叶变换这样一个数学工具在实用中开辟了一条宽阔的道路。因此，它不仅仅有理论价值，而且在某种意义上说它也有了更重要的实用价值。

离散傅里叶变换的定义 

如果x(n)为一数字序列，则其离散傅里叶正变换定义由下式来表示

傅里叶反变换定义由下式来表示

由(1)和(2)式可见，离散傅里叶变换是直接处理离散时间信号的傅里叶变换。如果要对一个连续信号进行计算机数字处理，那么就必须经过离散化处理。这样，对连续信号进行的傅里叶变换的积分过程就会自然地蜕变为求和过程。

快速傅里叶变换(FFT) 

   快速傅里叶变换并不是一种新的变换，它是离散傅里叶变换的一种算法。这种方法是 在分析离散傅里叶变换中的多余运算的基础上，进而消除这些重复工作的思想指导下得到的，所以在运算中大大节省了工作量，达到了快速运算的目的。扩大了傅里 叶变换的使用范围。下面我们从基本定义入手，讨论其原理。

对于一个有限长序列 {x(n)}(0 <= n < N-1) ,它的傅里叶变换由下式表示

令

因此，傅里叶变换对可写成下式

将正变换式（4）展开可得到如下算式

上面的方程式（6）可以用矩阵来表示

从上面的运算显然可以看出，要得到每一个频率分量，需进行N次乘法和N-1次加法运算。要完成整个变换需要N2次乘法和N(N-1)次加法运算。当序列较长时，必然要花费大量的时间。 

观察上述系数矩阵，发现W^mn是以为N周期的，即 

例如，当N=8时，其周期性如图3—6所示。由于 exp(iθ）= cosθ+isinθ,即exp(-iθ）= cos(-θ)+isin(-θ) = cosθ-isinθ,所以

当N=8时，可得：

 图 1,N=8 时W^mn的周期性和对称性 

 快速傅里叶变换简称FFT。算法根据分解的特点一般有两类，一类是按时间分解，一类是按频率分解。下面介绍一下ＦＦＴ的基本形式及运算蝶式流程图。 

把x(n)分成偶数点和奇数点,即:

由此，离散傅丽叶变换可写成如下的形式：

    (10)

这种算法的流程图如图2所示：图(a)输入为顺序的，运算结果是乱序的；图(b)输入为乱序的，运算结果是顺序的。

图2，FFT蝶式运算流程图 (按时间分解)

用计算机实现快速付傅里叶变换 

 利用FFT蝶式流程图算法在计算机上实现快速傅里叶变换必须解决如下问题：

1）、迭代次数r的确定；

由蝶式流程图可见，迭代次数r与N有关。值可由下式确定 

int r = (int)(Math.log10(n)/Math.log10(2)); //求迭代次数r

式中 N 是变换序列的长度。对于前述基数2的蝶式流程图是2的整数次幂。例如，序列长度为8则要三次迭代，序列长度为16时就要4次迭代等等。

2）、对偶节点的计算；

在流程图中把标有xl(k)的点称为节点。其中下标l为列数，也就是第几次迭代，例如，xl(k)则说明它是第一次迭代的结果。k代表流程图中的行数，也就是序列的序号数。其中每一节点的值均是用前一节点对计算得来的。例如，x₁(0)和 x₁(4)均是x(0)和x(4)计算得来的。在蝶式流程图中，把具有相同来源的一对节点叫做对偶节点。

    如: x₁(0)和 x₁(4)就是一对对偶节点，因为它们均来源于x(0)和x(4)。对偶节点的计算也就是求出在每次迭代中对偶节点的间隔或者节距。由流程图可见，第一次迭代的节距为N/2，第二次迭代的节距为N//2²，第三次迭代的节距为N//2³等等。由以上分析可得到如下对偶节点的计算方法。

如果某一节点为x_l(k)，那么，它的对偶节点为

x_l(k+N/2^l)                      (12)

if(k < n/Math.pow(2, l)) {
    x1 = k;
    x2 = x1 + (int)(n/Math.pow(2, l));
} else {
    x2 = k;
    x1 = x2 - (int)(n/Math.pow(2, l));
}

式中 l  是表明第几次迭代的数字，k 是序列的序号数，N  是序列长度。

例：如果序列长度N=8，求x2(1) 的对偶节点。 

x_l(k+N/2l) = x₂((1+8/22) = x2(3)

x₂((1) =x₁((1) + W80x₁(3)

x₂((3) = x₁((1) + W84x₁(3)

3）、加权系数 W_N^P的计算；

 W_N^P的计算主要是确定p值。p值可用下述方法求得。 

  （１）把k值写成r位的二进制数（k是序列的序号数，r是迭代次数）；

  （２）把这个二进制数右移r-l位，并把左边的空位补零（结果仍为r位）；

  （３）把这个右移后的二进制数进行比特倒转；

  （４）把这比特倒转后的二进制数翻成十进制数就得到p值。

例：求x2(2)的加权系数W₈^P。

  （１）因为k=2，所以写成二进制数为010；

  （２）r-l=3-2=1，把010右移一位得到001；

  （３）把001做位序颠倒，即做比特倒转，得到100；

  （４）把100译成十进制数，得到4，所以p=4，x2(2)的加权值为W_N^P。

 结合对偶节点的计算，可以看出W_N^P具有下述规律：如果某一节点上的加权系数为W_N^P，则其对偶节点的加权系数必然是W_N^p+N/2，而且W_N^P = -W_N^p+N/2 ，所以一对对偶节点可用下式计算 

/**
 * 求加权系数
 * 1.将数k写成r位的二进制数;2.将该二进制数向右移r-l位;3.将r位的二进制数比特倒转;4.求出倒置后的二进制数代表的十进制数;
 * @param k 要倒转的十进制数
 * @param l 下标值
 * @param r 二进制的位数
 * @return 加权系数
 */
private int getWeight(int k, int l, int r) {
    int d = r-l;    //位移量
    k = k>>d;
    return reverseRatio(k, r);
}

4）、重新排序问题。

由蝶式流程图可见，如果序列x(n)是按顺序排列的，经过蝶式运算后，其变换序列X(m)是非顺序排列的，即乱序的；反之，如果x(n)是乱序的，那么，就是顺序的。因此，为了便于输出使用，最好加入重新排序程序，以便保证x(n)与它的变换系数 X(m)的对应关系。具体排序方法如下： 

（１）将最后一次迭代结果 xl(k) 中的序号数k写成二进制数;

（２）将r位的二进制数比特倒转： 

（３）求出倒置后的二进制数代表的十进制数，就可以得到与x(k)相对应的X(m)的序号数。 

例如：

   N=8 的最后迭代结果：

     x3(0) →000→倒置→000→十进制（0）

     x3(1) →001→倒置→100→十进制（4）

     x3(2) →010→倒置→010→十进制（2）

     x3(3) →011→倒置→110→十进制（6）

     x3(4) →100→倒置→001→十进制（1）

     x3(5) →101→倒置→101→十进制（5）

     x3(6) →110→倒置→011→十进制（3）

     x3(7) →111→倒置→111→十进制（7）

/**
 * 将数进行二进制倒转， 如0101倒转至1010
 * 1.将数k写成r位的二进制数;2.将r位的二进制数比特倒转;3.求出倒置后的二进制数代表的十进制数;
 * @param k 要倒转的十进制数
 * @param r 二进制的位数
 * @return 倒转后的十进制数
 */
private int reverseRatio(int k, int r) {
    int n = 0;
    StringBuilder sb = new StringBuilder(Integer.toBinaryString(k));
    StringBuilder sb2 = new StringBuilder("");
    if(sb.length()<r) {
        n = r-sb.length();
        for(int i=0; i<n; i++) {
            sb.insert(0, "0");
        }
    }

    for(int i=0; i<sb.length(); i++) {
        sb2.append(sb.charAt(sb.length()-i-1));
    }
    return Integer.parseInt(sb2.toString(), 2);
}

一维傅丽叶变换的源码实现

里面用到了复数的计算，关于复数的计算的源码实现，请参考《 

模拟复数及其运算

》

/**
     * 一维快速傅里叶变换
     * @param values 一维复数集数组
     * @return 傅里叶变换后的数组集
     */
    public Complex[] fft(Complex[] values) {
        int n = values.length;
        int r = (int)(Math.log10(n)/Math.log10(2)); //求迭代次数r
        Complex[][] temp = new Complex[r+1][n]; //计算过程的临时矩阵
        Complex w = new Complex();  //权系数
        temp[0] = values;
        int x1, x2; //一对对偶结点的下标值
        int p, t;   //p表示加权系数Wpn的p值, t是重新排序后对应的序数值
        for(int l=1; l<=r; l++) {
            if(l != r) {
                for(int k=0; k<n; k++) {
                    if(k < n/Math.pow(2, l)) {
                        x1 = k;
                        x2 = x1 + (int)(n/Math.pow(2, l));
                    } else {
                        x2 = k;
                        x1 = x2 - (int)(n/Math.pow(2, l));
                    }
                    p = getWeight(k, l, r);
                    //xi(j) = temp[i-1][x1] + Wpn* temp[i-1][x2];
                    w.setA(Math.cos(-2*Math.PI*p/n));
                    w.setB(Math.sin(-2*Math.PI*p/n));
                    temp[l][k] = Complex.add(temp[l-1][x1] , Complex.multiply(w, temp[l-1][x2]) );

                }
            } else {
                for(int k=0; k<n/2; k++) {
                    x1 = 2*k;
                    x2 = 2*k+1;
                    //System.out.println("x1:" + x1 + "  x2:" + x2);
                    t = reverseRatio(2*k, r);
                    p = t;
                    w.setA(Math.cos(-2*Math.PI*p/n));
                    w.setB(Math.sin(-2*Math.PI*p/n));
                    temp[l][t] = Complex.add(temp[l-1][x1] , Complex.multiply(w, temp[l-1][x2]) );
                    t = reverseRatio(2*k+1, r);
                    p = t;
                    w.setA(Math.cos(-2*Math.PI*p/n));
                    w.setB(Math.sin(-2*Math.PI*p/n));
                    temp[l][t] = Complex.add(temp[l-1][x1] , Complex.multiply(w, temp[l-1][x2]) );
                }
            }
        }
        return temp[r];
    }

二维傅丽叶变换的源码实现

/**
     *  一维快速傅里叶变换
     * @param matrix 二维复数集数组
     * @param w 图像的宽
     * @param h 图像的高
    * @return 傅里叶变换后的数组集
     */
    public Complex[][] fft(Complex matrix[][], int w, int h) {
        double r1 = Math.log10(w)/Math.log10(2.0) - (int)(Math.log10(w)/Math.log10(2.0));
        double r2 = Math.log10(h)/Math.log10(2.0) - (int)(Math.log10(w)/Math.log10(2.0));
        if(r1 != 0.0 || r2 != 0.0) {
            System.err.println("输入的参数w或h不是2的n次幂！");
            return null;
        }
        int r = 0;
        r = (int)(Math.log10(w)/Math.log10(2));
        //进行行傅里叶变换
        for(int i=0; i<h; i++) {
            matrix[i] = fft(matrix[i]);
        }
        //进行列傅里叶变换
        int n = h;
        r = (int)(Math.log10(n)/Math.log10(2)); //求迭代次数r
        Complex tempCom[] = new Complex[h];
        for(int j=0; j<w; j++) {
            for(int i=0; i<h; i++) {
                tempCom[i] = matrix[i][j];
            }
            tempCom = fft(tempCom);
            for(int i=0; i<h; i++) {
                matrix[i][j] = tempCom[i];
            }
        }
        return matrix;
    }