Chris@32: /* 
Chris@32:  * FFT and convolution test (C)
Chris@32:  * 
Chris@32:  * Copyright (c) 2014 Project Nayuki
Chris@32:  * http://www.nayuki.io/page/free-small-fft-in-multiple-languages
Chris@32:  * 
Chris@32:  * (MIT License)
Chris@32:  * Permission is hereby granted, free of charge, to any person obtaining a copy of
Chris@32:  * this software and associated documentation files (the "Software"), to deal in
Chris@32:  * the Software without restriction, including without limitation the rights to
Chris@32:  * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
Chris@32:  * the Software, and to permit persons to whom the Software is furnished to do so,
Chris@32:  * subject to the following conditions:
Chris@32:  * - The above copyright notice and this permission notice shall be included in
Chris@32:  *   all copies or substantial portions of the Software.
Chris@32:  * - The Software is provided "as is", without warranty of any kind, express or
Chris@32:  *   implied, including but not limited to the warranties of merchantability,
Chris@32:  *   fitness for a particular purpose and noninfringement. In no event shall the
Chris@32:  *   authors or copyright holders be liable for any claim, damages or other
Chris@32:  *   liability, whether in an action of contract, tort or otherwise, arising from,
Chris@32:  *   out of or in connection with the Software or the use or other dealings in the
Chris@32:  *   Software.
Chris@32:  */
Chris@32: 
Chris@32: #include <math.h>
Chris@32: #include <stdio.h>
Chris@32: #include <stdlib.h>
Chris@32: #include <string.h>
Chris@32: #include <time.h>
Chris@32: #include "fft.h"
Chris@32: 
Chris@32: 
Chris@32: // Private function prototypes
Chris@32: static void test_fft(int n);
Chris@32: static void test_convolution(int n);
Chris@32: static void naive_dft(const double *inreal, const double *inimag, double *outreal, double *outimag, int inverse, int n);
Chris@32: static void naive_convolve(const double *xreal, const double *ximag, const double *yreal, const double *yimag, double *outreal, double *outimag, int n);
Chris@32: static double log10_rms_err(const double *xreal, const double *ximag, const double *yreal, const double *yimag, int n);
Chris@32: static double *random_reals(int n);
Chris@32: static void *memdup(const void *src, size_t n);
Chris@32: 
Chris@32: static double max_log_error = -INFINITY;
Chris@32: 
Chris@32: 
Chris@32: /* Main and test functions */
Chris@32: 
Chris@32: int main(int argc, char **argv) {
Chris@32: 	int i;
Chris@32: 	int prev;
Chris@32: 	srand(time(NULL));
Chris@32: 	
Chris@32: 	// Test power-of-2 size FFTs
Chris@32: 	for (i = 0; i <= 12; i++)
Chris@32: 		test_fft(1 << i);
Chris@32: 	
Chris@32: 	// Test small size FFTs
Chris@32: 	for (i = 0; i < 30; i++)
Chris@32: 		test_fft(i);
Chris@32: 	
Chris@32: 	// Test diverse size FFTs
Chris@32: 	prev = 0;
Chris@32: 	for (i = 0; i <= 100; i++) {
Chris@32: 		int n = (int)lround(pow(1500, i / 100.0));
Chris@32: 		if (n > prev) {
Chris@32: 			test_fft(n);
Chris@32: 			prev = n;
Chris@32: 		}
Chris@32: 	}
Chris@32: 	
Chris@32: 	// Test power-of-2 size convolutions
Chris@32: 	for (i = 0; i <= 12; i++)
Chris@32: 		test_convolution(1 << i);
Chris@32: 	
Chris@32: 	// Test diverse size convolutions
Chris@32: 	prev = 0;
Chris@32: 	for (i = 0; i <= 100; i++) {
Chris@32: 		int n = (int)lround(pow(1500, i / 100.0));
Chris@32: 		if (n > prev) {
Chris@32: 			test_convolution(n);
Chris@32: 			prev = n;
Chris@32: 		}
Chris@32: 	}
Chris@32: 	
Chris@32: 	printf("\n");
Chris@32: 	printf("Max log err = %.1f\n", max_log_error);
Chris@32: 	printf("Test %s\n", max_log_error < -10 ? "passed" : "failed");
Chris@32: 	return 0;
Chris@32: }
Chris@32: 
Chris@32: 
Chris@32: static void test_fft(int n) {
Chris@32: 	double *inputreal, *inputimag;
Chris@32: 	double *refoutreal, *refoutimag;
Chris@32: 	double *actualoutreal, *actualoutimag;
Chris@32: 	
Chris@32: 	inputreal = random_reals(n);
Chris@32: 	inputimag = random_reals(n);
Chris@32: 	
Chris@32: 	refoutreal = malloc(n * sizeof(double));
Chris@32: 	refoutimag = malloc(n * sizeof(double));
Chris@32: 	naive_dft(inputreal, inputimag, refoutreal, refoutimag, 0, n);
Chris@32: 	
Chris@32: 	actualoutreal = memdup(inputreal, n * sizeof(double));
Chris@32: 	actualoutimag = memdup(inputimag, n * sizeof(double));
Chris@32: 	transform(actualoutreal, actualoutimag, n);
Chris@32: 	
Chris@32: 	printf("fftsize=%4d  logerr=%5.1f\n", n, log10_rms_err(refoutreal, refoutimag, actualoutreal, actualoutimag, n));
Chris@32: 	
Chris@32: 	free(inputreal);
Chris@32: 	free(inputimag);
Chris@32: 	free(refoutreal);
Chris@32: 	free(refoutimag);
Chris@32: 	free(actualoutreal);
Chris@32: 	free(actualoutimag);
Chris@32: }
Chris@32: 
Chris@32: 
Chris@32: static void test_convolution(int n) {
Chris@32: 	double *input0real, *input0imag;
Chris@32: 	double *input1real, *input1imag;
Chris@32: 	double *refoutreal, *refoutimag;
Chris@32: 	double *actualoutreal, *actualoutimag;
Chris@32: 	
Chris@32: 	input0real = random_reals(n);
Chris@32: 	input0imag = random_reals(n);
Chris@32: 	input1real = random_reals(n);
Chris@32: 	input1imag = random_reals(n);
Chris@32: 	
Chris@32: 	refoutreal = malloc(n * sizeof(double));
Chris@32: 	refoutimag = malloc(n * sizeof(double));
Chris@32: 	naive_convolve(input0real, input0imag, input1real, input1imag, refoutreal, refoutimag, n);
Chris@32: 	
Chris@32: 	actualoutreal = malloc(n * sizeof(double));
Chris@32: 	actualoutimag = malloc(n * sizeof(double));
Chris@32: 	convolve_complex(input0real, input0imag, input1real, input1imag, actualoutreal, actualoutimag, n);
Chris@32: 	
Chris@32: 	printf("convsize=%4d  logerr=%5.1f\n", n, log10_rms_err(refoutreal, refoutimag, actualoutreal, actualoutimag, n));
Chris@32: 	
Chris@32: 	free(input0real);
Chris@32: 	free(input0imag);
Chris@32: 	free(input1real);
Chris@32: 	free(input1imag);
Chris@32: 	free(refoutreal);
Chris@32: 	free(refoutimag);
Chris@32: 	free(actualoutreal);
Chris@32: 	free(actualoutimag);
Chris@32: }
Chris@32: 
Chris@32: 
Chris@32: /* Naive reference computation functions */
Chris@32: 
Chris@32: static void naive_dft(const double *inreal, const double *inimag, double *outreal, double *outimag, int inverse, int n) {
Chris@32: 	double coef = (inverse ? 2 : -2) * M_PI;
Chris@32: 	int k;
Chris@32: 	for (k = 0; k < n; k++) {  // For each output element
Chris@32: 		double sumreal = 0;
Chris@32: 		double sumimag = 0;
Chris@32: 		int t;
Chris@32: 		for (t = 0; t < n; t++) {  // For each input element
Chris@32: 			double angle = coef * ((long long)t * k % n) / n;
Chris@32: 			sumreal += inreal[t]*cos(angle) - inimag[t]*sin(angle);
Chris@32: 			sumimag += inreal[t]*sin(angle) + inimag[t]*cos(angle);
Chris@32: 		}
Chris@32: 		outreal[k] = sumreal;
Chris@32: 		outimag[k] = sumimag;
Chris@32: 	}
Chris@32: }
Chris@32: 
Chris@32: 
Chris@32: static void naive_convolve(const double *xreal, const double *ximag, const double *yreal, const double *yimag, double *outreal, double *outimag, int n) {
Chris@32: 	int i;
Chris@32: 	for (i = 0; i < n; i++) {
Chris@32: 		double sumreal = 0;
Chris@32: 		double sumimag = 0;
Chris@32: 		int j;
Chris@32: 		for (j = 0; j < n; j++) {
Chris@32: 			int k = (i - j + n) % n;
Chris@32: 			sumreal += xreal[k] * yreal[j] - ximag[k] * yimag[j];
Chris@32: 			sumimag += xreal[k] * yimag[j] + ximag[k] * yreal[j];
Chris@32: 		}
Chris@32: 		outreal[i] = sumreal;
Chris@32: 		outimag[i] = sumimag;
Chris@32: 	}
Chris@32: }
Chris@32: 
Chris@32: 
Chris@32: /* Utility functions */
Chris@32: 
Chris@32: static double log10_rms_err(const double *xreal, const double *ximag, const double *yreal, const double *yimag, int n) {
Chris@32: 	double err = 0;
Chris@32: 	int i;
Chris@32: 	for (i = 0; i < n; i++)
Chris@32: 		err += (xreal[i] - yreal[i]) * (xreal[i] - yreal[i]) + (ximag[i] - yimag[i]) * (ximag[i] - yimag[i]);
Chris@32: 	
Chris@32: 	err /= n > 0 ? n : 1;
Chris@32: 	err = sqrt(err);  // Now this is a root mean square (RMS) error
Chris@32: 	err = err > 0 ? log10(err) : -99.0;
Chris@32: 	if (err > max_log_error)
Chris@32: 		max_log_error = err;
Chris@32: 	return err;
Chris@32: }
Chris@32: 
Chris@32: 
Chris@32: static double *random_reals(int n) {
Chris@32: 	double *result = malloc(n * sizeof(double));
Chris@32: 	int i;
Chris@32: 	for (i = 0; i < n; i++)
Chris@32: 		result[i] = (rand() / (RAND_MAX + 1.0)) * 2 - 1;
Chris@32: 	return result;
Chris@32: }
Chris@32: 
Chris@32: 
Chris@32: static void *memdup(const void *src, size_t n) {
Chris@32: 	void *dest = malloc(n);
Chris@32: 	if (dest != NULL)
Chris@32: 		memcpy(dest, src, n);
Chris@32: 	return dest;
Chris@32: }