diff --git a/audio/test/snr.c b/audio/test/snr.c index 4a9567b841..ab8b20c7dc 100644 --- a/audio/test/snr.c +++ b/audio/test/snr.c @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -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;