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
~Freq_Filt()
Destructor.
Definition: freq_filt.h:154
Various functions on vectors and matrices - header file.
Freq_Filt()
Constructor.
Definition: freq_filt.h:121
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition: mat.h:1582
int get_blk_size()
Return the data block size.
Definition: freq_filt.h:151
ITPP_EXPORT void ifft(const cvec &in, cvec &out)
Inverse Fast Fourier Transform.
ITPP_EXPORT void fft(const cvec &in, cvec &out)
Fast Fourier Transform.
Minimum and maximum functions on vectors and matrices.
int get_fft_size()
Return FFT size.
Definition: freq_filt.h:148
Freq_Filt Frequency domain filtering using the overlap-add techniqueThe Freq_Filt class implements an...
Definition: freq_filt.h:112
Definitions of converters between different vector and matrix types.
T min(const Vec< T > &in)
Minimum value of vector.
Definition: min_max.h:125
Fourier, Hadamard, Walsh-Hadamard, and 2D Hadamard transforms - header file.
svec to_svec(const Vec< T > &v)
Converts a Vec<T> to svec.
Definition: converters.h:65
Freq_Filt(const Vec< Num_T > &b, const int xlength)
Constructor with initialization of the impulse response b.
Definition: freq_filt.h:134
vec log2(const vec &x)
log-2 of the elements
Definition: log_exp.cpp:36
vec pow(const double x, const vec &y)
Calculates x to the power of y (x^y)
Definition: log_exp.h:176
vec to_vec(const Vec< T > &v)
Converts a Vec<T> to vec.
Definition: converters.h:93
itpp namespace
Definition: itmex.h:36
int length() const
The size of the vector.
Definition: vec.h:269
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition: converters.h:79
Definitions of special vectors and matrices.
int min_index(const Vec< T > &in)
Return the postion of the minimum element in the vector.
Definition: min_max.h:234
vec linspace(double from, double to, int points)
linspace (works in the same way as the MATLAB version)
Definition: specmat.cpp:106
Vector Class (Templated)
Definition: factory.h:42
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
Definition: converters.h:107
Vec< Num_T > filter(const Vec< Num_T > &x, const int strm=0)
Filter data in the input vector x.
Definition: freq_filt.h:230
Elementary mathematical functions - header file.
vec real(const cvec &data)
Real part of complex values.
Definition: elem_math.cpp:157
vec ceil(const vec &x)
Round to nearest upper integer.
Definition: converters.h:335
Templated Vector Class Definitions.
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
Definition: array.h:486

Generated on Thu Jun 21 2018 16:06:18 for IT++ by Doxygen 1.8.13