• API Main Page
  • Documentation
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

src/distances.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  *   Copyright (C) 2008/2009 by Matthias Ihrke   *
00003  *   mihrke@uni-goettingen.de   *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019  ***************************************************************************/
00020 
00021 #include "distances.h"
00022 #include "optarg.h"
00023 #include <gsl/gsl_statistics.h>
00024 
00036 Array* matrix_distmatrix( VectorDistanceFunction f, 
00037                                   const Array *X, Array *D, 
00038                                   OptArgList *optargs ){
00039   int i,j;
00040   int idx;
00041   int N,p;
00042   void *ptr;
00043   ProgressBarFunction progress=NULL; 
00044   if( optarglist_has_key( optargs, "progress" ) ){
00045      ptr = optarglist_ptr_by_key( optargs, "progress" );
00046      if( ptr ) progress = (ProgressBarFunction)ptr;
00047   }
00048   
00049   bool ismatrix;
00050   matrix_CHECK( ismatrix, X );
00051   if( !ismatrix ) return NULL;
00052   N = X->size[0];
00053   p = X->size[1];
00054   if( D==ALLOC_IN_FCT ){
00055      D = array_new2( DOUBLE, 2, N, N );
00056   } else {
00057      matrix_CHECK( ismatrix, D );
00058      if( !ismatrix ) return NULL;
00059      if( D->size[0]<N || D->size[1]<N ){
00060         errprintf("Matrix not large enough, need %i x %i, have %i x %i\n",
00061                      N, N, D->size[0], D->size[1] );
00062      }
00063   }
00064 
00065   if( progress ){
00066      progress( PROGRESSBAR_INIT, N*(N-1)/2 );
00067   }
00068   
00069   idx=1;
00070   for( i=0; i<N; i++ ){
00071      for( j=i+1; j<N; j++ ){
00072         if( progress ){
00073           progress( PROGRESSBAR_CONTINUE_LONG, idx );
00074         }
00075         mat_IDX( D, i, j) = f( (double*)array_INDEXMEM2(X, i, 0), 
00076                                       (double*)array_INDEXMEM2(X, j, 0), p, optargs );
00077         mat_IDX(D, j,i) = mat_IDX( D, i,j);
00078         idx++;
00079      }
00080   }
00081   if( progress ){
00082      progress( PROGRESSBAR_FINISH, 0 );
00083   }
00084 
00085   return D;
00086 }
00087 
00098 Array* distmatrix_signaldist( VectorDistanceFunction f, const Array *s1, 
00099                                         const Array *s2, Array *out,    OptArgList *optargs ){
00100   bool ismatrix;
00101   int i,j;
00102 
00103   if( s1->dtype!=DOUBLE || s2->dtype!=DOUBLE ){
00104      errprintf("Input signals must be double-arrays\n");
00105      return NULL;
00106   }
00107   if( s1->ndim>2 || s2->ndim>2 || s1->ndim != s2->ndim ){
00108      warnprintf("Input signals should be 1D or 2D... continue at your own risk (%i,%i)\n",
00109                     s1->ndim, s2->ndim );
00110   }
00111 
00112   /* thin matrix-wrapper (in case of 1D data */
00113   Array *s1m, *s2m;
00114   s1m = array_fromptr2( DOUBLE, 2, s1->data, s1->size[0], (s1->ndim>1)?(s1->size[1]):1 );
00115   s2m = array_fromptr2( DOUBLE, 2, s2->data, s2->size[0], (s2->ndim>1)?(s2->size[1]):1 );
00116 
00117   if( s1m->size[1] != s2m->size[1] ){
00118      errprintf("Signals must have the same dimension 2\n");
00119      array_free( s1m );
00120      array_free( s2m );
00121      return NULL; 
00122   }
00123 
00124   int N1 = s1m->size[0];
00125   int N2 = s2m->size[0];
00126   int p = s1m->size[1];
00127   if( out ){
00128      matrix_CHECK( ismatrix, out );
00129      if( !ismatrix ) return NULL;
00130      if( out->size[0] < N1 || out->size[1] < N2 ){
00131         errprintf("output matrix is not %ix%i: %i,%i\n",
00132                      N1,N2,out->size[0],out->size[1]);
00133         array_free( s1m );
00134         array_free( s2m );
00135         return NULL;
00136      }
00137   } else {
00138      out = array_new2( DOUBLE, 2, N1, N2 );
00139   }
00140 
00141   /* --------------------- computation ------------------------*/
00142   for( i=0; i<N1; i++ ){
00143      for( j=0; j<N2; j++ ){
00144         mat_IDX( out, i, j ) = f( (double*)array_INDEXMEM2( s1m, i, 0 ),
00145                                           (double*)array_INDEXMEM2( s2m, j, 0 ), p, optargs );
00146      }
00147   }
00148   /* --------------------- /computation ------------------------*/
00149 
00150   array_free( s1m );
00151   array_free( s2m );
00152   return out; 
00153 }
00154 
00160 double vectordist_euclidean( const double *x1, const double *x2, int p, OptArgList *optargs ){
00161   double d;
00162   int i;
00163 
00164   d=0.0;
00165   for( i=0; i<p; i++ ){
00166      d += SQR( x1[i]-x2[i] );
00167   }
00168   d = sqrt( d );
00169 
00170   return d;
00171 }
00182 double   vectordist_euclidean_normalized( const double *x1, const double *x2, int p,
00183                                                         OptArgList *optargs ){
00184   int i;
00185   double d;
00186   double avg1, avg2,
00187      rms1=0, rms2=0;
00188   double norm1,norm2; /* normalized signal at time t */
00189      
00190   avg1 = gsl_stats_mean(x1, 1, p);
00191   avg2 = gsl_stats_mean(x2, 1, p);
00192   
00193   for(i=0; i<p; i++){
00194      rms1 += SQR( x1[i] );
00195      rms2 += SQR( x2[i] );
00196   }
00197   
00198   rms1 = sqrt(rms1/(double)p);
00199   rms2 = sqrt(rms2/(double)p);
00200 
00201   d=0.0;
00202   for( i=0; i<p; i++ ){
00203      norm1 = (x1[i]-avg1)/rms1;
00204      norm2 = (x2[i]-avg2)/rms2;
00205 
00206      d += SQR( norm1-norm2 );
00207   }
00208   d = sqrt( d );
00209 
00210   return d;
00211 }
00221 double   vectordist_dtw( const double *x1, const double *x2, int p, 
00222                                  OptArgList *optargs ){
00223   double d;
00224   double **D;
00225   PointwiseDistanceFunction f;
00226   void *tmp;
00227   WarpPath *P;
00228   double diag1[2], diag2[2], path[2];
00229   int i;
00230 
00231   f = signaldist_euclidean;
00232   if( optarglist_has_key( optargs, "distfct" ) ){
00233      tmp = optarglist_ptr_by_key( optargs, "distfct" );
00234      if( tmp )
00235         f = (PointwiseDistanceFunction) tmp;
00236   }
00237 
00238   D = f( (double*)x1, p, (double*)x2, p, NULL, optargs );
00239   dtw_cumulate_matrix( D, p, p, NULL );
00240   P = dtw_backtrack( (const double**)D, p, p, NULL );
00241 
00242   d = 0;
00243   diag1[0] = 0;   
00244   diag1[1] = 0;
00245   diag2[0] = p;
00246   diag2[0] = p; /* line through (0,0), (p,p) */
00247   for( i=0; i<P->n; i++ ){
00248      path[0] = (double)P->t1[i];
00249      path[1] = (double)P->t2[i];
00250      d += dist_point_line( path, diag1, diag2 );
00251   }
00252 
00253   dblpp_free( D, p );
00254   free_warppath( P );
00255 
00256   return d;
00257 }
00258 
00269 double** eeg_distmatrix( EEG *eeg, VectorDistanceFunction f, double **d, OptArgList *optargs ){
00270 #ifdef FIXEEG
00271   int i,j;
00272   int c;
00273   void *tmp;
00274   ProgressBarFunction progress=NULL;
00275 
00276   if( d==ALLOC_IN_FCT ){
00277      warnprintf(" allocating matrix in fct\n");
00278      d = dblpp_init( eeg->ntrials, eeg->ntrials );
00279   }
00280 
00281   if( optarglist_has_key( optargs, "progress" ) ){
00282      tmp = optarglist_ptr_by_key( optargs, "progress" );
00283      if( tmp ) progress = (ProgressBarFunction) tmp;
00284   }
00285   dprintf("progress=%p\n", progress );
00286 
00287   if( progress ){
00288      progress( PROGRESSBAR_INIT, eeg->ntrials );
00289   }
00290 
00291   for( i=0; i<eeg->ntrials; i++ ){
00292      d[i][i] = 0.0;  
00293      if( progress ){
00294         progress( PROGRESSBAR_CONTINUE_LONG, i );
00295      }
00296      for( j=i+1; j<eeg->ntrials; j++ ){
00297         d[i][j] = 0.0;
00298         for( c=0; c<eeg->nbchan; c++ ){
00299           if( progress ){
00300              progress( PROGRESSBAR_CONTINUE_SHORT, 0 );
00301           }
00302           d[i][j] += f( eeg->data[c][i], 
00303                              eeg->data[c][j],
00304                              eeg->n, optargs );
00305           if( isnan( d[i][j] ) ){
00306              errprintf("D[%i][%i] is nan\n", i, j);
00307           }
00308         }
00309         d[i][j] /= (double)eeg->nbchan;
00310         d[j][i] = d[i][j];
00311      }
00312   }
00313   if( progress ){
00314      progress( PROGRESSBAR_FINISH, 0 );
00315   }
00316   return d;
00317 #endif
00318 }
00319 
00320 
00330 double   dist_point_line(double *p, double *x, double *y){
00331   double det=0.0;
00332   double d=0.0;
00333   det = ((y[0]-x[0])*(p[1]-x[1])) - ((y[1]-x[1])*(p[0]-x[0]));
00334   d = (double)ABS( det ) / (double)sqrt( SQR(y[0]-x[0]) + SQR(y[1]-x[1]) );
00335 
00336   return d;
00337 }
00338 
00351 double   pathdist_euclidean_dt(WarpPath *p1, WarpPath *p2){
00352 #ifdef LIBEEGTOOLS_FIXME
00353   int **I;
00354   int n,m,
00355      i,j;
00356   double dist;
00357   double **d;
00358   
00359   n = p1->t1[p1->n-1]+1;
00360   m = p1->t2[p1->n-1]+1;
00361 
00362   dprintf("calculate distance with (%i x %i) points\n", n, m );
00363 
00364   I = (int**)malloc( n*sizeof(int*) );
00365   for( i=0; i<n; i++){
00366      I[i] = (int*)malloc( m*sizeof(int) );
00367      for( j=0; j<m; j++ ){
00368         I[i][j] = 0;
00369      }
00370   }
00371 
00372   for( i=0; i<p1->n; i++ ){
00373      //dprintf(" (%i , %i)\n", p1->t1[i], p1->t2[i] );
00374      I[p1->t1[i]][p1->t2[i]]=1;
00375   }
00376   
00377   dprintf(" init done \n");
00378      
00379   d = disttransform_deadreckoning( I, n, m, ALLOC_IN_FCT );
00380   dblpp_divide_scalar( d, n, m, dblpp_max( (const double**)d, n, m, NULL, NULL ) );
00381   dprintf(" dt done \n");
00382 
00383   dist = 0.0;
00384   for( i=0; i<p2->n; i++ ){
00385      //  dprintf(" (%i , %i)\n", p2->t1[i], p2->t2[i] );
00386      if(  p2->t1[i]>=n ||  p2->t2[i]>=n ){
00387         warnprintf("points in P2 do not fall in P1,  (%i , %i) -> ignored\n", p2->t1[i], p2->t2[i] );
00388      } else {
00389         dist += d[p2->t1[i]][p2->t2[i]];
00390      }
00391   }
00392   
00393   //  dist /= (double)n*m;
00394   dprintf(" dist=%f compute done \n", dist);
00395  
00396   for( i=0; i<n; i++ ){
00397      free( I[i] );
00398   }
00399   free(I); 
00400   //dblpp_free( d, m );
00401 
00402   return dist;
00403 #endif
00404 }
00405 
00406 /**************************************************************
00407 GOING TO BE OBSOLETE 
00408 **************************************************************/
00412 double** signaldist_euclidean( double *s1, int n1, double *s2, int n2,
00413                                          double **d, OptArgList *optargs ){
00414   int i,j;
00415   if( d==ALLOC_IN_FCT ){
00416      d=dblpp_init( n1, n2 );
00417   }
00418 
00419   for( i=0; i<n1; i++){
00420      for( j=0; j<n2; j++ ){
00421         d[i][j] = SQR( fabs( s1[i] - s2[j] ) );
00422      }
00423   }
00424   return d;
00425 }
00426 
00427 
00445 double** signaldist_euclidean_derivative( double *s1, int n1, double *s2, int n2, double **d, OptArgList *optargs ){
00446   int i, j, k;
00447   double avg1, avg2,
00448      rms1=0, rms2=0;
00449   double norm1,norm2, /* normalized signal at time t */
00450      norm1p,norm2p;    /* normalized signal at time t-1 */
00451   double theta1,theta2;
00452   double x;
00453 
00454   theta1=1.0;
00455   theta2=1.0;
00456 
00457   if( optarglist_has_key( optargs, "theta1" ) ){
00458      x = optarglist_scalar_by_key( optargs, "theta1" );
00459      if( !isnan( x ) )
00460         theta1=x;
00461   }
00462   if( optarglist_has_key( optargs, "theta2" ) ){
00463      x = optarglist_scalar_by_key( optargs, "theta2" );
00464      if( !isnan( x ) )
00465         theta2=x;
00466   }
00467   dprintf("Using theta=(%f,%f)\n", theta1, theta2 );
00468 
00469 
00470   if( d==ALLOC_IN_FCT ){
00471      d=dblpp_init( n1, n2 );
00472   }
00473 
00474   avg1 = gsl_stats_mean(s1, 1, n1);
00475   avg2 = gsl_stats_mean(s2, 1, n2);
00476   for(i=0; i<MAX(n1, n2); i++){
00477      if( i<n1 )
00478         rms1 += SQR( s1[i] );
00479      if( i<n2 )
00480         rms2 += SQR( s2[i] );
00481   }
00482 
00483   rms1 = sqrt(rms1/(double)n1);
00484   rms2 = sqrt(rms2/(double)n2);
00485 
00486   for( j=0; j<n1; j++ ){
00487      for( k=0; k<n2; k++ ){
00488       norm1 = (s1[j]-avg1)/rms1;
00489       norm2 = (s2[k]-avg2)/rms2;
00490       (j==0) ? (norm1p = 0) : (norm1p = ((s1[j-1]-avg1)/rms1));
00491       (k==0) ? (norm2p = 0) : (norm2p = ((s2[k-1]-avg2)/rms2));
00492         d[j][k] = theta1*fabs(norm1 - norm2) + theta2*fabs( (norm1-norm1p) - (norm2-norm2p) ); /* (1) */
00493      }
00494   }
00495   return d;
00496 }
00497 
00498 
00523 double** signaldist_stft( double *s1, int n1, double *s2, int n2, double **d, OptArgList *optargs ){
00524 #ifdef LIBEEGTOOLS_FIXME
00525   int i, j;
00526   Spectrogram *sp1, *sp2;
00527   double *window;
00528   double *d1, *d2,
00529      mean1, mean2;
00530   double x;
00531   void *ptr;
00532 
00533   double sample_frequency;
00534   WindowFunction winfct; 
00535   int winlength;
00536   int N_freq;
00537   int N_time;
00538   double corner_freqs[2];
00539 
00540   /* defaults */
00541   sample_frequency = 500.0;
00542   winfct = window_hanning;
00543   winlength =  MAX( SQR( sqrt(next_pow2( n1 ))-3 ), 5 );
00544   N_freq = winlength*4;
00545   N_time = n1;
00546   corner_freqs[0] = 0.0;
00547   corner_freqs[1] = 250.0;
00548   
00549   /* get params */
00550   if( optarglist_has_key( optargs, "sample_frequency" ) ){
00551      x = optarglist_scalar_by_key( optargs, "sample_frequency" );
00552      if( !isnan( x ) ) sample_frequency=x;
00553   }
00554   if( optarglist_has_key( optargs, "winfct" ) ){
00555      ptr = optarglist_ptr_by_key( optargs, "winfct" );
00556      if( ptr ) winfct = (WindowFunction)ptr;
00557   }
00558   if( optarglist_has_key( optargs, "winlength" ) ){
00559      x = optarglist_scalar_by_key( optargs, "winlength" );
00560      if( !isnan( x ) ) winlength=(int)x;
00561   }
00562   if( optarglist_has_key( optargs, "N_freq" ) ){
00563      x = optarglist_scalar_by_key( optargs, "N_freq" );
00564      if( !isnan( x ) ) N_freq=(int)x;
00565   }
00566   if( optarglist_has_key( optargs, "N_time" ) ){
00567      x = optarglist_scalar_by_key( optargs, "N_time" );
00568      if( !isnan( x ) ) N_time=(int)x;
00569   }
00570   if( optarglist_has_key( optargs, "corner_freqs" ) ){
00571      ptr = optarglist_ptr_by_key( optargs, "corner_freqs" );
00572      if( ptr ){
00573         corner_freqs[0] = ((double*)ptr)[0];
00574         corner_freqs[1] = ((double*)ptr)[1];
00575      }
00576   }
00577 
00578   dprintf("got sfreq=%f, winfct=%p, winlength=%i, N_f=%i, N_t=%i, cf=(%f,%f)\n",
00579              sample_frequency, winfct, winlength, N_freq, N_time, corner_freqs[0], corner_freqs[1] );
00580 
00581   if( d==ALLOC_IN_FCT ){
00582      warnprintf("allocating matrix in fct\n");
00583      d=dblpp_init( n1, n2 );
00584   }
00585 
00586   /* remove mean from the signals */
00587   dprintf("Normalization\n");
00588   d1 = dblp_init( NULL, n1, 0.0 );
00589   d2 = dblp_init( NULL, n2, 0.0 );
00590   mean1 = gsl_stats_mean( s1, 1, n1 );
00591   mean2 = gsl_stats_mean( s2, 1, n2 );
00592   for( i=0; i<MAX(n1,n2); i++ ){      
00593      if( n1<i )
00594         d1[i] = s1[i]-mean1;
00595      if( n2<i )
00596         d2[i] = s2[i]-mean2;
00597   }
00598 
00599   window = winfct(ALLOC_IN_FCT, winlength);
00600   sp1 = spectrogram_stft( d1, n1, sample_frequency,
00601                                   window, winlength, N_freq, N_time, 
00602                                   corner_freqs, NULL, ALLOC_IN_FCT );
00603   sp2 = spectrogram_stft( d2, n2, sample_frequency,
00604                                   window, winlength, N_freq, N_time,
00605                                   corner_freqs, NULL, ALLOC_IN_FCT );
00606   if( !sp1 || !sp2 ){
00607      errprintf("Spectrogram broken\n");
00608      return NULL;
00609   }
00610   spectrogram_compute_powerspectrum( sp1 );
00611   spectrogram_compute_powerspectrum( sp2 );
00612 
00613   dprintf(" Compute distances (%i vectors)\n", sp1->N_freq);
00614   for( i=0; i<N_time; i++ ){
00615      for( j=0; j<N_time; j++ ){
00616         d[i][j] = dblp_euclidean_distance( sp1->powerspect[i], sp2->powerspect[j], sp1->N_freq );
00617      }
00618   }
00619   dprintf("\nDone\n");
00620 
00621   /* free */
00622   spectrogram_free( sp1 );
00623   spectrogram_free( sp2 );
00624   free( d1 );
00625   free( d2 );
00626   free( window );
00627   
00628   return d;
00629 #endif
00630 }
00631 
00642 double** vectordist_distmatrix( VectorDistanceFunction f, 
00643                                           const double **X, int n, int p, double **D, 
00644                                           ProgressBarFunction progress,
00645                                           OptArgList *optargs ){
00646   int i,j;
00647   int idx;
00648   
00649   dprintf("data dim is (%i,%i)\n", n, p);
00650   if( D==ALLOC_IN_FCT ){
00651      warnprintf(" allocating matrix in fct\n");
00652      D = dblpp_init( n, n );
00653   }
00654   if( progress ){
00655      progress( PROGRESSBAR_INIT, n*(n-1)/2 );
00656   }
00657   
00658   idx=1;
00659   for( i=0; i<n; i++ ){
00660      for( j=i+1; j<n; j++ ){
00661         if( progress ){
00662           progress( PROGRESSBAR_CONTINUE_LONG, idx );
00663         }
00664         D[i][j] = f( (double*)X[i], (double*)X[j], p, optargs );
00665         D[j][i] = D[i][j];
00666 
00667 
00668         if( isnan( D[j][i] ) ){
00669           errprintf("D[%i][%i] is nan\n", j, i);
00670         }
00671         idx++;
00672      }
00673   }
00674   if( progress ){
00675      progress( PROGRESSBAR_FINISH, 0 );
00676   }
00677 
00678   return D;
00679 }

Generated on Fri Jun 25 2010 14:10:17 for libeegtools by  doxygen 1.7.0