时间戳

首页 > 解决方案 > 实施IFFT

c - 实施IFFT

问题描述

尝试实施IFFT。如果我将输入设置为此：

{0, 1, 0, 0, 0, 0, 0, 1}

它应该将其转换为具有最低波数的余弦。相反，我得到了这个奇怪的模式：

bin 0: 2.000000     I 0.000000
bin 1: 0.000000     I 1.414214
bin 2: 0.000000     I 0.000000
bin 3: 0.000000     I 1.414214
bin 4: -2.000000    I 0.000000
bin 5: 0.000000     I -1.414214
bin 6: 0.000000     I 0.000000
bin 7: -0.000000    I -1.414214

我只是采用了 FFT 代码和共轭旋转因子。我告诉这应该可以正常工作。（除了比例因子）我不明白出了什么问题。

主程序

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "utils.h"

COMPLEX_FLOAT * coefficients;

void calculate_coefficients(int N)
{
    COMPLEX_FLOAT * co = (COMPLEX_FLOAT *)malloc(N/2 * sizeof(COMPLEX_FLOAT));
    coefficients = co;
    for(int i = 0; i < N/2; i++)
    {
        co[i] = cf_exp(-2.0f*PI*((float)i)/((float)N));
    }
}

COMPLEX_FLOAT* ifftr(COMPLEX_FLOAT* input, int N, int multiplier, int offset)
{
    COMPLEX_FLOAT *output;
    COMPLEX_FLOAT *E, *O;

    output = (COMPLEX_FLOAT *)malloc(N * sizeof(COMPLEX_FLOAT));

    if(N == 1)
    {
        output[0] = input[offset];
    }
    else
    {
        E = fftr(input, N/2, multiplier*2, offset);
        O = fftr(input, N/2, multiplier*2, multiplier + offset);
        for(int i = 0; i < N/2; i++)
        {
            int index = i * multiplier;
            COMPLEX_FLOAT tmp = cmul(conjugate(coefficients[index]),O[i]);
            output[i] = cadd(E[i], tmp);
            output[i + N/2] = csub(E[i], tmp);
        }
        free(E);
        free(O);
    }
    return output;
}

void ifft(COMPLEX_FLOAT* input, COMPLEX_FLOAT* output, int N)
{
    COMPLEX_FLOAT * out;
    calculate_coefficients(N);
    out = ifftr(input,N,1,0);
    for(int i = 0; i < N; i++)
    {
        output[i] = out[i];
    }
    free(out);
    free(coefficients);
    return;
}

int main()
{
    int size = 8;
    COMPLEX_FLOAT dummy[size];
    COMPLEX_FLOAT frq[size];

    for(int i = 0; i < size; i++)
    {
       frq[i].imag = 0;
       if(i == 2)
       {
            frq[i].real = 1;
       }
       else if(size - i == 1)
       {
            frq[i].real = 1;
       }
       else
       {
            frq[i].real = 0;
       }
    }

    ifft(frq, dummy, size);

    for(int i = 0; i < size; i++)
    {
        printf("bin %d: %f\t I %f\n", i, dummy[i].real, dummy[i].imag);
    }

    return 0;
}

实用程序.c


#include "utils.h"
#include "math.h"

COMPLEX_FLOAT cadd(COMPLEX_FLOAT a, COMPLEX_FLOAT b)
{
    COMPLEX_FLOAT out;
    out.real = a.real + b.real;
    out.imag = a.imag + b.imag;
    return out;
}

COMPLEX_FLOAT csub(COMPLEX_FLOAT a, COMPLEX_FLOAT b)
{
    COMPLEX_FLOAT out;
    out.real = a.real - b.real;
    out.imag = a.imag - b.imag;
    return out;
}

COMPLEX_FLOAT cmul(COMPLEX_FLOAT a, COMPLEX_FLOAT b)
{
    COMPLEX_FLOAT out;
    out.real = a.real * b.real - a.imag * b.imag;
    out.imag = a.real * b.imag + b.real * a.imag;
    return out;
}

COMPLEX_FLOAT cf_exp(float a)
{
    COMPLEX_FLOAT out;
    out.real = (float)cos(a);
    out.imag = (float)sin(a);
    return out;
}

COMPLEX_FLOAT conjugate(COMPLEX_FLOAT a)
{
    COMPLEX_FLOAT out;
    out.real = a.real;
    out.imag = -a.imag;
    return out;
}

COMPLEX_FLOAT real_num(float n)
{
    COMPLEX_FLOAT out;
    out.real = n;
    out.imag = 0;
    return out;
}

标签： cfft

解决方案

我调用了错误的递归函数fftr而不是ifftr

COMPLEX_FLOAT* ifftr(COMPLEX_FLOAT* input, int N, int multiplier, int offset)
{
    COMPLEX_FLOAT *output;
    COMPLEX_FLOAT *E, *O;

    output = (COMPLEX_FLOAT *)malloc(N * sizeof(COMPLEX_FLOAT));

    if(N == 1)
    {
        output[0] = input[offset];
    }
    else
    {
        E = ifftr(input, N/2, multiplier*2, offset);
        O = ifftr(input, N/2, multiplier*2, multiplier + offset);
        for(int i = 0; i < N/2; i++)
        {
            int index = i * multiplier;
            COMPLEX_FLOAT tmp = cmul(conjugate(coefficients[index]),O[i]);
            output[i] = cadd(E[i], tmp);
            output[i + N/2] = csub(E[i], tmp);
        }
        free(E);
        free(O);
    }
    return output;
}

推荐阅读