IT++ Logo Newcom Logo

freq_filt.h

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

Generated on Sat Aug 25 23:40:24 2007 for IT++ by Doxygen 1.5.2