#include <math.h>
#include <complex.h>
#include <stdbool.h>
#include <stdint.h>
#include <string.h>
void fft (double complex * out, const double complex * in, int n, bool inverse, int skip) { // use skip=1 for initial calling. internal parameter for recursion.
	if (n == 1) {
		out[0] = in[0];
		return;
	}
	double complex omega = cpow(M_E, -2*M_PI*I/n); // nth root of unity (omega^(n-1)=1)
	if (inverse)
		omega = conj(omega);
	// fprintf(stderr, "omega je %lf+%lfi, n je %d\n", creal(omega), cimag(omega), n);
	fft(out, in, n/2, inverse, skip*2);
	fft(out+n/2, in+skip, n/2, inverse, skip*2);
	for (int i = 0; i < n/2; i++) {
		double complex sod = out[i];
		double complex lih = out[n/2+i];
		out[i] = (sod + cpow(omega, i)*lih)/(skip == 1 && inverse ? n : 1);
		out[n/2+i] = (sod - cpow(omega, i)*lih)/(skip == 1 && inverse ? n : 1);
	}
}
double complex qam (int stopnja /* koren orderja */ , int simbol) {
	double x = simbol%stopnja-(stopnja-1.0)/2.0;
	double y = simbol/stopnja-(stopnja-1.0)/2.0;
	return y*I+x;
}
int qam_symbols (int stopnja) { // how many possible symbols does this qam configuration let you use
	return stopnja*stopnja;
}
int ofdm_columns (int vzorcev, int skip /* 1 za 0 razmika */) {
	return 1+((vzorcev-1)/skip);
}
void ofdm (double complex * out /* space for vzorcev values */, int vzorcev /* must be power of 2 */, int skip, int stopnja, uint64_t simbol) {
	memset(out, 0, sizeof out[0] * vzorcev);
	for (int i = 0; i < ofdm_columns(vzorcev, skip); i++) {
		out[i*skip] = qam(stopnja, simbol % qam_symbols(stopnja));
		simbol /= qam_symbols(stopnja);
	}
} // returns frequency domain, run ifft to get complex time domain, which then has to be moduliran
void moduliraj (double * out, double complex * in, int insize, int faktor) {
	for (int i = 0; i < insize*faktor; i++)
		out[i] = in[i%insize]*sin(M_PI*2*i/(insize*faktor))+in[i%insize]*cos(M_PI*2*i/(insize*faktor));
}
#ifdef MODEMTEST
#include <stdio.h>
#include <error.h>
#include <stdlib.h>
#include <unistd.h>
int main (int argc, char ** argv) {
	if (argc != 5)
		error(1, 0, "%s log_2(vzorcev) ofdm_skip stopnja faktor_modulacije", argv[0]);
	int vzorcev = 1 << atoi(argv[1]);
	int skip = atoi(argv[2]);
	int stopnja = atoi(argv[3]);
	int faktor = atoi(argv[4]);
	int simbolov = pow(qam_symbols(stopnja), ofdm_columns(vzorcev, skip));
	fprintf(stderr, "s temi nastavitvami je stoplcev %d, vsak nosi %d simbolov, skupaj je torej na voljo %d simbolov\n", ofdm_columns(vzorcev, skip), qam_symbols(stopnja), simbolov);
	while (true)
		for (int simbol = 0; simbol < simbolov; simbol++) {
			double complex frequency[vzorcev];
			double complex time[vzorcev];
			double modulirano[vzorcev*faktor];
			ofdm(frequency, vzorcev, skip, stopnja, simbol);
			fft(time, frequency, vzorcev, true, 1);
			moduliraj(modulirano, time, vzorcev, faktor);
			write(STDOUT_FILENO, modulirano, sizeof modulirano);
		}
}
#endif
#ifdef FFTTEST
#include <stdio.h>
int main () {
	printf("fft test.\n");
	double complex sinusoid[128];
	for (int i = 0; i < 128; i++) {
		sinusoid[i] = cpow(M_E, 2*M_PI*I/8*i);
		for (int j = 0; j < 16+creal(sinusoid[i])*16; j++) printf("#");
		printf("\n");
	}
	double complex freq[128];
	printf("fft:\n");
	fft(freq, sinusoid, 128, false, 1);
	for (int i = 0; i < 128; i++) {
		for (int j = 0; j < 16+creal(freq[i]); j++) printf("#");
		printf("\n");
	}
	printf("ifft:\n");
	fft(sinusoid, freq, 128, true, 1);
	for (int i = 0; i < 128; i++) {
		for (int j = 0; j < 16+creal(sinusoid[i])*16; j++) printf("#");
		printf("\n");
	}
	double complex testdata[] = {1, 2, 3, 4};
	double complex output[4];
	fft(output, testdata, 4, false, 1);
	fft(testdata, output, 4, true, 1);
	printf("fftd:\n");
	for (int i = 0; i < 4; i++)
		printf("%lf+%lfi\n", creal(output[i]), cimag(output[i]));
	printf("ifftd:\n");
	for (int i = 0; i < 4; i++)
		printf("%lf+%lfi\n", creal(testdata[i]), cimag(testdata[i]));
}
#endif