IT++ Logo

freq_filt.h

Go to the documentation of this file.
00001 
00030 #ifndef FREQ_FILT_H
00031 #define FREQ_FILT_H
00032 
00033 #include <itpp/base/vec.h>
00034 #include <itpp/base/converters.h>
00035 #include <itpp/base/math/elem_math.h>
00036 #include <itpp/base/matfunc.h>
00037 #include <itpp/base/specmat.h>
00038 #include <itpp/base/math/min_max.h>
00039 
00040 
00041 namespace itpp {
00042 
00108   template<class Num_T>
00109     class Freq_Filt {
00110     public:
00117     Freq_Filt() {}
00118 
00130     Freq_Filt(const Vec<Num_T> &b, const int xlength) {init(b,xlength);}
00131 
00141     Vec<Num_T> filter(const Vec<Num_T> &x, const int strm = 0);
00142 
00144     int get_fft_size() { return fftsize; }
00145 
00147     int get_blk_size() { return blksize; }
00148 
00150     ~Freq_Filt() {}
00151 
00152     private:
00153     int fftsize, blksize;
00154     cvec B;     // FFT of impulse vector
00155     Vec<Num_T> impulse;
00156     Vec<Num_T> old_data;
00157     cvec zfinal;
00158 
00159     void init(const Vec<Num_T> &b, const int xlength);
00160     vec overlap_add(const vec &x);
00161     svec overlap_add(const svec &x);
00162     ivec overlap_add(const ivec &x);
00163     cvec overlap_add(const cvec &x);
00164     void overlap_add(const cvec &x, cvec &y);
00165   };
00166 
00167 
00168   // Initialize impulse rsponse, FFT size and block size
00169   template <class Num_T>
00170     void Freq_Filt<Num_T>::init(const Vec<Num_T> &b, const int xlength)
00171     {
00172       // Set the impulse response
00173       impulse = b;
00174 
00175       // Get the length of the impulse response
00176       int num_elements = impulse.length();
00177 
00178       // Initizlize old data
00179       old_data.set_size(0,false);
00180 
00181       // Initialize the final state
00182       zfinal.set_size(num_elements-1,false);
00183       zfinal.zeros();
00184 
00185       vec Lvec;
00186 
00187       /*
00188        * Compute the FFT size and the data block size to use.
00189        * The FFT size is N = L + M -1, where L corresponds to
00190        * the block size (to be determined) and M corresponds
00191        * to the impulse length. The method used here is designed
00192        * to minimize the (number of blocks) * (number of flops per FFT)
00193        * Use the FFTW flop computation of 5*N*log2(N).
00194        */
00195       vec n = linspace(1,20,20);
00196       n = pow(2.0,n);
00197       ivec fftflops = to_ivec(elem_mult(5.0*n,log2(n)));
00198 
00199       // Find where the FFT sizes are > (num_elements-1)
00200       //ivec index = find(n,">", static_cast<double>(num_elements-1));
00201       ivec index(n.length());
00202       int cnt = 0;
00203       for(int ii=0; ii<n.length(); ii++)
00204         {
00205           if( n(ii) > (num_elements-1) )
00206             {
00207               index(cnt) = ii;
00208               cnt += 1;
00209             }
00210         }
00211       index.set_size(cnt,true);
00212 
00213       fftflops = fftflops(index);
00214       n = n(index);
00215 
00216       // Minimize number of blocks * number of flops per FFT
00217       Lvec = n - (double)(num_elements - 1);
00218       int min_ind = min_index(elem_mult(ceil((double)xlength/Lvec),to_vec(fftflops)));
00219       fftsize = static_cast<int>(n(min_ind));
00220       blksize = static_cast<int>(Lvec(min_ind));
00221 
00222       // Finally, compute the FFT of the impulse response
00223       B = fft(to_cvec(impulse),fftsize);
00224     }
00225 
00226   // Filter data
00227   template <class Num_T>
00228     Vec<Num_T> Freq_Filt<Num_T>::filter(const Vec<Num_T> &input, const int strm)
00229     {
00230       Vec<Num_T> x,tempv;
00231       Vec<Num_T> output;
00232 
00233       /*
00234        * If we are not in streaming mode we want to process all of the data.
00235        * However, if we are in streaming mode we want to process the data in
00236        * blocks that are commensurate with the designed 'blksize' parameter.
00237        * So, we do a little book keeping to divide the iinput vector into the
00238        * correct size.
00239        */
00240       if(!strm)
00241         {
00242           x = input;
00243           zfinal.zeros();
00244           old_data.set_size(0,false);
00245         }
00246       else { // we aare in streaming mode
00247         tempv = concat(old_data,input);
00248         if(tempv.length() <= blksize)
00249           {
00250             x = tempv;
00251             old_data.set_size(0,false);
00252           }
00253         else
00254           {
00255             int end = tempv.length();
00256             int numblks = end / blksize;
00257             if( (end % blksize) )
00258               {
00259                 x = tempv(0,blksize*numblks-1);
00260                 old_data = tempv(blksize*numblks,end-1);
00261               }
00262             else
00263               {
00264                 x = tempv(0,blksize*numblks-1);
00265                 old_data.set_size(0,false);
00266               }
00267           }
00268       }
00269       output = overlap_add(x);
00270 
00271       return output;
00272     }
00273 
00274 } // namespace itpp
00275 
00276 #endif // #ifndef FREQ_FILT_H
SourceForge Logo

Generated on Mon Jan 7 22:28:59 2008 for IT++ by Doxygen 1.5.4