IFFT error

[Apriqi]

编译的.so .a文件，进行测试，正变换结果正确，逆变换结果错误。 N_fft=8192 16个正确值+16个错误值+16个正确值+16个错误值+16个正确值+16个错误值+...........

[Apriqi]

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include "fftw3.h"
#include <iostream>

// 定义 WAV 文件头结构
typedef struct
{
    char chunkId[4];
    uint32_t chunkSize;
    char format[4];
    char subchunk1Id[4];
    uint32_t subchunk1Size;
    uint16_t audioFormat;
    uint16_t numChannels;
    uint32_t sampleRate;
    uint32_t byteRate;
    uint16_t blockAlign;
    uint16_t bitsPerSample;
    char subchunk2Id[4];
    uint32_t subchunk2Size;
} WavHeader;

// 从 WAV 文件中读取数据
int32_t **readWavData(const char *filename, WavHeader *header)
{
    FILE *file = fopen(filename, "rb");
    if (!file)
    {
        printf("Failed to open file.\n");
        return NULL;
    }

    // 读取 WAV 文件头
    fread(header, sizeof(WavHeader), 1, file);

    // 计算每个通道的样本数
    int samplesPerChannel = header->subchunk2Size / (header->numChannels * sizeof(int32_t));

    // 分配存储数据的内存
    int32_t **data = (int32_t **)malloc(header->numChannels * sizeof(int32_t *));
    for (int i = 0; i < header->numChannels; ++i)
    {
        data[i] = (int32_t *)malloc(samplesPerChannel * sizeof(int32_t));
    }

    // 读取音频数据
    for (size_t i = 0; i < samplesPerChannel; ++i)
    {
        for (int j = 0; j < header->numChannels; ++j)
        {
            fread(&data[j][i], sizeof(int32_t), 1, file);
        }
    }

    fclose(file);
    return data;
}

int main()
{
    const char *filename = "AKM_processed.wav";
    WavHeader header;
    int32_t **data = readWavData(filename, &header);
    if (!data)
    {
        std::cout << "Failed to open wave file." << std::endl;
        return 1;
    }

    // 打开输出文件
    FILE *inputFile = fopen("fft_input.txt", "w");
    FILE *outputFile = fopen("fft_output.txt", "w");
    FILE *out_inputFile = fopen("ifft_output.txt", "w");

    // 创建 FFTW 输入数组
    int N = header.subchunk2Size / (header.numChannels * sizeof(int32_t)); // 数据长度

    double *in = (double *)fftw_malloc(sizeof(double) * N);
    fftw_complex *out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (N / 2 + 1));

    fftw_plan dft = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
    fftw_plan idft = fftw_plan_dft_c2r_1d(N, out, in, FFTW_ESTIMATE); // 逆变换结果不对,16个对，16个不对，交叉出现对错

    for (int i = 0; i < header.numChannels; i++)
    {
        for (int j = 0; j < N; j++)
        {
            in[j] = (double)data[i][j];
        }

        for (int j = 0; j < N; j++)
        {
            fprintf(inputFile, "%f", in[j]);
            fprintf(inputFile, "\n");
        }

        fftw_execute(dft);
        fftw_execute(idft);

        for (int j = 0; j < N / 2 + 1; j++)
        {
            fprintf(outputFile, "%.6f + %.6fi ", out[j][0], out[j][1]);
            fprintf(outputFile, "\n");
        }

        for (int j = 0; j < N; j++)
        {
            fprintf(out_inputFile, "%f", in[j] / N);
            fprintf(out_inputFile, "\n");
        }
    }

    // 关闭输出文件
    fclose(inputFile);
    fclose(outputFile);
    fclose(out_inputFile);

    // 释放内存和销毁计划
    fftw_destroy_plan(dft);
    fftw_destroy_plan(idft);
    fftw_free(in);
    fftw_free(out);

    // 释放 WAV 数据的内存
    for (int i = 0; i < header.numChannels; ++i)
    {
        free(data[i]);
    }
    free(data);

    return 0;
}

[rdolbeau]

Could it just be the scaling ? http://fftw.org/faq/section3.html#whyscaled

[stevengj]

The code looks like it has the correct 1/N scaling.

According to Google translate, @Apriqi is seeing:

the inverse transformation result is wrong. N_fft=8192 16 correct values + 16 false values + 16 correct values + 16 false values + 16 correct values + 16 false values +...........

Are you expecting the values to match exactly? They will differ slightly from the original inputs due to floating-point roundoff errors.