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>
35#include <itpp/base/matfunc.h>
36#include <itpp/base/specmat.h>
40#include <itpp/itexports.h>
41
42
43namespace itpp
44{
45
111template<class Num_T>
113{
114public:
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
156private:
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
173template <class Num_T>
174void 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
229template <class Num_T>
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
272template<class Num_T>
273void 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
296template<class Num_T>
297vec 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
304template<class Num_T>
305svec 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
312template<class Num_T>
313ivec 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
320template<class Num_T>
321cvec 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
334ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<double>;
335ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<std::complex<double> >;
336ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<short>;
337ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<int>;
338
340
341} // namespace itpp
342
343#endif // #ifndef FREQ_FILT_H
General array class.
Definition array.h:105
void set_size(int n, bool copy=false)
Resizing an Array<T>.
Definition array.h:257
friend const Array< T > concat(const Array< T > &a1, const T &e)
Append element e to the end of the Array a.
Definition array.h:486
int length() const
Returns the number of data elements in the array object.
Definition array.h:157
Freq_Filt Frequency domain filtering using the overlap-add technique.
Definition freq_filt.h:113
int get_blk_size()
Return the data block size.
Definition freq_filt.h:151
Freq_Filt(const Vec< Num_T > &b, const int xlength)
Constructor with initialization of the impulse response b.
Definition freq_filt.h:134
int get_fft_size()
Return FFT size.
Definition freq_filt.h:148
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
Freq_Filt()
Constructor.
Definition freq_filt.h:121
~Freq_Filt()
Destructor.
Definition freq_filt.h:154
Definitions of converters between different vector and matrix types.
Elementary mathematical functions - header file.
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.
vec pow(const double x, const vec &y)
Calculates x to the power of y (x^y)
Definition log_exp.h:176
vec log2(const vec &x)
log-2 of the elements
Definition log_exp.cpp:36
int min_index(const Vec< T > &in)
Return the postion of the minimum element in the vector.
Definition min_max.h:234
vec real(const cvec &data)
Real part of complex values.
vec linspace(double from, double to, int points)
linspace (works in the same way as the MATLAB version)
Definition specmat.cpp:106
vec impulse(int size)
Impulse vector.
Definition specmat.cpp:98
Various functions on vectors and matrices - header file.
Minimum and maximum functions on vectors and matrices.
itpp namespace
Definition itmex.h:37
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
Definition converters.h:107
svec to_svec(const Vec< T > &v)
Converts a Vec<T> to svec.
Definition converters.h:65
vec ceil(const vec &x)
Round to nearest upper integer.
Definition converters.h:335
vec to_vec(const Vec< T > &v)
Converts a Vec<T> to vec.
Definition converters.h:93
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition converters.h:79
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
Definitions of special vectors and matrices.
Fourier, Hadamard, Walsh-Hadamard, and 2D Hadamard transforms - header file.
Templated Vector Class Definitions.

Generated on Tue Aug 17 2021 10:59:15 for IT++ by Doxygen 1.9.8