bpmrf/rf_butterworthlowpass.c

Go to the documentation of this file.

#include <bpm/bpm_interface.h>
#include <bpm/bpm_rf.h>
#include <bpm/bpm_alloc.h>
#include <bpm/bpm_nr.h>

int rf_butterworthlowpass( double **RF, int order, double fc ) {


  int i = 0;
  double f;

  if ( ! RF ) {
    bpm_error( "Invalid RF pointer in rf_butterworthlowpass(...)",
               __FILE__, __LINE__ );
    return BPM_FAILURE;
  }

  if ( fc <= 0 ) {
    bpm_error( "Cut off frequency =< 0 in rf_butterworthlowpass(...)",
               __FILE__, __LINE__);
    return BPM_FAILURE;
  }

  if ( order <= 0 ) {
    bpm_error( "Order of filter =< 0 in rf_butterworthlowpass(...)",
               __FILE__, __LINE__);
    return BPM_FAILURE;
  }

  // Allocate space for the FFT
  double **FFT = alloc_complex_wave_double( rf_nsamples );

  // perform FFT
  rf_complexFFT( RF, FFT, NR_FFTFORWARD );
  
  // apply filter
  for( i = 0; i < rf_nsamples; i++ ) {

    // calculate the frequency for the sample...
    // translate the 2nd nyquist band to zero frequency to get
    // mirror image, so automatically filtering twice...
    f = rf_samplefreq * 1.e6 / rf_nsamples * ( i + 0.5 );
    if ( f > rf_samplefreq * 1.e6 / 2. ) f -= rf_samplefreq * 1.e6; 
    FFT[i][Re] *= 1. / (  1. + pow( f / ( fc * 1.e6 ), 2.* order ) );
    FFT[i][Im] *= 1. / (  1. + pow( f / ( fc * 1.e6 ), 2.* order ) );
  }  

  // backward FFT
  rf_complexFFT( FFT, RF, NR_FFTBACKWARD );

  // Apply weighting factor after fft
  for( i = 0; i < rf_nsamples; i++ ) {
    RF[i][Re] *= 2. / (double) rf_nsamples;
    RF[i][Im] *= 2. / (double) rf_nsamples;
  }

  // Clean up
  free_complex_wave_double( FFT, rf_nsamples );

  return BPM_SUCCESS;
}
// end of file

---