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
Generated on Thu Apr 19 14:43:43 2007 for IT++ by Doxygen 1.5.1