实序列快速傅里叶变换(一)

一、功能

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

二、方法简介

实序列\(x(n)\)的离散傅立叶变换为
\[ X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk} \ , \ k=0,1,...,N-1 \]
上式可用复序列FFT算法进行计算。但考虑到\(x(n)\)是实数,为进一步提高计算效率,需要对按时间抽取的基2复序列FFT算法进行一定的修改。

如果序列\(x(n)\)是实数,那么其傅立叶变换\(X(k)\)一般是复数,但其实部是偶对称,虚部是奇对称,即\(X(k)\)具有如下共辄对称性: \(X(0)\)和\(X(N/2)\)都是实数,且有
\[ X(k)=X^{*}(N-k) \ , \ 1 \leqslant k \leqslant \frac{N}{2} - 1 \]
在计算离散傅立叶变换时,利用这种共辄对称性,我们就可以不必计算与存储\(X(k)(N/2 + 1 \leqslant k \leqslant N — 1)\)以及\(X(0)\)和\(X(N/2)\)的虚部,而仅需计算\(X(0)\)到\(X(N/2)\)即可。此处我们选择的是计算\(X(0)\)到\(X(N/4)\)和\(X(N/2)\)到\(X(3N/4)\), 这样做可以恰好利用复序列FFT 算法的前\((N/4)+1\)个复数蝶形。这就是按时间抽取的基2实序列FFT算法,它比复序列FFT算法大约可减少一半的运算量和存储量。

三、使用方法

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

/************************************
    x       ----长度为n。开始时存放要变换的实数据,最后存放变换结果的前n/2+1个值,
                其存储顺序为[Re(0),Re(1),...,Re(n/2),Im(n/2-1),...,Im(1)]。
                其中Re(0)=X(0),Re(n/2)=X(n/2)。根据X(k)的共轭对称性,很容易写
                出后半部分的值。
    n       ----数据长度,必须是2的整数次幂,即n=2^m。
************************************/
#include "math.h"

void rfft(double *x, int n)
{
    int i, j, k, m, il, i2, i3, i4, nl, n2, n4;
    double a, e, cc, ss, xt, tl, t2;

    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) {
            xt = x[j];
            x[j] = x[i];
            x[i] = xt;
        }
        k = n / 2;
        while(k < (j + 1)) {
            j = j - k;
            k = k / 2;
        }
        j = j + k;
    }
    for(i = 0; i < n; i += 2) {
        xt = x[i];
        x[i] = xt + x[i + 1];
        x[i + 1] = xt - x[i + 1];
    }
    n2 = 1;
    for(k = 2; k <= m; k++) {
        n4 = n2;
        n2 = 2 * n4;
        n1 = 2 * n2;
        e = 6.28318530718 / nl;
        for(i = 0; i < n; i += n1) {
            xt = x[i];
            x[i] = xt + x[i + n2];
            x[i + n2] = xt - x[i + n2];
            x[i + n2 + n4] = -x[i + n2 + n4];
            a = e;
            for(j = 1; j <= (n4-1); j++) {
                i1 = i + j;
                i2 = i - j + n2;
                i3 = i + j + n2;
                i4 = i - j + n1;
                cc = cos(a);
                ss = sin(a);
                a = a + e;
                t1 = cc * x[i3] + ss * x[i4];
                t2 = ss * x[i3] - cc * x[i4];
                x[i4] = x[i2] - t2;
                x[i3] = -x[i2] - t2;
                x[i2] = x[i1] - t1;
                x[i1] = x[i1] + t1;
            }
        }
    }
}

声明:该文章系转载,转载该文章的目的在于更广泛的传递信息,并不代表本网站赞同其观点,文章内容仅供参考。

本站是一个个人学习和交流平台,网站上部分文章为网站管理员和网友从相关媒体转载而来,并不用于任何商业目的,内容为作者个人观点, 并不代表本网站赞同其观点和对其真实性负责。

我们已经尽可能的对作者和来源进行了通告,但是可能由于能力有限或疏忽,导致作者和来源有误,亦可能您并不期望您的作品在我们的网站上发布。我们为这些问题向您致歉,如果您在我站上发现此类问题,请及时联系我们,我们将根据您的要求,立即更正或者删除有关内容。本站拥有对此声明的最终解释权。