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 {
00043 
00109 template<class Num_T>
00110 class Freq_Filt
00111 {
00112 public:
00119   Freq_Filt() {}
00120 
00132   Freq_Filt(const Vec<Num_T> &b, const int xlength) {init(b, xlength);}
00133 
00143   Vec<Num_T> filter(const Vec<Num_T> &x, const int strm = 0);
00144 
00146   int get_fft_size() { return fftsize; }
00147 
00149   int get_blk_size() { return blksize; }
00150 
00152   ~Freq_Filt() {}
00153 
00154 private:
00155   int fftsize, blksize;
00156   cvec B; // FFT of impulse vector
00157   Vec<Num_T> impulse;
00158   Vec<Num_T> old_data;
00159   cvec zfinal;
00160 
00161   void init(const Vec<Num_T> &b, const int xlength);
00162   vec overlap_add(const vec &x);
00163   svec overlap_add(const svec &x);
00164   ivec overlap_add(const ivec &x);
00165   cvec overlap_add(const cvec &x);
00166   void overlap_add(const cvec &x, cvec &y);
00167 };
00168 
00169 
00170 // Initialize impulse rsponse, FFT size and block size
00171 template <class Num_T>
00172 void Freq_Filt<Num_T>::init(const Vec<Num_T> &b, const int xlength)
00173 {
00174   // Set the impulse response
00175   impulse = b;
00176 
00177   // Get the length of the impulse response
00178   int num_elements = impulse.length();
00179 
00180   // Initizlize old data
00181   old_data.set_size(0, false);
00182 
00183   // Initialize the final state
00184   zfinal.set_size(num_elements - 1, false);
00185   zfinal.zeros();
00186 
00187   vec Lvec;
00188 
00189   /*
00190    * Compute the FFT size and the data block size to use.
00191    * The FFT size is N = L + M -1, where L corresponds to
00192    * the block size (to be determined) and M corresponds
00193    * to the impulse length. The method used here is designed
00194    * to minimize the (number of blocks) * (number of flops per FFT)
00195    * Use the FFTW flop computation of 5*N*log2(N).
00196    */
00197   vec n = linspace(1, 20, 20);
00198   n = pow(2.0, n);
00199   ivec fftflops = to_ivec(elem_mult(5.0 * n, log2(n)));
00200 
00201   // Find where the FFT sizes are > (num_elements-1)
00202   //ivec index = find(n,">", static_cast<double>(num_elements-1));
00203   ivec index(n.length());
00204   int cnt = 0;
00205   for (int ii = 0; ii < n.length(); ii++) {
00206     if (n(ii) > (num_elements - 1)) {
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     x = input;
00242     zfinal.zeros();
00243     old_data.set_size(0, false);
00244   }
00245   else { // we aare in streaming mode
00246     tempv = concat(old_data, input);
00247     if (tempv.length() <= blksize) {
00248       x = tempv;
00249       old_data.set_size(0, false);
00250     }
00251     else {
00252       int end = tempv.length();
00253       int numblks = end / blksize;
00254       if ((end % blksize)) {
00255         x = tempv(0, blksize * numblks - 1);
00256         old_data = tempv(blksize * numblks, end - 1);
00257       }
00258       else {
00259         x = tempv(0, blksize * numblks - 1);
00260         old_data.set_size(0, false);
00261       }
00262     }
00263   }
00264   output = overlap_add(x);
00265 
00266   return output;
00267 }
00268 
00269 } // namespace itpp
00270 
00271 #endif // #ifndef FREQ_FILT_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sun Dec 20 07:05:57 2009 for IT++ by Doxygen 1.6.1