IT++ Logo
freq_filt.h
Go to the documentation of this file.
1 
29 #ifndef FREQ_FILT_H
30 #define FREQ_FILT_H
31 
32 #include <itpp/base/vec.h>
33 #include <itpp/base/converters.h>
35 #include <itpp/base/matfunc.h>
36 #include <itpp/base/specmat.h>
37 #include <itpp/base/math/min_max.h>
38 #include <itpp/signal/transforms.h>
40 #include <itpp/itexports.h>
41 
42 
43 namespace itpp
44 {
45 
111 template<class Num_T>
113 {
114 public:
122 
134  Freq_Filt(const Vec<Num_T> &b, const int xlength) {init(b, xlength);}
135 
145  Vec<Num_T> filter(const Vec<Num_T> &x, const int strm = 0);
146 
148  int get_fft_size() { return fftsize; }
149 
151  int get_blk_size() { return blksize; }
152 
155 
156 private:
157  int fftsize, blksize;
158  cvec B; // FFT of impulse vector
159  Vec<Num_T> impulse;
160  Vec<Num_T> old_data;
161  cvec zfinal;
162 
163  void init(const Vec<Num_T> &b, const int xlength);
164  vec overlap_add(const vec &x);
165  svec overlap_add(const svec &x);
166  ivec overlap_add(const ivec &x);
167  cvec overlap_add(const cvec &x);
168  void overlap_add(const cvec &x, cvec &y);
169 };
170 
171 
172 // Initialize impulse rsponse, FFT size and block size
173 template <class Num_T>
174 void Freq_Filt<Num_T>::init(const Vec<Num_T> &b, const int xlength)
175 {
176  // Set the impulse response
177  impulse = b;
178 
179  // Get the length of the impulse response
180  int num_elements = impulse.length();
181 
182  // Initizlize old data
183  old_data.set_size(0, false);
184 
185  // Initialize the final state
186  zfinal.set_size(num_elements - 1, false);
187  zfinal.zeros();
188 
189  vec Lvec;
190 
191  /*
192  * Compute the FFT size and the data block size to use.
193  * The FFT size is N = L + M -1, where L corresponds to
194  * the block size (to be determined) and M corresponds
195  * to the impulse length. The method used here is designed
196  * to minimize the (number of blocks) * (number of flops per FFT)
197  * Use the FFTW flop computation of 5*N*log2(N).
198  */
199  vec n = linspace(1, 20, 20);
200  n = pow(2.0, n);
201  ivec fftflops = to_ivec(elem_mult(5.0 * n, log2(n)));
202 
203  // Find where the FFT sizes are > (num_elements-1)
204  //ivec index = find(n,">", static_cast<double>(num_elements-1));
205  ivec index(n.length());
206  int cnt = 0;
207  for (int ii = 0; ii < n.length(); ii++) {
208  if (n(ii) > (num_elements - 1)) {
209  index(cnt) = ii;
210  cnt += 1;
211  }
212  }
213  index.set_size(cnt, true);
214 
215  fftflops = fftflops(index);
216  n = n(index);
217 
218  // Minimize number of blocks * number of flops per FFT
219  Lvec = n - (double)(num_elements - 1);
220  int min_ind = min_index(elem_mult(ceil((double)xlength / Lvec), to_vec(fftflops)));
221  fftsize = static_cast<int>(n(min_ind));
222  blksize = static_cast<int>(Lvec(min_ind));
223 
224  // Finally, compute the FFT of the impulse response
225  B = fft(to_cvec(impulse), fftsize);
226 }
227 
228 // Filter data
229 template <class Num_T>
230 Vec<Num_T> Freq_Filt<Num_T>::filter(const Vec<Num_T> &input, const int strm)
231 {
232  Vec<Num_T> x, tempv;
233  Vec<Num_T> output;
234 
235  /*
236  * If we are not in streaming mode we want to process all of the data.
237  * However, if we are in streaming mode we want to process the data in
238  * blocks that are commensurate with the designed 'blksize' parameter.
239  * So, we do a little book keeping to divide the iinput vector into the
240  * correct size.
241  */
242  if (!strm) {
243  x = input;
244  zfinal.zeros();
245  old_data.set_size(0, false);
246  }
247  else { // we are in streaming mode
248  tempv = concat(old_data, input);
249  if (tempv.length() <= blksize) {
250  x = tempv;
251  old_data.set_size(0, false);
252  }
253  else {
254  int end = tempv.length();
255  int numblks = end / blksize;
256  if ((end % blksize)) {
257  x = tempv(0, blksize * numblks - 1);
258  old_data = tempv(blksize * numblks, end - 1);
259  }
260  else {
261  x = tempv(0, blksize * numblks - 1);
262  old_data.set_size(0, false);
263  }
264  }
265  }
266  output = overlap_add(x);
267 
268  return output;
269 }
270 
271 // Overlap-add routine
272 template<class Num_T>
273 void Freq_Filt<Num_T>::overlap_add(const cvec&x, cvec &y)
274 {
275  int nb = impulse.length();
276  int nx = x.length();
277 
278  y.set_size(nx, false);
279  y.zeros();
280  cvec X, Y;
281  int istart = 0;
282  int L = blksize;
283  while (istart < nx) {
284  int iend = std::min(istart + L - 1, nx - 1);
285 
286  X = fft(x(istart, iend), fftsize);
287  Y = ifft(elem_mult(X, B));
288  Y.set_subvector(0, Y(0, nb - 2) + zfinal);
289  int yend = std::min(nx - 1, istart + fftsize - 1);
290  y.set_subvector(istart, Y(0, yend - istart));
291  zfinal = Y(fftsize - (nb - 1), fftsize - 1);
292  istart += L;
293  }
294 }
295 
296 template<class Num_T>
297 vec Freq_Filt<Num_T>::overlap_add(const vec &x)
298 {
299  cvec y; // Size of y is set later
300  overlap_add(to_cvec(x), y);
301  return real(y);
302 }
303 
304 template<class Num_T>
305 svec Freq_Filt<Num_T>::overlap_add(const svec &x)
306 {
307  cvec y; // Size of y is set later
308  overlap_add(to_cvec(x), y);
309  return to_svec(real(y));
310 }
311 
312 template<class Num_T>
313 ivec Freq_Filt<Num_T>::overlap_add(const ivec &x)
314 {
315  cvec y; // Size of y is set later
316  overlap_add(to_cvec(x), y);
317  return to_ivec(real(y));
318 }
319 
320 template<class Num_T>
321 cvec Freq_Filt<Num_T>::overlap_add(const cvec &x)
322 {
323  cvec y; // Size of y is set later
324  overlap_add(x, y);
325  return y;
326 }
327 
329 
330 // ----------------------------------------------------------------------
331 // Instantiations
332 // ----------------------------------------------------------------------
333 
334 ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<double>;
335 ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<std::complex<double> >;
336 ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<short>;
337 ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<int>;
338 
340 
341 } // namespace itpp
342 
343 #endif // #ifndef FREQ_FILT_H
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:25 for IT++ by Doxygen 1.8.2