orsa_fft.cc File Reference

#include <iostream>
#include <complex>
#include <fftw.h>
#include <gsl/gsl_heapsort.h>
#include "orsa_fft.h"
#include "orsa_common.h"

Include dependency graph for orsa_fft.cc:

Go to the source code of this file.

Namespaces

namespace  orsa

Functions

double norm (const fftw_complex z)
double norm_sq (const fftw_complex z)
fftw_complex phi (double omega, fftw_complex in[], const int size)
 Discrete Fourier Transform.
fftw_complex phi_Hanning (double omega, fftw_complex in[], const int size)
 Discrete Fourier Transform with Hanning windowing.
double phi_amp (double omega, fftw_complex in[], const int size)
 amplitude for spectrum, without windowing
double phi_Hanning_amp (double omega, fftw_complex in[], const int size)
 amplitude for spectrum, with Hanning windowing
double phi_gsl (double x, void *params)
double phi_gsl_two (double x, void *params)
double phi_Hanning_gsl (double x, void *params)
int compare_binamp (const binamp *a, const binamp *b)
 sort binamp struct from the bigger to the smaller...
double psd_max_again (const fftw_complex *transformed_signal, const int size)
void psd_max_again_many (const fftw_complex *transformed_signal, const int size, vector< double > &candidate, const unsigned int nfreq)
double psd_max (const fftw_complex *transformed_signal, const int size)
void apply_window (fftw_complex *signal_win, fftw_complex *signal, int size)
void amph (double *amp, double *phase, fftw_complex *phiR, fftw_complex *phiL, double freq, fftw_complex *in, int size)
double accurate_peak (double left, double center, double right, fftw_complex *in, int size)
double dQ (double y)


Function Documentation

double accurate_peak double  left,
double  center,
double  right,
fftw_complex *  in,
int  size
 

Definition at line 479 of file orsa_fft.cc.

References orsa::norm(), orsa::phi(), orsa::phi_gsl_two(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.

00479                                                                                              {
00480     
00481     // check values 
00482     // if ( (center<0) || (center>0.5) ) return center; // ERROR...
00483     if ( (center<(-0.5)) || (center>0.5) ) { 
00484       cerr << "warning!! Peak out of range!" << endl;
00485       return center; // ERROR...
00486     }
00487     
00488     // int    c_iter=0, l_iter=0, r_iter=0, max_iter = 100;
00489     int    max_iter = 100;
00490     double center_amp, left_amp, right_amp;
00491     
00492     left_amp   = norm(phi(left,  in,size));
00493     right_amp  = norm(phi(right, in,size));
00494     center_amp = norm(phi(center,in,size));
00495     
00496     // TODO: check for gsl peak search conditions!!!
00497     
00498     // gsl stuff
00499     gsl_d1_minimization_parameters par;
00500     gsl_function F;
00501     gsl_min_fminimizer *s;
00502     const gsl_min_fminimizer_type *T;
00503     //
00504     T = gsl_min_fminimizer_goldensection;  
00505     s = gsl_min_fminimizer_alloc(T);
00506     //
00507     par.size = size;
00508     par.pointer_points_sequence = &in[0];
00509     //
00510     F.function = &phi_gsl_two;
00511     F.params   = &par;
00512     //
00513     gsl_min_fminimizer_set(s,&F,center,left,right);
00514     
00515     double epsrel = 1.0e-10; // for the d=1 minimization
00516     double epsabs = 0.0;     // for the d=1 minimization
00517     
00518     int iter = 0;
00519     int status;
00520     do { 
00521       iter++;
00522       
00523       status = gsl_min_fminimizer_iterate (s);
00524       //
00525       center = gsl_min_fminimizer_minimum (s);
00526       left   = gsl_min_fminimizer_x_lower (s);
00527       right  = gsl_min_fminimizer_x_upper (s);
00528       // delta     = range_stop - range_start;
00529       // minimize!
00530       status = gsl_min_test_interval(left,right,epsrel,epsabs);
00531       
00532       // post-minimization...
00533       /* printf("%5d [%f, %f] %.7f %.7e  bin:%i\n", 
00534          iter,range_start,range_stop,test_peak.frequency,range_stop-range_start,bin);
00535       */
00536     } while ((status == GSL_CONTINUE) && (iter < max_iter));
00537     
00538     return (center);
00539   }

Here is the call graph for this function:

void amph double *  amp,
double *  phase,
fftw_complex *  phiR,
fftw_complex *  phiL,
double  freq,
fftw_complex *  in,
int  size
 

Definition at line 467 of file orsa_fft.cc.

References orsa::norm(), orsa::phi(), and orsa::secure_atan2().

00467                                                                                                                          {
00468     
00469     *phiR = phi(freq, in,size);
00470     *phiL = phi(-freq,in,size);
00471     
00472     *amp   = norm(*phiR);
00473     *phase = secure_atan2(phiR->im,phiR->re);
00474     
00475     // cerr << "amph() DUMP:   amp = " << *amp << "  phase = " << *phase << endl;
00476     
00477   }

Here is the call graph for this function:

void apply_window fftw_complex *  signal_win,
fftw_complex *  signal,
int  size
 

Definition at line 449 of file orsa_fft.cc.

References orsa::cos(), and orsa::twopi.

00449                                                                               {
00450     
00451     int j;
00452     double win_coeff,arg;
00453     
00454     for(j=0;j<size;j++) {
00455       
00456       arg = twopi*j/size;
00457       win_coeff = (1. - cos(arg));
00458       
00459       signal_win[j].re = signal[j].re * win_coeff;
00460       signal_win[j].im = signal[j].im * win_coeff;
00461       
00462     }
00463     
00464   }

Here is the call graph for this function:

int compare_binamp const binamp *  a,
const binamp *  b
 

sort binamp struct from the bigger to the smaller...

Definition at line 153 of file orsa_fft.cc.

Referenced by orsa::psd_max_again_many().

00153                                                        {
00154     
00155     if (( (*a).amp - (*b).amp) < 0)
00156       return  1;
00157     if (( (*a).amp - (*b).amp) > 0)
00158       return -1;
00159     
00160     return 0;
00161   }

double dQ double  y  ) 
 

Definition at line 1319 of file orsa_fft.cc.

References orsa::cos(), orsa::pi, orsa::pisq, and orsa::sin().

01319                       {
01320     
01321     // y==0
01322     if (fabs(y) < 1.0e-12) return (0.0);
01323     
01324     // y==PI
01325     if (fabs(y-pi) < 1.0e-12) return (-3.0/(4*pi));
01326     
01327     // y==-PI
01328     if (fabs(y+pi) < 1.0e-12) return ( 3.0/(4*pi));
01329     
01330     double ysq=y*y;
01331     return (pisq/(y*(pisq-ysq))*(cos(y)+(sin(y)/y)*((3*ysq-pisq)/(pisq-ysq))));
01332     
01333   }

Here is the call graph for this function:

double norm const fftw_complex  z  ) 
 

Definition at line 38 of file orsa_fft.cc.

Referenced by orsa::accurate_peak(), orsa::amph(), orsa::phi_gsl_two(), orsa::psd_max(), orsa::psd_max_again(), and orsa::psd_max_again_many().

00038                                       {
00039     return sqrt( z.re * z.re + z.im * z.im );
00040   }

double norm_sq const fftw_complex  z  ) 
 

Definition at line 42 of file orsa_fft.cc.

Referenced by orsa::phi_amp(), orsa::phi_Hanning_amp(), and orsa::psd_max().

00042                                          {
00043     return ( z.re * z.re + z.im * z.im );
00044   }

fftw_complex phi double  omega,
fftw_complex  in[],
const int  size
 

Discrete Fourier Transform.

Definition at line 49 of file orsa_fft.cc.

References orsa::cos(), orsa::sin(), and orsa::twopi.

Referenced by orsa::accurate_peak(), orsa::amph(), orsa::phi_amp(), and orsa::phi_gsl_two().

00049                                                                     {
00050     
00051     fftw_complex result;
00052     result.re = 0;
00053     result.im = 0;
00054     
00055     double arg,c,s;
00056     int k;
00057     
00058     for(k=0;k<size;k++) {  
00059       arg = twopi*k*omega;
00060       c = cos(arg);
00061       s = sin(arg);
00062       result.re += in[k].re * c + in[k].im * s;
00063       result.im -= in[k].re * s - in[k].im * c;            
00064     }
00065     
00066     double scale = (double)size;
00067     result.re /= scale;
00068     result.im /= scale;
00069     
00070     // OUT HERE!
00071     return result; 
00072   }

Here is the call graph for this function:

double phi_amp double  omega,
fftw_complex  in[],
const int  size
 

amplitude for spectrum, without windowing

Definition at line 106 of file orsa_fft.cc.

References orsa::norm_sq(), and orsa::phi().

Referenced by orsa::phi_gsl().

00106                                                                   {
00107     return sqrt( norm_sq(phi( omega,&in[0],size)) + 
00108                  norm_sq(phi(-omega,&in[0],size)) );
00109   }

Here is the call graph for this function:

double phi_gsl double  x,
void *  params
 

Definition at line 119 of file orsa_fft.cc.

References orsa::params, orsa::phi_amp(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.

00119                                            {
00120     
00121     struct gsl_d1_minimization_parameters * p = (gsl_d1_minimization_parameters *) params;
00122     
00123     // return  (-phi_Hanning_amp(x, p->pointer_points_sequence, p->size));
00124     return  (-phi_amp(x, p->pointer_points_sequence, p->size));
00125     
00126   }

Here is the call graph for this function:

double phi_gsl_two double  x,
void *  params
 

Definition at line 128 of file orsa_fft.cc.

References orsa::norm(), orsa::params, orsa::phi(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.

Referenced by orsa::accurate_peak().

00128                                                {
00129     
00130     struct gsl_d1_minimization_parameters * p = (gsl_d1_minimization_parameters *) params;
00131     
00132     // return  (-phi_Hanning_amp(x, p->pointer_points_sequence, p->size));
00133     return  (-norm(phi(x, p->pointer_points_sequence, p->size)));
00134   }

Here is the call graph for this function:

fftw_complex phi_Hanning double  omega,
fftw_complex  in[],
const int  size
 

Discrete Fourier Transform with Hanning windowing.

Definition at line 77 of file orsa_fft.cc.

References orsa::cos(), orsa::sin(), and orsa::twopi.

Referenced by orsa::phi_Hanning_amp().

00077                                                                             {
00078     
00079     fftw_complex result;
00080     result.re = 0;
00081     result.im = 0;
00082     
00083     double arg,c,s;
00084     int k;
00085     double window_factor;
00086     
00087     for(k=0;k<size;k++) {  
00088       arg = twopi*k*omega;
00089       c = cos(arg);
00090       s = sin(arg);
00091       window_factor = (1-cos(k*twopi/size));
00092       result.re += window_factor * (in[k].re * c + in[k].im * s);
00093       result.im -= window_factor * (in[k].re * s - in[k].im * c);
00094     }
00095     
00096     double scale = (double)size;
00097     result.re /= scale;
00098     result.im /= scale;
00099     
00100     // OUT HERE!
00101     return result; 
00102   }

Here is the call graph for this function:

double phi_Hanning_amp double  omega,
fftw_complex  in[],
const int  size
 

amplitude for spectrum, with Hanning windowing

Definition at line 113 of file orsa_fft.cc.

References orsa::norm_sq(), and orsa::phi_Hanning().

Referenced by orsa::phi_Hanning_gsl().

00113                                                                           {
00114     return sqrt( norm_sq(phi_Hanning( omega,&in[0],size)) + 
00115                  norm_sq(phi_Hanning(-omega,&in[0],size)) );
00116   }

Here is the call graph for this function:

double phi_Hanning_gsl double  x,
void *  params
 

Definition at line 136 of file orsa_fft.cc.

References orsa::params, orsa::phi_Hanning_amp(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.

00136                                                    {
00137     
00138     struct gsl_d1_minimization_parameters * p = (gsl_d1_minimization_parameters *) params;
00139     
00140     return  (-phi_Hanning_amp(x, p->pointer_points_sequence, p->size));
00141     // return  (-phi_amp(x, p->pointer_points_sequence, p->size));
00142     
00143   }

Here is the call graph for this function:

double psd_max const fftw_complex *  transformed_signal,
const int  size
 

Definition at line 375 of file orsa_fft.cc.

References orsa::norm(), and orsa::norm_sq().

00375                                                                          {
00376     
00377     vector<double> psd;
00378     
00379     int k;
00380     
00381     psd.resize(size/2+1);
00382     psd[0] = norm(transformed_signal[0])/size;
00383     for (k=1;k<(size+1)/2;k++)
00384       psd[k] = sqrt(norm_sq(transformed_signal[k])+norm_sq(transformed_signal[size-k]))/size;
00385     if ( (size%2) == 0) // size is even
00386       psd[size/2] = norm(transformed_signal[size/2])/size;  // Nyquist freq.
00387     
00388     /* cerr << " === DUMP PSD ===" << endl;
00389        for (k=0;psd.size();++k) {
00390        printf("PSD[%05i]=%g\n",k,psd[k]);
00391        }
00392     */
00393     
00394     double maxamp=0;
00395     int bin=0;
00396     unsigned int l;
00397     //
00398     double found=false;
00399     if ( (psd[0] > psd[1]) &&
00400          (psd[0] > maxamp) ) {
00401       maxamp = psd[0];
00402       bin    = 0;
00403       found  = true;
00404     }
00405     //
00406     for (l=1;l<(psd.size()-1);l++) {
00407       if ( ( (psd[l] > psd[l-1]) && (psd[l] > psd[l+1]) ) &&
00408            (psd[l] > maxamp) ) {
00409         maxamp = psd[l];
00410         bin    = l;
00411         found  = true;
00412       }
00413     }
00414     
00415     /* for (l=1;l<psd.size();l++) {
00416        if (psd[l] > maxamp) {
00417        maxamp = psd[l];
00418        bin    = l;
00419        }
00420        }
00421     */ 
00422     
00423     // DEBUG: DUMP DATA TO FILE
00424     if (0) {
00425       int filenum=0;
00426       char filename[20];
00427       sprintf(filename,"dump_psd_%02i.dat",filenum);
00428       cerr << "data dump on file " << filename << endl;
00429       FILE *fp = fopen(filename,"w");
00430       if (fp!=0) {
00431         for (l=0;l<psd.size();l++) {
00432           fprintf(fp,"%i  %g\n",l,psd[l]);
00433         }
00434       }
00435       fclose(fp);
00436       filenum++;
00437     }
00438     
00439     // cerr << " psd_max: returning " << ((double)bin/(double)size) << "   amp: " << maxamp << "        bin: " << bin << "  found: " << found << endl;
00440     
00441     if (found) return ((double)bin/(double)size);
00442     else { 
00443       cerr << "\n        *\n        **\n        ***\n        ****\n*************\n**************  WARNING!!! peak don't found... returning (-1)!!\n*************\n        ****\n        ***\n        **\n        *\n" << endl;
00444       return (-1);
00445     }
00446   }

Here is the call graph for this function:

double psd_max_again const fftw_complex *  transformed_signal,
const int  size
 

Definition at line 164 of file orsa_fft.cc.

References orsa::norm().

00164                                                                                {
00165     
00166     vector<double> psd_plus((size-1)/2), psd_minus((size-1)/2);
00167     double Nyquist=0;
00168     
00169     int k;
00170     
00171     const double DC = norm(transformed_signal[0])/size;
00172     for (k=0;k<(size-1)/2;k++)
00173       psd_plus[k]  = norm(transformed_signal[k+1])/size;
00174     for (k=0;k<(size-1)/2;k++)
00175       psd_minus[k] = norm(transformed_signal[size-(k+1)])/size;
00176     if ( (size%2) == 0) // size is even
00177       Nyquist = norm(transformed_signal[size/2])/size;  // Nyquist freq.
00178     
00179     // cerr << "...some test data: " << DC << "  " << psd_plus[0]  << "  " << psd_plus[1]  << "  " << psd_plus[2]  << "  " << psd_plus[3]  << "  " << psd_plus[4] << endl; 
00180     
00181     // DEBUG: DUMP DATA TO FILE
00182     if (0) {
00183       static int filenum=0;
00184       char filename[20];
00185       sprintf(filename,"dump_psd_%02i.dat",filenum);
00186       cerr << "data dump on file " << filename << endl;
00187       FILE *fp = fopen(filename,"w");
00188       if (fp!=0) {
00189         int l;
00190         for (l=psd_minus.size()-1;l>=0;l--) fprintf(fp,"%i  %g\n",-(l+1),psd_minus[l]);
00191         fprintf(fp,"%i  %g\n",0,DC);
00192         for (l=0;l<(int)psd_plus.size(); l++) fprintf(fp,"%i  %g\n", (l+1),psd_plus[l]);
00193         if ( (size%2) == 0) fprintf(fp,"%i  %g\n",size/2,Nyquist);
00194       }
00195       fclose(fp);
00196       filenum++;
00197     }
00198     
00199     // psd.resize(size);
00200     //
00201     /* psd[0] = norm(transformed_signal[0])/size;
00202        for (k=1;k<(size-1)/2;k++)
00203        psd[k] = norm(transformed_signal[k])/size;
00204        for (k=(size/2)+1;k<size;k++)
00205        psd[k] = norm(transformed_signal[size-k])/size;
00206        if ( (size%2) == 0) // size is even
00207        psd[size/2] = norm(transformed_signal[size/2])/size;  // Nyquist freq.
00208     */
00209     
00210     double maxpow=DC;
00211     int bin=0;
00212     
00213     for (k=1;k<(size-3)/2;k++) {
00214       if ( (psd_plus[k] > psd_plus[k-1]) && 
00215            (psd_plus[k] > psd_plus[k+1]) &&
00216            (psd_plus[k] > maxpow) ) {
00217         maxpow = psd_plus[k];
00218         bin    = k+1;
00219       }
00220     }
00221     
00222     if ( (psd_plus[0] > DC) && 
00223          (psd_plus[0] > psd_plus[1]) &&
00224          (psd_plus[0] > maxpow) ) {
00225       maxpow = psd_plus[0];
00226       bin    = 1;
00227     }
00228     
00229     for (k=1;k<(size-3)/2;k++) {
00230       if ( (psd_minus[k] > psd_minus[k-1]) && 
00231            (psd_minus[k] > psd_minus[k+1]) &&
00232            (psd_minus[k] > maxpow) ) {
00233         maxpow = psd_minus[k];
00234         bin    = -(k+1);
00235       }
00236     }
00237     
00238     if ( (psd_minus[0] > DC) && 
00239          (psd_minus[0] > psd_minus[1]) &&
00240          (psd_minus[0] > maxpow) ) {
00241       maxpow = psd_minus[0];
00242       bin    = -1;
00243     }
00244     
00245     if ( (DC > psd_plus[0]) && 
00246          (DC > psd_minus[0]) &&
00247          (DC > maxpow) ) {
00248       maxpow = DC;
00249       bin    = 0;
00250     }
00251     
00252     // TODO: Nyquist...
00253     
00254     // cerr << " psd_max_again: returning " << ((double)bin/(double)size) << "   amp: " << maxpow << "        bin: " << bin << endl;
00255     
00256     // DEBUG: DUMP DATA TO FILE
00257     if (0) {
00258       int filenum=0;
00259       char filename[20];
00260       cerr << "DC:   " << DC << endl;
00261       cerr << "size: " << size << endl;
00262       sprintf(filename,"dump_psd_%02i.dat",filenum);
00263       cerr << "data dump on file " << filename << endl;
00264       FILE *fp = fopen(filename,"w");
00265       if (fp!=0) {
00266         int l;
00267         for (l=psd_minus.size()-1;l>=0;l--) fprintf(fp,"%i  %g\n",-(l+1),psd_minus[l]);
00268         fprintf(fp,"%i  %g\n",0,DC);
00269         for (l=0;l<(int)psd_plus.size(); l++) fprintf(fp,"%i  %g\n", (l+1),psd_plus[l]);
00270         if ( (size%2) == 0) fprintf(fp,"%i  %g\n",size/2,Nyquist);
00271       }
00272       fclose(fp);
00273       filenum++;
00274     }
00275     
00276     if (maxpow>0) return ((double)bin/(double)size);
00277     /* else { 
00278        cerr << "\n        *\n        **\n        ***\n        ****\n*************\n**************  WARNING!!! peak don't found... returning (-1)!!\n*************\n        ****\n        ***\n        **\n        *\n" << endl;
00279        return (-1);
00280        }
00281     */
00282     
00283     return ((double)bin/(double)size);
00284   }

Here is the call graph for this function:

void psd_max_again_many const fftw_complex *  transformed_signal,
const int  size,
vector< double > &  candidate,
const unsigned int  nfreq
 

Definition at line 286 of file orsa_fft.cc.

References orsa::compare_binamp(), and orsa::norm().

00286                                                                                                                                        {
00287     
00288     vector<double> psd_plus((size-1)/2), psd_minus((size-1)/2);
00289     double Nyquist=0;
00290     
00291     int k;
00292     
00293     const double DC = norm(transformed_signal[0])/size;
00294     for (k=0;k<(size-1)/2;k++)
00295       psd_plus[k]  = norm(transformed_signal[k+1])/size;
00296     for (k=0;k<(size-1)/2;k++)
00297       psd_minus[k] = norm(transformed_signal[size-(k+1)])/size;
00298     if ( (size%2) == 0) // size is even
00299       Nyquist = norm(transformed_signal[size/2])/size;  // Nyquist freq.
00300     
00301     vector<binamp> vec_binamp;
00302     binamp ba;
00303     
00304     for (k=1;k<(size-3)/2;k++) {
00305       if ( (psd_plus[k] > psd_plus[k-1]) && 
00306            (psd_plus[k] > psd_plus[k+1]) ) {
00307         // candidate.push_back((double)(k+1)/(double)size);
00308         ba.bin = k+1;
00309         ba.amp = psd_plus[k];
00310         vec_binamp.push_back(ba);
00311       }
00312     }
00313     
00314     if ( (psd_plus[0] > DC) && 
00315          (psd_plus[0] > psd_plus[1]) ) {
00316       // candidate.push_back(1.0/(double)size);
00317       ba.bin = 1;
00318       ba.amp = psd_plus[0];
00319       vec_binamp.push_back(ba); 
00320     }
00321     
00322     for (k=1;k<(size-3)/2;k++) {
00323       if ( (psd_minus[k] > psd_minus[k-1]) && 
00324            (psd_minus[k] > psd_minus[k+1]) ) {
00325         // candidate.push_back((double)-(k+1)/(double)size);
00326         ba.bin = -(k+1);
00327         ba.amp = psd_minus[k];
00328         vec_binamp.push_back(ba);  
00329       }
00330     }
00331     
00332     if ( (psd_minus[0] > DC) && 
00333          (psd_minus[0] > psd_minus[1]) ) {
00334       // candidate.push_back(-1.0/(double)size);
00335       ba.bin = -1;
00336       ba.amp = psd_minus[0];
00337       vec_binamp.push_back(ba);
00338     }
00339     
00340     if ( (DC > psd_plus[0]) && 
00341          (DC > psd_minus[0]) ) {
00342       // candidate.push_back(0.0);
00343       ba.bin = 0;
00344       ba.amp = DC;
00345       vec_binamp.push_back(ba);
00346     }
00347     
00348     
00349     unsigned int bsize = vec_binamp.size();
00350     binamp *ptr_binamp;
00351     ptr_binamp = (binamp *) malloc(bsize*sizeof(binamp));
00352     { 
00353       unsigned int k;
00354       for (k=0;k<bsize;k++) {
00355         ptr_binamp[k].bin = vec_binamp[k].bin;
00356         ptr_binamp[k].amp = vec_binamp[k].amp;
00357       }
00358     }
00359     gsl_heapsort(ptr_binamp,bsize,sizeof(binamp),(gsl_comparison_fn_t)compare_binamp);
00360     
00361     candidate.clear();
00362     
00363     {
00364       unsigned int k;
00365       for (k=0;k<nfreq;k++) {
00366         candidate.push_back((double)ptr_binamp[k].bin/(double)size);
00367       } 
00368     }
00369     
00370     free(ptr_binamp);
00371   }

Here is the call graph for this function:


Generated on Tue Jan 11 15:26:48 2005 for liborsa by  doxygen 1.4.0