From 2b96ef209016cd93cf40af9ac9bf9a0aed5c7362 Mon Sep 17 00:00:00 2001 From: lgrezesb Date: Fri, 19 Feb 2010 14:01:42 +0000 Subject: [PATCH] fftw --- common/CMakeLists.txt | 3 +- common/clitkSignal.cxx | 190 ++++++++++++++++++++--------------------- common/clitkSignal.h | 16 ++-- 3 files changed, 105 insertions(+), 104 deletions(-) diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index 9f4700e..38eb7e1 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -26,7 +26,8 @@ SET(clitkCommon_SRC clitkOrientation.cxx vvImage.cxx clitkImageToImageGenericFilter.cxx - ) + clitkSignal.cxx +) ADD_LIBRARY(clitkCommon STATIC ${clitkCommon_SRC}) diff --git a/common/clitkSignal.cxx b/common/clitkSignal.cxx index 9f65bbf..25a7db0 100644 --- a/common/clitkSignal.cxx +++ b/common/clitkSignal.cxx @@ -315,117 +315,117 @@ namespace clitk { //--------------------------------------------------------------------- - Signal Signal::HighPassFilter (double sampPeriod, double cutOffFrequency ) - { - //output - Signal temp(m_SamplingPeriod); - temp.resize(size()); - - //the fft - SIGNAL_FFT_TYPE fft; - - //calculate the cut off frequency - unsigned int samp=lrint(cutOffFrequency*static_cast(size())*sampPeriod); - - //forward fft with empty fft - if(fft.size()==0) OneDForwardFourier(*this, fft); - - //remove the frequencies - for(unsigned int i=0;i(0.,0.); - - //backward with remaining frequencies - OneDBackwardFourier(fft,temp); - return temp; - } +// Signal Signal::HighPassFilter (double sampPeriod, double cutOffFrequency ) +// { +// //output +// Signal temp(m_SamplingPeriod); +// temp.resize(size()); +// +// //the fft +// SIGNAL_FFT_TYPE fft; +// +// //calculate the cut off frequency +// unsigned int samp=lrint(cutOffFrequency*static_cast(size())*sampPeriod); +// +// //forward fft with empty fft +// if(fft.size()==0) OneDForwardFourier(*this, fft); +// +// //remove the frequencies +// for(unsigned int i=0;i(0.,0.); +// +// //backward with remaining frequencies +// OneDBackwardFourier(fft,temp); +// return temp; +// } //--------------------------------------------------------------------- //--------------------------------------------------------------------- - Signal Signal::LowPassFilter (double sampPeriod, double cutOffFrequency ) - { - //output - Signal temp(m_SamplingPeriod); - temp.resize(size()); - - //the fft - SIGNAL_FFT_TYPE fft; - - //calculate the cut off frequency - unsigned int samp=lrint(cutOffFrequency*static_cast(size())*sampPeriod); - - //forward fft with empty fft - if(fft.size()==0) OneDForwardFourier(*this, fft); - unsigned int fsize=fft.size(); - - //remove the frequencies - unsigned int limit=min (samp, fsize); - for(unsigned int i=limit;i(0.,0.); - - //backward with remaining frequencies - OneDBackwardFourier(fft,temp); - return temp; - } +// Signal Signal::LowPassFilter (double sampPeriod, double cutOffFrequency ) +// { +// //output +// Signal temp(m_SamplingPeriod); +// temp.resize(size()); +// +// //the fft +// SIGNAL_FFT_TYPE fft; +// +// //calculate the cut off frequency +// unsigned int samp=lrint(cutOffFrequency*static_cast(size())*sampPeriod); +// +// //forward fft with empty fft +// if(fft.size()==0) OneDForwardFourier(*this, fft); +// unsigned int fsize=fft.size(); +// +// //remove the frequencies +// unsigned int limit=min (samp, fsize); +// for(unsigned int i=limit;i(0.,0.); +// +// //backward with remaining frequencies +// OneDBackwardFourier(fft,temp); +// return temp; +// } //--------------------------------------------------------------------- //--------------------------------------------------------------------- - void Signal::OneDForwardFourier(const Signal& input, SIGNAL_FFT_TYPE & fft) - { - //Create output array - fft.resize(input.size()/2+1); - //Temp copy - double *tempCopy=new double[size()]; - copy(begin(), end(), tempCopy); - - //Forward Fourier Transform - fftw_plan p; - p=fftw_plan_dft_r2c_1d(size(),tempCopy,reinterpret_cast(&(fft[0])),FFTW_ESTIMATE); - fftw_execute(p); - fftw_destroy_plan(p); - //delete tempCopy; - return; - } +// void Signal::OneDForwardFourier(const Signal& input, SIGNAL_FFT_TYPE & fft) +// { +// //Create output array +// fft.resize(input.size()/2+1); +// //Temp copy +// double *tempCopy=new double[size()]; +// copy(begin(), end(), tempCopy); +// +// //Forward Fourier Transform +// fftw_plan p; +// p=fftw_plan_dft_r2c_1d(size(),tempCopy,reinterpret_cast(&(fft[0])),FFTW_ESTIMATE); +// fftw_execute(p); +// fftw_destroy_plan(p); +// //delete tempCopy; +// return; +// } //--------------------------------------------------------------------- //--------------------------------------------------------------------- - void Signal::OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output) - { - - //Backward - fftw_plan p; - p=fftw_plan_dft_c2r_1d(output.size(),reinterpret_cast(&(fft[0])),&(output[0]),FFTW_ESTIMATE); - fftw_execute(p); - fftw_destroy_plan(p); - - vector::iterator it=output.begin(); - while(it!=output.end()){ - *it /= (double)output.size(); - it++; - } - return; - } +// void Signal::OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output) +// { +// +// //Backward +// fftw_plan p; +// p=fftw_plan_dft_c2r_1d(output.size(),reinterpret_cast(&(fft[0])),&(output[0]),FFTW_ESTIMATE); +// fftw_execute(p); +// fftw_destroy_plan(p); +// +// vector::iterator it=output.begin(); +// while(it!=output.end()){ +// *it /= (double)output.size(); +// it++; +// } +// return; +// } //--------------------------------------------------------------------- //--------------------------------------------------------------------- - double Signal::MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft) - { - - if(fft.size()==0) OneDForwardFourier(sig,fft); - int posMax=1; - double amplitude, amplitudeMax=abs(fft[1]); - for(unsigned int i=1;iamplitudeMax){ - posMax=i; - amplitudeMax=amplitude; - } - } - return ((double)(posMax)/((double)sig.size()*sig.GetSamplingPeriod())); - } +// double Signal::MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft) +// { +// +// if(fft.size()==0) OneDForwardFourier(sig,fft); +// int posMax=1; +// double amplitude, amplitudeMax=abs(fft[1]); +// for(unsigned int i=1;iamplitudeMax){ +// posMax=i; +// amplitudeMax=amplitude; +// } +// } +// return ((double)(posMax)/((double)sig.size()*sig.GetSamplingPeriod())); +// } //--------------------------------------------------------------------- diff --git a/common/clitkSignal.h b/common/clitkSignal.h index c955936..fc57b78 100644 --- a/common/clitkSignal.h +++ b/common/clitkSignal.h @@ -7,8 +7,8 @@ #include "clitkIO.h" //include external library -#include -#include +//#include +//#include //itk include #include "itkImage.h" @@ -29,7 +29,7 @@ class Signal{ typedef vector< SignalValueType > SignalType; typedef SignalType::iterator iterator; typedef SignalType::const_iterator const_iterator; - typedef vector< complex > SIGNAL_FFT_TYPE; + //typedef vector< complex > SIGNAL_FFT_TYPE; typedef itk::Image ImageType; typedef itk::Vector VectorType; @@ -84,11 +84,11 @@ class Signal{ Signal MovingAverageFilter ( unsigned int length); Signal GaussLikeFilter (); Signal NormalizeMeanStdDev(double newMean=0.5,double newStdDev=0.5); - Signal HighPassFilter (double sampPeriod, double cutOffFrequency ); - Signal LowPassFilter (double sampPeriod, double cutOffFrequency ); - double MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft); - void OneDForwardFourier(const Signal& input,SIGNAL_FFT_TYPE & fft); - void OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output); + //Signal HighPassFilter (double sampPeriod, double cutOffFrequency ); + //Signal LowPassFilter (double sampPeriod, double cutOffFrequency ); + //double MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft); + //void OneDForwardFourier(const Signal& input,SIGNAL_FFT_TYPE & fft); + //void OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output); Signal DetectLocalExtrema(unsigned int width); Signal LimPhase(); Signal MonPhase(); -- 2.45.1