Test for IDFT.

This commit is contained in:
Themaister 2012-02-28 00:52:49 +01:00
parent c8173aaf62
commit 18aebfedc2

View File

@ -21,6 +21,7 @@
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <string.h>
#include <assert.h>
#include <stdbool.h>
@ -76,7 +77,7 @@ static void interleave(complex double *butterfly_buf, size_t samples)
static complex double gen_phase(double index)
{
return cexp(-M_PI * I * index);
return cexp(M_PI * I * index);
}
static void butterfly(complex double *a, complex double *b, complex double mod)
@ -88,11 +89,11 @@ static void butterfly(complex double *a, complex double *b, complex double mod)
*b = b_;
}
static void butterflies(complex double *butterfly_buf, size_t step_size, size_t samples)
static void butterflies(complex double *butterfly_buf, double phase_dir, size_t step_size, size_t samples)
{
for (unsigned i = 0; i < samples; i += 2 * step_size)
for (unsigned j = i; j < i + step_size; j++)
butterfly(&butterfly_buf[j], &butterfly_buf[j + step_size], gen_phase((double)(j - i) / step_size));
butterfly(&butterfly_buf[j], &butterfly_buf[j + step_size], gen_phase((phase_dir * (j - i)) / step_size));
}
static void calculate_fft(const float *data, complex double *butterfly_buf, size_t samples)
@ -108,15 +109,35 @@ static void calculate_fft(const float *data, complex double *butterfly_buf, size
// Fly, lovely butterflies! :D
for (unsigned step_size = 1; step_size < samples; step_size *= 2)
butterflies(butterfly_buf, step_size, samples);
butterflies(butterfly_buf, -1.0, step_size, samples);
}
// We only have real data.
for (unsigned i = 1; i < samples / 2; i++)
butterfly_buf[i] += conj(butterfly_buf[samples - i]);
static void calculate_fft_adjust(complex double *butterfly_buf, double gain, bool merge_high, size_t samples)
{
if (merge_high)
{
for (unsigned i = 1; i < samples / 2; i++)
butterfly_buf[i] += conj(butterfly_buf[samples - i]);
}
// Normalize amplitudes.
for (unsigned i = 0; i < samples / 2; i++)
butterfly_buf[i] /= (double)samples;
for (unsigned i = 0; i < samples; i++)
butterfly_buf[i] *= gain;
}
static void calculate_ifft(complex double *butterfly_buf, size_t samples, bool normalize)
{
// Enforce POT.
assert((samples & (samples - 1)) == 0);
interleave(butterfly_buf, samples);
// Fly, lovely butterflies! In opposite direction! :D
for (unsigned step_size = 1; step_size < samples; step_size *= 2)
butterflies(butterfly_buf, 1.0, step_size, samples);
if (normalize)
calculate_fft_adjust(butterfly_buf, 1.0 / samples, false, samples);
}
static void test_fft(void)
@ -124,6 +145,7 @@ static void test_fft(void)
fprintf(stderr, "Sanity checking FFT ...\n");
float signal[32];
complex double butterfly_buf[16];
complex double buf_tmp[16];
const float cos_freqs[] = {
1.0, 4.0, 6.0,
@ -143,11 +165,25 @@ static void test_fft(void)
}
calculate_fft(signal, butterfly_buf, 16);
memcpy(buf_tmp, butterfly_buf, sizeof(buf_tmp));
calculate_fft_adjust(buf_tmp, 1.0 / 16, true, 16);
printf("FFT: { ");
for (unsigned i = 0; i < 7; i++)
printf("(%4.2lf, %4.2lf), ", creal(butterfly_buf[i]), cimag(butterfly_buf[i]));
printf("(%4.2lf, %4.2lf) }\n", creal(butterfly_buf[7]), cimag(butterfly_buf[7]));
printf("(%4.2lf, %4.2lf), ", creal(buf_tmp[i]), cimag(buf_tmp[i]));
printf("(%4.2lf, %4.2lf) }\n", creal(buf_tmp[7]), cimag(buf_tmp[7]));
calculate_ifft(butterfly_buf, 16, true);
printf("Original: { ");
for (unsigned i = 0; i < 15; i++)
printf("%5.2f, ", signal[2 * i]);
printf("%5.2f }\n", signal[2 * 15]);
printf("FFT => IFFT: { ");
for (unsigned i = 0; i < 15; i++)
printf("%5.2lf, ", creal(butterfly_buf[i]));
printf("%5.2lf }\n", creal(butterfly_buf[15]));
}
// This doesn't yet take account for slight phase distortions,
@ -158,6 +194,7 @@ static void calculate_snr(struct snr_result *res,
{
samples >>= 1;
calculate_fft(resamp, butterfly_buf, samples);
calculate_fft_adjust(butterfly_buf, 1.0 / samples, true, samples);
double signal = cabs(butterfly_buf[in_rate] * butterfly_buf[in_rate]);
butterfly_buf[in_rate] = 0.0;