/* Gaussian random number sequence generator by Fourier transform (FT) */
/* uniform distribution ---> Guassian distribution                     */
/*           xr[], xi[] ---> xr[], xi[] (inplace processing)           */
/*                                                                     */
/* (c) 2025 cepstrum.co.jp                                             */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>

#define PI     3.14159265359    // pi
#define FFTLEN 1048576          // 1048576=2^20

/*** you can use better FFT routin (FFTW etc) ***/
#include "fft.inc"              // include FFT routine

float xr[FFTLEN];               // FFT work space (real part)
float xi[FFTLEN];               // FFT work space (imaginary part)

#define SEED_VAL  0x5a5aa5a5    // xorshift seed value
#define OUT_FNAME "ft_out.txt"  // output data filne name

//================================================================
// xorshift random number generator (32bit)

uint32_t uint_xorshift32bit(uint32_t seed) {
  static uint32_t y=0x0f0fa5a5;

  if (seed!=0) {
    y=seed;
    return seed;
  }

  y^=(y<<5);
  y^=(y>>13);
  y^=(y<<6);
  return y;
}

//================================================================

int main() {
  FILE     *fp;
  unsigned i;
  float    offset, maxval;

  uint_xorshift32bit(SEED_VAL);      // set xorshift seed

  fft(xr, xi, -1, 1);                // initialize FFT routine

  // generate uniform distribution signal
  for (i=0; i<FFTLEN; i=i+1) xr[i]=uint_xorshift32bit(0);
  for (i=0; i<FFTLEN; i=i+1) xi[i]=uint_xorshift32bit(0);
  // remove DC offset
  offset=0.0;
  for (i=0; i<FFTLEN; i=i+1) offset=offset+xr[i];
  for (i=0; i<FFTLEN; i=i+1) xr[i]=xr[i]-offset/(float)FFTLEN;
  offset=0.0;
  for (i=0; i<FFTLEN; i=i+1) offset=offset+xi[i];
  for (i=0; i<FFTLEN; i=i+1) xi[i]=xi[i]-offset/(float)FFTLEN;
  // normalize [-1.0, 1.0]
  maxval=0.0;
  for (i=0; i<FFTLEN; i=i+1) {
    if (maxval<fabs(xr[i])) maxval=fabs(xr[i]);
  }
  for (i=0; i<FFTLEN; i=i+1) xr[i]=xr[i]/maxval;
  maxval=0.0;
  for (i=0; i<FFTLEN; i=i+1) {
    if (maxval<fabs(xi[i])) maxval=fabs(xi[i]);
  }
  for (i=0; i<FFTLEN; i=i+1) xi[i]=xi[i]/maxval;

  fft(xr, xi, -1, 0);      // FFT

  // output Gaussian random sequence
  fp=fopen(OUT_FNAME, "w");
  for (i=0; i<FFTLEN; i=i+1) fprintf(fp, "%e %e\n", xr[i], xi[i]);
  fclose(fp);

}