src/filter.c

Go to the documentation of this file.

/* **************************************************************************
 *   Copyright (C) 2009 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 "filter.h"
#include "eeg.h"
#include "distances.h"
#include "fidlib/fidlib.h"

/* ---------------------------------------------------------------------------- 
   -- Denoising routines                                                     -- 
   ---------------------------------------------------------------------------- */
EEG* eeg_filter_running_median(EEG *eeg, int win, bool alloc){
#ifdef FIXEEG
  EEG *s;
  int c,i;
  if( alloc ){
     s = eeg_clone( eeg, EEG_CLONE_ALL  );
  } else {
     s = eeg;
  }
  for(c=0; c<s->nbchan; c++){
     for(i=0; i<s->ntrials; i++){
        running_median(s->data[c][i], s->n, win);
     }
  }
  return s;
#endif
}

double* running_median(double *d, int n, int win){
  double *tmp, *org;
  int i, awin, ewin;
  double med;
  dprintf("called with %p, of length %i, win=%i\n", d,n,win);
  org = (double*)malloc(n*sizeof(double));
  org = (double*)memcpy(org, d, n*sizeof(double));
  tmp = (double*)malloc((2*win+1)*sizeof(double));


  for(i=0; i<n; i++){
     awin = i-win;
     ewin = i+win;
     if(awin<0) awin=0;
     if(ewin>n-1) ewin=n;
    
    tmp = (double*)memcpy(tmp, &(org[awin]), sizeof(double)*(ewin-awin));
    gsl_sort(tmp, 1, ewin-awin);
     if( (ewin-awin) % 2==0 ){
        med = (tmp[(ewin-awin)/2]+tmp[(ewin-awin)/2+1])/2.0;
     } else {
        med = tmp[(ewin-awin)/2+1];
     }
    d[i] = med;
  }
  free(tmp);
  free(org);
  return d;
}


double* moving_average(double *s, int n, int win){
  int i, awin, ewin;
  double m;
  double *org;

  org = (double*)malloc(n*sizeof(double));
  org = (double*)memcpy(org, s, n*sizeof(double));

  for(i=0; i<n; i++){
     awin = i-win;
     ewin = i+win;
     if(awin<0) awin=0;
     if(ewin>n-1) ewin=n;
     m = dblp_mean( &(org[awin]), ewin-awin);
     s[i] = m;
  }
  free(org);
  return s;
}

EEG* eeg_filter_weighted_running_median( EEG *eeg, int win, bool alloc){
#ifdef FIXEEG
  EEG *s;
  int c,i;
  if( alloc ){
     s = eeg_clone( eeg, EEG_CLONE_ALL );
  } else {
     s = eeg;
  }
  for(c=0; c<s->nbchan; c++){
     for(i=0; i<s->ntrials; i++){
        weighted_running_median(s->data[c][i], s->n, win );
     }
  }
  return s;
#endif
}


double* weighted_running_median(double *d, int n, int win ){
  double *pptr, *wptr;
  int i, j;
  int awin, ewin;
  double med;
  
  pptr = (double*)malloc((2*win+1)*sizeof(double));
  wptr = (double*)malloc((2*win+1)*sizeof(double));

  for(i=0; i<n; i++){
     awin = i-win;
     ewin = i+win;
     if(awin<0) awin=0;
     if(ewin>n-1) ewin=n;

    memcpy(pptr, &(d[awin]), sizeof(double)*(ewin-awin));
    for(j=0; j<(ewin-awin); j++){
      wptr[j] = fabs(i-awin+j);
    }
/*   dprintf("i=%i, awin=%i, ewin=%i\n", i, awin, ewin); */
    med = weighted_median_from_unsorted(pptr, wptr, ewin-awin);
     
    d[i] = med;
  }
  free(pptr);
  free(wptr);

  return d;
}



EEG* eeg_filter_fidlib( EEG *eeg, const char *spec, bool alloc ){
#ifdef FIXEEG
  int len;
  int c,i,j;
  FidFilter *filt;
  char *filtspec, *orgspec;
  FidRun *run;
  FidFunc *funcp;
  void **fbuf;

  EEG *eeg_out;
  if( alloc ){
     eeg_out = eeg_clone( eeg, EEG_CLONE_ALL );
  } else {
     eeg_out = eeg;
  }

  fbuf = (void**) malloc( eeg->ntrials*sizeof(void*));

  filtspec = (char*) malloc( (strlen(spec)+10)*sizeof(char));
  orgspec = filtspec;
  memset( filtspec, 0, (strlen(spec)+10)*sizeof(char));
  strcpy( filtspec, spec );
  fid_parse( eeg_out->sampling_rate, &filtspec, &filt);
  run= fid_run_new(filt, &funcp);

  len= fid_run_bufsize(run);
  for( i=0; i<eeg_out->ntrials; i++ ){
     fbuf[i]= fid_run_newbuf(run);
  }
  dprintf("Buffer initialized\n");\
  for( c=0; c<eeg_out->nbchan; c++ ){
     for( i=0; i<eeg_out->ntrials; i++ ){
        for( j=0; j<eeg_out->n; j++ ){
          //        dprintf("Sample %i, channel %i\n", i, j );
          eeg_out->data[c][i][j] = funcp( fbuf[i], eeg_out->data[c][i][j] );
        }
        fid_run_zapbuf( fbuf[i] );
     }
  }
  for( i=0; i<eeg_out->ntrials; i++ ){
     fid_run_freebuf(fbuf[i]);
  }
  free( fbuf );
  fid_run_free(run);
  free( filt );
  free( orgspec ) ;

  return eeg_out;
#endif
}