/* From: "A SURVEY OF COMPUTATIONAL PHYSICS" 
   by RH Landau, MJ Paez, and CC BORDEIANU 
   Copyright Princeton University Press, Princeton, 2008.
   Electronic Materials copyright: R Landau, Oregon State Univ, 2008;
   MJ Paez, Univ Antioquia, 2008; and CC BORDEIANU, Univ Bucharest, 2008.
   Support by National Science Foundation                              
   */
// Noise.java:  Power spectrum of function + noise  
import java.io.*; 

public class Noise  {
  static final int max =1000;                                  
  static double array[] = new double[max], ps[] = new double[max]; 
  static double Autocorr[] = new double[max], twopi = 2*Math.PI; 
  static double step = 4*Math.PI/1000;    
    
  public static void main(String[] argv)  {
    double dftreal[] = new double[max], dftimag[] = new double[max]; 
    discrete(array, max); 
    autocorr(array); 
    fourier(dftreal, dftimag); 
    filter1();       
  }
      
  public static void fourier(double dftreal[], double dftimag[])  {
    double real, imag; 
    int n, k; 
    for ( n = 0;  n < max;  n++ )  {                         // Loop over frequencies
      real = imag = 0. ;                                           // Clear variables
      for ( k = 0;  k < max;  k++ )  {                               // Loop for sums
        real += Autocorr[k]*Math.cos((twopi*k*n)/max); 
        imag += Autocorr[k]*Math.sin((twopi*k*n)/max);  
      }
      dftreal[n] = (1/Math.sqrt(twopi))*real; 
      dftimag[n] = (1/Math.sqrt(twopi))*imag; 
      ps[n] = Math.abs(dftreal[n]*dftreal[n]+dftimag[n]*dftimag[n])/max; 
    }
  }
  
  public static double func(double t)
    { return  1./(1. + 0.9*Math.sin(t)); }
    
 // Setup, plot function + - random noise  
  public static void discrete(double [] array, int max)  {
    double f[] = new double[max + 1], t = 0.; 
    int i;
    for ( i=0;  i < max;  i++ )  {                              // Function + - noise
      f[i] = func(t); 
      array[i] = f[i]+0.5*(2*Math.random()-1); 
      t += step;                                   
    }
  }
                                     // Low-pass windowed-sinc filter, max-100 output
  public static void filter1()  {
    double y[] = new double[max], h[] = new double[100], fc = 0.07;  // Output,filter
    int i, n, m = 100;                              // Set filter length (101 points)
                           // Set cutoff frequency (between 0 & 0.5)
    for ( i = 0;  i < 100; i++ )  {                    // Calc low-pass filter kernel
       // inverse DFT of rectangular pulse, f(t) = sinc = sin x/x
     if ((i-(m/2)) == 0) h[i] = twopi*fc; 
     if ((i-(m/2)) != 0) h[i] = Math.sin(twopi*fc*(i-m/2)) / (i-m/2); 
     h[i] = h[i]*(0.54 - 0.46*Math.cos(twopi*i/m) );                // Hamming window
    }
    double sum = 0. ;                            // Normalize low-pass filter kernel
    for ( i = 0; i < 100; i++ )   sum = sum + h[i];
    for ( i = 0;  i <  100; i++ ) h[i] = h[i] / sum; 
    for ( n = 100 ; n <  max-1; n++ )  {                        // Convolute & filter
      y[n] = 0. ;                   
      for ( i = 0;  i <  100; i++ )  y[n] = y[n] + array[n-i] * h[i]; 
    }
  }
                                                // Calculate autocorrelation function
  public static void autocorr(double[]array)  {
    int i, n;   
    double tau = 0., t, sum; 
    for ( i=0; i < max;  i++ )  {
      sum = 0. ; 
      t = 0. ; 
      for ( n=0; n < max-1; n++ )  {                         // Trap integration rule
        sum  += (func(t)*func(t +tau) + func(t +step)*func(t + step + tau))*step/2.;
        t  += step;         
      }   
      Autocorr[i] = sum; 
      tau += step;         
} } }