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
Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4