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
Generated on Sun Jul 26 08:36:50 2009 for IT++ by Doxygen 1.5.9