IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
39 
40 namespace itpp
41 {
42 
108 template<class Num_T>
110 {
111 public:
119 
131  Freq_Filt(const Vec<Num_T> &b, const int xlength) {init(b, xlength);}
132 
142  Vec<Num_T> filter(const Vec<Num_T> &x, const int strm = 0);
143 
145  int get_fft_size() { return fftsize; }
146 
148  int get_blk_size() { return blksize; }
149 
152 
153 private:
154  int fftsize, blksize;
155  cvec B; // FFT of impulse vector
156  Vec<Num_T> impulse;
157  Vec<Num_T> old_data;
158  cvec zfinal;
159 
160  void init(const Vec<Num_T> &b, const int xlength);
161  vec overlap_add(const vec &x);
162  svec overlap_add(const svec &x);
163  ivec overlap_add(const ivec &x);
164  cvec overlap_add(const cvec &x);
165  void overlap_add(const cvec &x, cvec &y);
166 };
167 
168 
169 // Initialize impulse rsponse, FFT size and block size
170 template <class Num_T>
171 void Freq_Filt<Num_T>::init(const Vec<Num_T> &b, const int xlength)
172 {
173  // Set the impulse response
174  impulse = b;
175 
176  // Get the length of the impulse response
177  int num_elements = impulse.length();
178 
179  // Initizlize old data
180  old_data.set_size(0, false);
181 
182  // Initialize the final state
183  zfinal.set_size(num_elements - 1, false);
184  zfinal.zeros();
185 
186  vec Lvec;
187 
188  /*
189  * Compute the FFT size and the data block size to use.
190  * The FFT size is N = L + M -1, where L corresponds to
191  * the block size (to be determined) and M corresponds
192  * to the impulse length. The method used here is designed
193  * to minimize the (number of blocks) * (number of flops per FFT)
194  * Use the FFTW flop computation of 5*N*log2(N).
195  */
196  vec n = linspace(1, 20, 20);
197  n = pow(2.0, n);
198  ivec fftflops = to_ivec(elem_mult(5.0 * n, log2(n)));
199 
200  // Find where the FFT sizes are > (num_elements-1)
201  //ivec index = find(n,">", static_cast<double>(num_elements-1));
202  ivec index(n.length());
203  int cnt = 0;
204  for (int ii = 0; ii < n.length(); ii++) {
205  if (n(ii) > (num_elements - 1)) {
206  index(cnt) = ii;
207  cnt += 1;
208  }
209  }
210  index.set_size(cnt, true);
211 
212  fftflops = fftflops(index);
213  n = n(index);
214 
215  // Minimize number of blocks * number of flops per FFT
216  Lvec = n - (double)(num_elements - 1);
217  int min_ind = min_index(elem_mult(ceil((double)xlength / Lvec), to_vec(fftflops)));
218  fftsize = static_cast<int>(n(min_ind));
219  blksize = static_cast<int>(Lvec(min_ind));
220 
221  // Finally, compute the FFT of the impulse response
222  B = fft(to_cvec(impulse), fftsize);
223 }
224 
225 // Filter data
226 template <class Num_T>
227 Vec<Num_T> Freq_Filt<Num_T>::filter(const Vec<Num_T> &input, const int strm)
228 {
229  Vec<Num_T> x, tempv;
230  Vec<Num_T> output;
231 
232  /*
233  * If we are not in streaming mode we want to process all of the data.
234  * However, if we are in streaming mode we want to process the data in
235  * blocks that are commensurate with the designed 'blksize' parameter.
236  * So, we do a little book keeping to divide the iinput vector into the
237  * correct size.
238  */
239  if (!strm) {
240  x = input;
241  zfinal.zeros();
242  old_data.set_size(0, false);
243  }
244  else { // we aare in streaming mode
245  tempv = concat(old_data, input);
246  if (tempv.length() <= blksize) {
247  x = tempv;
248  old_data.set_size(0, false);
249  }
250  else {
251  int end = tempv.length();
252  int numblks = end / blksize;
253  if ((end % blksize)) {
254  x = tempv(0, blksize * numblks - 1);
255  old_data = tempv(blksize * numblks, end - 1);
256  }
257  else {
258  x = tempv(0, blksize * numblks - 1);
259  old_data.set_size(0, false);
260  }
261  }
262  }
263  output = overlap_add(x);
264 
265  return output;
266 }
267 
268 } // namespace itpp
269 
270 #endif // #ifndef FREQ_FILT_H
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:14 for IT++ by Doxygen 1.8.1.2