ExampleFFT.cpp

Example for FFT usage

Author

MIP

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/** @example ExampleFFT.cpp
    @relates FFT
    @ingroup g_examples    
    @brief  Example for FFT usage
    @author MIP
*/



#include <MathAlgo/FFT.hh>
#include <Base/Math/Random.hh>
#include <iostream>
#include <iomanip>

using namespace BIAS;
using namespace std;

int main(int argc, char *argv[])
{
  FFT fft;
  Random ran;
  bool debug=false;

  int data_length=100;
  // generate some data by multiplying some sin of random freq.
  int num_sin=5;
  double const_offset=0.0;
  double *sin_freq=new double[num_sin], *phase_offs=new double[num_sin];
  ran.GetUniformDistributed(0.0, 0.5, num_sin, sin_freq);
  ran.GetUniformDistributed(0.0, 2.0*M_PI, num_sin, phase_offs);

  cout << "#using sin waves with freq and phase offset : "<<endl;
  for (int i=0; i<num_sin; i++) 
    cout <<"#freq: "<< sin_freq[i]<<"  \tphase offset: "<<phase_offs[i]<<endl;

  // initialize FFT
  fft.Init(data_length);

  // now generate data
  double *in=fft.GetInput();
  double *p=in;
  for (int l=0; l<data_length; l++, p++){
    *p=const_offset;
    for (int i=0; i<num_sin; i++)
      *p+=sin(2*M_PI*sin_freq[i]*l+phase_offs[i]);
    if (debug)
      cout <<"sample "<< setw(4) << l <<" : "<<setw(15)<<*p<<endl; 
  }

  // now calculate fft
  fft.Compute();
  
  // now get result
  int res_length=(data_length>>1)+1;
  double *phase=new double[res_length], *mag=new double[res_length];
  
  fft.GetPhase(phase);
  fft.GetMagnitude(mag);

  cout << "#result of fft: "<<endl;
  cout <<setw(4)<<"#num"<<setw(15)<<"freq"<<setw(15)<<"mag"<<setw(15)
       <<"phase"<<endl;
  for (int i=0; i<res_length; i++){
    cout <<setw(4)<<i
         <<setw(15)<<(double)i/(double)data_length
         <<setw(15)<<mag[i]
         <<setw(15)<<phase[i]<<"\n";
  }

  if (argc>0)
    cout << "# pipe outputt into file: "<<argv[0]<<" > tmp.dat"<<endl
    << "# and view with gnuplot : echo \"plot \\\"tmp.dat\\\" u 2:3 t \\\"magnitude\\\" w l, \\\"tmp.dat\\\" u 2:4 t \\\"phase\\\" w l; pause 99999;\" |gnuplot"<<endl;

  // or easier
  fft.DumpResult();

  delete[] sin_freq;
  delete[] phase_offs;
  delete[] phase;
  delete[] mag;

  return 0;
}