快速傅里叶变换

2019/10/22 21:13
阅读数 25

一、功能

计算复序列的快速傅里叶变换。

二、方法简介

序列$x(n)(n=0,1,...,N-1)$的离散傅里叶变换定义为 $$ X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk}, \qquad k=0,1,...,N-1 $$ 其中$W_{N}^{nk}=e^{-j\frac{2\pi nk}{N}}$,将序列$x(n)$按序号$n$的奇偶分成两组,即 $$ \left.\begin{matrix}\begin{align*}x_{1}(n)=&x(2n)\ x_{2}(n)=&x(2n+1)\end{align*}\end{matrix}\right} \qquad n=0,1,...,\frac{N}{2}-1 $$ 因此,$x(n)$的傅里叶变换可写成 $$ \begin{align*}X(k) &= \sum_{n=0}^{N/2-1}x(2n)W^{2nk}{N} + \sum{n=0}^{N/2-1}x(2n+1)W^{(2n+1)k}{N}\&= \sum{n=0}^{N/2-1}x_{1}(n)W^{nk}{N/2} + W{N}^{k}\sum_{n=0}^{N/2-1}x_{2}(n)W^{nk}{N/2}\end{align*} $$ 由此可得$X(k)=X{1}(k)+W_{N}^{k}X_{2}(k), \qquad k = 0,1,...,\frac{N}{2}$,式中 $$ \begin{align*}X_{1}(k)&=\sum_{n=0}^{N/2-1}x(2n)W^{2nk}{N}\X{2}(k)&=\sum_{n=0}^{N/2-1}x(2n+1)W^{(2n+1)k}{N}\end{align*} $$ 他们分别是$x_1(n)$和$x_2(n)$的$N/2$点DFT。上面的推导表明:一个$N$点DFT被分解为两个$N/2$点DFT,这两个$N/2$点DFT又可合成一个$N$点DFT。但上面给出的公式仅能得到$X(k)$的前$N/2$点的值,要用$X{1}(k)$和$X_{2}(k)$来表达$X(k)$的后半部分的值,还必须运用权系数$W_N$的周期性与对称性,即 $$ W_{N/2}^{n(k+N/2)}=W_{N/2}^{nk}, \quad W_{N}^{(k+N/2)}=-W_{N}^{k} $$ 因此,$X(k)$的后$N/2$点的值可表示为 $$ \begin{align*}X(k+\frac{N}{2})&=X_{1}(k+\frac{N}{2})+W_{N}^{k+N/2}X_{2}(k+\frac{N}{2})\&=X_{1}(k)-W_{N}^{k}X_{2}(k), \ k=0,1,...,\frac{N}{2}-1\end{align*} $$ 通过上面的推导可以看出,$N$点的DFT可以分解为两个$N/2$点DFT,每个$N/2$点DFT又可以分解为两个$N/4$点DFT。依此类推,当$N$为2的整数次幂时($N=2^M$),由于每分解一次降低一阶幂次,所以通过$M$次分解,最后全部成为一系列2点DFT运算。以上就是按时间抽取的快速傅里叶变换(FFT)算法。

序列$X(k)$的离散傅里叶反变换定义为 $$ x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X(k)W_{N}^{-nk}, \qquad n=0,1,...,N-1 $$ 它与离散傅里叶正变换的区别在于将$W_N$改变为$W_N^{-1}$,并多了一个除以$N$的运算。因为$W_N$和$W_N^{-1}$对于推导按时间抽取的快速傅里叶变换算法并无实质性区别,因此可将FFT和快速傅里叶反变换(IFFT)算法合并在同一程序中。

三、使用说明

是用C语言实现快速傅里叶变换(FFT)的方法如下:

/************************************
	x       ---一维数组,长度为n,开始时存放要变换数据的实部,最后存放变换结果的实部。
	y       ---一维数组,长度为n,开始时存放要变换数据的虚部,最后存放变换结果的虚部。
	n 		---数据长度,必须是2的整数次幂。
	sign 	---当sign=1时,子函数计算离散傅里叶正变换;当sign=-1时,子函数计算离散傅里叶反变换
************************************/
#include "math.h"

void fft(double *x, double *y, int n, int sign)
{
	int i, j, k, l, m, n1, n2;
	double c, c1, e, s, s1, t, tr, ti;
	for(j = 1, i=1; i < 16; i++) {
		m = i;
		j = 2 * j;
		if(j == n)
			break;
	}
	n1 = n - 1;
	for(j = 0, i = 0; i < n1; i++) {
		if(i < j) {
			tr = x[j];
			ti = j[j];
			x[j] = x[i];
			y[j] = j[i];
			x[i] = tr;
			y[i] = ti;
		}
		k = n / 2;
		while(k < (j + 1)) {
			j = j - k;
			k = k / 2;
		}
		j = j + k;
	}
	n1 = 1;
	for(l = 1; l <= m; l++) {
		n1 = 2 * n1;
		n2 = n1 / 2;
		e = 3.14159265359 / n2;
		c = 1.0;
		s = 0.0;
		c1 = cos(e);
		s1 = -sign * sin(e);
		for(j = 0; j < n2; j++) {
			for(i = j; i < n; i += n1) {
				k = i + n2;
				tr = c * x[k] - s * y(k);
				ti = c * y[k] + s * x[k];
				x[k] = x[i] - tr;
				y[k] = y[i] - ti;
				x[i] = x[i] + tr;
				y[i] = y[i] + ti;
			}
			t = c;
			c = c * c1 - s * s1;
			s = t * s1 + s * c1;
		}
	}
	if(sign == -1) {
		for(i = 0; i < n; i++) {
			x[i] /= n;
			y[i] /= n;
		}
	}

}
展开阅读全文
打赏
0
0 收藏
分享
加载中
更多评论
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部