src/time_frequency.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2008 by Matthias Ihrke   *
 *   mihrke@uni-goettingen.de   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program 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 General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "time_frequency.h"
#include "linalg.h"

Spectrogram* spectrogram( const Array *sig, Spectrogram *spectgram, 
                                  OptArgList *optargs ){
  double x;
  void *ptr;
  double sample_frequency;
  WindowFunction winfct; 
  int winlength;
  double winparam;
  int N_freq;
  int N_time;
  double corner_freqs[2];
  int *timepoints;

  /* check input */
  bool isvec;
  vector_CHECK( isvec, sig );
  if( !isvec ) return NULL;

  int n=sig->size[0];
  /* defaults */
  sample_frequency = 500.0;
  winfct = window_hanning;
  winlength =  MAX( SQR( sqrt(next_pow2( n ))-3 ), 5 );
  winparam=0.0;
  N_freq = winlength*4;
  N_time = n;
  corner_freqs[0] = 0.0;
  corner_freqs[1] = 250.0;
  timepoints = NULL;
  bool N_time_set=FALSE;

  /* get params */
  optarg_PARSE_SCALAR( optargs, "sample_frequency", sample_frequency, double, x );
  optarg_PARSE_PTR( optargs, "winfct", winfct, WindowFunction, ptr );
  optarg_PARSE_SCALAR( optargs, "winlength", winlength, int, x );
  optarg_PARSE_SCALAR( optargs, "winparam", winparam, double, x );
  optarg_PARSE_SCALAR( optargs, "N_freq", N_freq, int, x );
  if( optarglist_has_key( optargs, "N_time" ) ){
     x = optarglist_scalar_by_key( optargs, "N_time" );
     if( !isnan( x ) ){
        N_time=(int)x;
        N_time_set=TRUE;
     }
  }
  if( optarglist_has_key( optargs, "corner_freqs" ) ){
     ptr = optarglist_ptr_by_key( optargs, "corner_freqs" );
     if( ptr ){
        corner_freqs[0] = ((double*)ptr)[0];
        corner_freqs[1] = ((double*)ptr)[1];
     }
  }
  if( optarglist_has_key( optargs, "timepoints" ) ){
     ptr = optarglist_ptr_by_key( optargs, "timepoints" );
     if( ptr && N_time_set ){
        timepoints = (int*)ptr;
     }
  }

  dprintf("got sfreq=%f, winfct=%p, winlength=%i, winparam=%f, N_f=%i,"
             "N_t=%i, cf=(%f,%f), timepoints=%p\n",
             sample_frequency, winfct, winlength, winparam, N_freq, N_time, 
             corner_freqs[0], corner_freqs[1], timepoints );

  Array *window=winfct( winlength, winparam );

  spectgram = spectrogram_stft( sig, sample_frequency, 
                                          window, N_freq, N_time, corner_freqs, 
                                          timepoints, spectgram );
  return spectgram;
}


Spectrogram* spectrogram_stft(const Array* sig, double sampling_rate,
                                        const Array *Window, int N_freq, int N_time, 
                                        double corner_freqs[2], 
                                        int *timepoints, Spectrogram *spectgram){

  int            Nfft, timepoint, frequency, time;
  int            taumin, taumax, tau;
  int            half_Window_Length;
  double        *wind_sig_complex;      /* windowed signal */
  double         normh;
  double         inter;
  int            corner_freqs_idx[2]; /* nearest sample in discretized setting */
  int            N_freq_spect;    /* number of frequencies in spectrgram */
  int            Window_Length;

  /* check for 1D and DOUBLE */
  bool isvec;
  vector_CHECK( isvec, sig );
  if( !isvec ) return NULL;
  vector_CHECK( isvec, Window );
  if( !isvec ) return NULL;
  
  int n=sig->size[0];
  Window_Length=Window->size[0];

  /* normalizing to [0, 100] */
  half_Window_Length = (Window_Length - 1) / 2;
  inter = 0.0;
  for (frequency = 0; frequency <Window_Length; frequency++){
     inter += SQR (vec_IDX(Window,frequency));
  }
  normh = sqrt (inter);

  /* FFT-settings */
  Nfft = pow( 2, next_pow2( N_freq ));    /* num of samples for FFT */
  if( corner_freqs[1]<=corner_freqs[0] || 
        corner_freqs[0]<0 || corner_freqs[1]>sampling_rate/2.0 ){
     errprintf("corner_freqs funny: (%f,%f)\n", 
                  corner_freqs[0], corner_freqs[1] );
     return NULL;
  }
  corner_freqs_idx[0] = (int)round( (corner_freqs[0]/(sampling_rate/2.0)) * (N_freq-1) );
  corner_freqs_idx[1] = (int)round( (corner_freqs[1]/(sampling_rate/2.0)) * (N_freq-1) );
  dprintf(" Corner-freqs=(%.2f,%.2f) in idx=(%i,%i)\n", corner_freqs[0], corner_freqs[1],
             corner_freqs_idx[0], corner_freqs_idx[1] );
  dprintf(" Returning actual frequencies: (%.3f, %.3f)\n", 
             (corner_freqs_idx[0]/(double)N_freq)*sampling_rate/2.0, 
             (corner_freqs_idx[1]/(double)N_freq)*sampling_rate/2.0 ); 
  N_freq_spect = corner_freqs_idx[1]-corner_freqs_idx[0];
  

  /*                memory allocation for the windowed signal           */
  wind_sig_complex = (double *) calloc (2*Nfft, sizeof (double));

  if( spectgram==NULL ){
     spectgram = spectrogram_init( N_freq_spect, N_time );
  }
  spectgram->low_corner_freq = corner_freqs[0];
  spectgram->up_corner_freq  = corner_freqs[1];

  /* computation of the fft for the current windowed signal */
  for (timepoint = 0; timepoint < N_time; timepoint++) {
     /* initialization of the intermediary vector */
     for (frequency = 0; frequency < 2*Nfft; frequency++){
        wind_sig_complex[frequency] = 0.0;
     }

     /* current time to compute the stft */
     if( timepoints ){
        time = timepoints[timepoint];
     } else {
        time = timepoint*(int)(n/N_time);
     }
     /* dprintf( "Spectrum at sample '%i'\n", time ); */

     /* the signal is multipied by the window between the instants
         time-taumin and time+taumax 
        when the window is wider than the number of desired frequencies (tfr.N_freq),
         the range is limited to tfr.N_freq */
     taumin = MIN (N_freq / 2, half_Window_Length);
     taumin = MIN (taumin, time);

     taumax = MIN ((N_freq / 2 - 1), half_Window_Length);
     taumax = MIN (taumax, (n - time - 1));

     /* dprintf(" tau = (%i, %i)\n", taumin, taumax); */
     /* The signal is windowed around the current time */
     frequency=0;
     for (tau = -taumin; tau <= taumax; tau++){
        wind_sig_complex[2*frequency  ] = vec_IDX(sig,time + tau) *
          vec_IDX( Window, half_Window_Length + tau ) / normh;
        wind_sig_complex[2*frequency+1] = 0.0;
        frequency++;
     }
     /* fft of the windowed signal */
     fft( wind_sig_complex, Nfft, 1 );  
     /* gsl_fft_complex_radix2_forward ( wind_sig_complex, 1, Nfft );  */

     /* the first half of the fft is put in the stft matrix  */
     for (frequency = corner_freqs_idx[0]; frequency < corner_freqs_idx[1]; frequency++){
        spectgram->sgram[timepoint][frequency].re = wind_sig_complex[2*frequency  ];
        spectgram->sgram[timepoint][frequency].im = wind_sig_complex[2*frequency+1];
     }
  } /* for t */

 /*--------------------------------------------------------------------*/
  free (wind_sig_complex);

  return spectgram;
}


Array* spectrogram_powerspectrum( const Spectrogram *s ){
  int i,j;
  Array *p = array_new2( DOUBLE, 2, s->N_time, s->N_freq );
  
  for( i=0; i<s->N_time; i++ ){
     for( j=0; j<s->N_freq; j++ ){
        mat_IDX( p, i,j) = SQR( complex_abs( s->sgram[i][j] ) );
     }
  }
  return p;
}

Spectrogram* spectrogram_init( int N_freq, int N_time ){
  Spectrogram *s;
  int i;
  
  dprintf("N_freq=%i, N_time=%i\n", N_freq, N_time );
  s = (Spectrogram*) malloc( sizeof(Spectrogram) );
  s->samplingrate = 500.0;
  s->N_freq = N_freq;
  s->N_time = N_time;
  s->low_corner_freq=0;
  s->up_corner_freq=0;
  s->sgram = (Complex**) malloc( N_time*sizeof(Complex*) );
  for( i=0; i<N_time; i++ ){
     s->sgram[i] = (Complex*) malloc( N_freq*sizeof(Complex) );
  }
  return s;
}

void spectrogram_free( Spectrogram *s ){
  int i;

  for( i=0; i<s->N_time;  i++ ){
     free( s->sgram[i] );
  }
  free( s->sgram );
  free( s );
}

/*----------------- WINDOW-FUNCTIONS ------------------ */

Array* window_dirichlet(  int n, double noparam ){
  int i;

  if( !ISODD( n ) ){
     errprintf("window size must be ODD! Adding 1, n=%i!\n", n+1);
     n++;
  }
  Array *window = array_new2( DOUBLE, 1, n );
  for( i=0; i<n; i++ ){
     vec_IDX( window,i ) = 1.0;
  }

  return window;
}

Array* window_gaussian( int n, double sigma ){
  int i;

  if( !ISODD( n ) ){
     errprintf("window size must be ODD! Adding 1, n=%i!\n", n+1);
     n++;
  }

  Array *window = array_new2( DOUBLE, 1, n );
  for( i=0; i<n; i++ ){
     vec_IDX(window,i) = exp(-0.5*SQR( (i-(n-1)/2.0)/(double)(sigma*(n-1)/2.0)) );
  }

  return window;
}

Array* window_hamming( int n, double noparam ){
  int i;

  if( !ISODD( n ) ){
     errprintf("window size must be ODD! Adding 1, n=%i!\n", n+1);
     n++;
  }

  Array *window = array_new2( DOUBLE, 1, n );
  for( i=0; i<n; i++ ){
     vec_IDX(window,i) =  0.53836 - 0.46164*cos( (2* PI *i)/(double)(n-1) );
  }

  return window;
}
Array* window_hanning( int n, double noparam ){
  int i;

  if( !ISODD( n ) ){
     errprintf("window size must be ODD! Adding 1, n=%i!\n", n+1);
     n++;
  }
  Array *window = array_new2( DOUBLE, 1, n );
  for( i=0; i<n; i++ ){
     vec_IDX( window,i) =  0.5*(1 - cos( (2* PI *i)/(double)(n-1) ));
  }

  return window;
}


Array* window_kaiser( int n, double alpha ){
  int i;
  double bessel_alpha;

  Array *window = array_new2( DOUBLE, 1, n );

  bessel_alpha=gsl_sf_bessel_I0( alpha );
  for( i=0; i<n-1; i++ ){
     vec_IDX(window,i) = gsl_sf_bessel_I0( alpha*sqrt(1 - SQR( (2*i)/(double)(n-1) - 1  ) ) ) / bessel_alpha;
  }
  vec_IDX(window,n-1)=0;

  return window;
}