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

Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4