radix-4 FFT 原理和C语言代码实现

目录

参考书籍：

1、按照频率抽取

2、按照时间抽取：

3、C代码实现

---

参考书籍：

《数字信号处理》 邓小玲 徐梅宣等主编

1、按照频率抽取

 

2、按照时间抽取：

 

3、C代码实现

//radix-4 FFT 按照频率抽取
/***************************************************************
  *  @brief     按照频率抽取的radix-4 FFT
  *  @param x:  长度n的double型一维数组，开始存储要变换数据的实部，最后存变换结果的实部 
  *  @param y:  长度n的double型一维数组，开始存储要变换数据的虚部，最后存变换结果的虚部 
  *  @param n:  n=4^m, m是正整数
  *  @note     注意基四的FFT要求n是4的整数次幂。radix-4比radix-2乘法量减少25%，加法量略减少；
  *  @Sample usage: 直接调用
 **************************************************************/
void radix_4_FFT(x, y, n)
{
    int n;
    double x[], y[];
    int i, j, k, m, i1, i2, i3, n1, n2;
    double a, b, c, e, r1, r2, r3, r4, s1, s2, s3, s4;
    double co1, co2, co3, si1, si2, si3;
    for(j = 1, i=1; i < 10; ++i)
    {
        m = i;
        j = 4*j;
        if(j == n) break;
	}
    n2 = n;
    for(k = 1; k <= m; ++k)
    {
        n1 = n2;
        n2 = n2 / 4;
        e = 6.28318530718 / n1;
        a = 0;
        for(j = 0; j < n2; ++j)
        {
            b = a + a;
            c = a + b;
            co1 = cos(a);
            co2 = cos(b);
            co3 = cos(c);
            si1 = sin(a);
            si2 = sin(b);
            si3 = sin(c);
            a = (j + 1) * e;
            for(i = j; i < n; i = i + n1)
            {
                i1 = i + n2;
                i2 = i1 + n2;
                i3 = i2 + n2;
                r1 = x[i] + x[i2];
                r3 = x[i] - x[i2];
                s1 = y[i] + y[i2];
                s3 = y[i] - y[i2];
                r2 = x[i1] + x[i3];
                r4 = x[i1] - x[i3];
                s2 = y[i1] + y[i3];
                s4 = y[i1] - y[i3];
                x[i] = r1 - r2;
                r2 = r1 - r2;
                r1 = r3 - s4;
                r3  = r3 + s4;
                y[i] = s1 + s2;
                s2 = s2 - s2;
                s1 = s3 + r4;
                s3 = s3 - r4;
                x[i1] = co1 * r3 + si1 * s3;
                y[i1] = col * s3 - si1 * r3;
                x[i2] = co2 * r3 + si2 * s3;
                y[i2] = co2 * s3 - si2 * r3;
                x[i3] = co3 * r3 + si3 * s3;
                y[i3] = co3 * s3 - si3 * r3;
            }
         }
    }
    n1 = n - 1;
    for(j = 0; i = 0; i < n1; ++i)
    {
        if(i < j)
        {
            r1 = x[j];
            s1 = y[j];
            x[j] = x[i];
            y[j] = y[i];
            x[i] = r1;
            y[i] = s1;
        }
        k = n / 4;
        while(3 * k < (j + 1))
        {
            j = j - 3*k;
            k = k / 4;
        }
        j = j + k;
	}
}