IFFT error
Apriqi opened this issue · 3 comments
Apriqi commented
编译的.so .a文件,进行测试,正变换结果正确,逆变换结果错误。
N_fft=8192
16个正确值+16个错误值+16个正确值+16个错误值+16个正确值+16个错误值+...........
Apriqi commented
#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 commented
Could it just be the scaling ? http://fftw.org/faq/section3.html#whyscaled
stevengj commented
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.