傅丽叶变换（二）

傅丽叶变换（二）

——(java)算法实现

离散傅里叶变换

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

离散傅里叶变换的定义

如果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;
	}