/* 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                              
   */
 
// Filter.java Low-pass windowed - sinc filter 
import java.io.*;
import ptolemy.plot.*;

public class Filter  {
  static final int max = 800;                                 
  static double array[] = new double[max];
  static double ps[]=new double[max];
  static double y[]=new double[max];//output
  static double h[]=new double[800];//filter kernel
  public static void main(String[] argv)  {
	  int i, n, k;
    double step = 2*Math.PI/1000;
    double dftreal[] = new double[max];
    double dftimag[] = new double[max];
    function();
    fourier(dftreal,dftimag);
    filter1(h);
  }
    
  public static void filter1(double h[])  {
  //Low-pass windowed - sinc filter 
  //Filters max samples with a 101 point windowed-sinc filter, 
  //resulting in max 100 samples of filtered data.          
    double step = 2*Math.PI/1000;
    int i, j;
    Plot filter = new Plot();
    //double fc = .14;                // Set the cutoff frequency (between 0 and 0.5)
    double fc = .07;  
    int m = 100;                                    // Set filter length (101 points)
    for( i = 0; i<100;i++) {                  // Calculate the low-pass filter kernel
      // inverse DFT of rectangular pulse = sinc function sinx/x
      // h[i]=Math.sin(2.0*Math.PI*fc*i)/i*Math.PI;
      if ((i-(m/2)) == 0)  h[i] = 2*Math.PI*fc;
      if ((i-(m/2))!= 0)   h[i] = Math.sin(2*Math.PI*fc * (i-m/2)) / (i-m/2);
      h[i] = h[i] * (0.54 - 0.46*Math.cos(2*Math.PI*i/m) );// Hamming window
    }
    double sum = 0.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( j = 100 ;j<max-1;j++) {             // Convolve input signal & filter kernel
      y[j] = 0.0;                  
      for (i = 0; i<100;i++)  y[j]+=h[i];
    }
    for (j = 0; j < max-1; j++) filter.addPoint(0, j, h[j], true); 
    filter.setTitle("Windowed sinc-filter kernel");
    filter.setXLabel("Sample Number");
    filter.setYLabel("Amplitude");
    filter.setXRange(0.0, 150.0);
    PlotApplication app4 = new PlotApplication(filter);
    }
    
  public static void fourier(double dftreal[],double dftimag[])  {
    double real, imag;
    int i,j, k, m = 100;                            // Set filter length (101 points)
    Plot powerspec = new Plot();
    double step = 2*Math.PI/1000, fc = .07;  
    for ( i = 0; i<100;i++)  {                // Calculate the low-pass filter kernel
      // inverse DFT of a rectangular pulse corresponds to a sinc function sinx/x
      // h[i]=Math.sin(2.0*Math.PI*fc*i)/i*Math.PI;
      if ((i-(m/2)) == 0)   h[i] = 2*Math.PI*fc;
      if ((i-(m/2)) != 0)   h[i] = Math.sin(2*Math.PI*fc * (i-m/2)) / (i-m/2);
      h[i] = h[i] * (0.54 - 0.46*Math.cos(2*Math.PI*i/m) );         // Hamming window
    }
    double sum = 0.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 (j = 0; j<100; j++)   {                           // loop for frequency index
      real = imag = 0.0;                                           // Clear variables
      for(k = 0; k<100; k++) {                                       // Loop for sums
        real += h[k]*Math.cos((2*Math.PI*k*j)/max);
        imag += h[k]*Math.sin((2*Math.PI*k*j)/max);
      }
		dftreal[j] = (1/Math.sqrt(2*Math.PI))*real;
		dftimag[j] = (1/Math.sqrt(2*Math.PI))*imag;
		ps[j]=(Math.sqrt(dftreal[j]*dftreal[j]+ dftimag[j]*dftimag[j]));
		powerspec.addPoint(0, (0.5*j-1.9)/100, 2.5*ps[j], false);
		}
    powerspec.setTitle("Power Spectrum");
    powerspec.setXLabel("Frequency");
    powerspec.setYLabel("P(w)");
    // powerspec.setXRange(-1.0, 45.0);
    powerspec.setImpulses(true, 0);
    PlotApplication app1 = new PlotApplication(powerspec);
  }
  
  public static void function()  {
    double f[] = new double[max+1];
    int i;
    double step = 2*Math.PI/10000, x = 0.;
    Plot initfunc = new Plot();
    Plot func = new Plot();
    for (i=0; i<max; i++)  {                                      // Initial function
      if ((i>=0)  && (i<500))  f[i]=1.0;
      if ((i>500) && (i<max))  f[i]=0.0;
      if(i==500) f[i]=0.0;
      initfunc.addPoint(0, x, f[i], true);
      x+=step;
    }
    initfunc.setTitle("Ideal frequency response");
    initfunc.setXLabel("Frequency");
    initfunc.setYLabel("Amplitude");
    PlotApplication app3 = new PlotApplication(initfunc);
} }