00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
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
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
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;
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;
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
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
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
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
00401
00402 return dist;
00403 #endif
00404 }
00405
00406
00407
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,
00450 norm1p,norm2p;
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) );
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
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
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
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
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 }
