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