Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

IFFT error #351

Open
Apriqi opened this issue Apr 19, 2024 · 3 comments
Open

IFFT error #351

Apriqi opened this issue Apr 19, 2024 · 3 comments

Comments

@Apriqi
Copy link

Apriqi commented Apr 19, 2024

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

@Apriqi
Copy link
Author

Apriqi commented Apr 19, 2024

#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
Copy link
Contributor

rdolbeau commented Jul 3, 2024

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

@stevengj
Copy link
Contributor

stevengj commented Jul 3, 2024

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants