Nsound  0.9.4
FilterLowPassFIR.cc
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------
2 //
3 // $Id: FilterLowPassFIR.cc 874 2014-09-08 02:21:29Z weegreenblobbie $
4 //
5 // Copyright (c) 2006 to Present Nick Hilton
6 //
7 // weegreenblobbie_yahoo_com (replace '_' with '@' and '.')
8 //
9 //-----------------------------------------------------------------------------
10 
11 //-----------------------------------------------------------------------------
12 //
13 // This program is free software; you can redistribute it and/or modify
14 // it under the terms of the GNU General Public License as published by
15 // the Free Software Foundation; either version 2 of the License, or
16 // (at your option) any later version.
17 //
18 // This program is distributed in the hope that it will be useful,
19 // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 // GNU Library General Public License for more details.
22 //
23 // You should have received a copy of the GNU General Public License
24 // along with this program; if not, write to the Free Software
25 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
26 //
27 //-----------------------------------------------------------------------------
28 
29 #include <Nsound/AudioStream.h>
30 #include <Nsound/Buffer.h>
32 #include <Nsound/Plotter.h>
33 
34 #include <cmath>
35 #include <cstring>
36 #include <cstdio>
37 #include <iostream>
38 
39 using namespace Nsound;
40 
41 using std::cerr;
42 using std::endl;
43 
44 #define CERR_HEADER __FILE__ << ":" << __LINE__ << ": "
45 
46 //-----------------------------------------------------------------------------
49  const float64 & sample_rate,
50  uint32 kernel_size,
51  const float64 & cutoff_frequency_Hz)
52  :
53  Filter(sample_rate),
54  b_(NULL),
55  window_(NULL),
56  x_history_(NULL),
57  x_ptr_(NULL),
58  x_end_ptr_(NULL),
59  frequency_1_Hz_(cutoff_frequency_Hz),
60  lp_cache_()
61 {
62  kernel_size_ = kernel_size;
63 
64  // Keep kernel_size_ odd.
65  if(kernel_size_ % 2 == 0)
66  {
67  ++kernel_size_;
68  }
69 
70  x_history_ = new float64[kernel_size_ + 1];
73 
74  // Create the Blackman window.
75  float64 ks = static_cast<float64>(kernel_size_);
77  for(uint32 i = 0; i < kernel_size_; ++i)
78  {
79  float64 di = static_cast<float64>(i);
80 
81  window_[i] = 0.420
82  - 0.500 * std::cos((di + 0.5) * 2.0 * M_PI / ks)
83  + 0.080 * std::cos((di + 0.5) * 4.0 * M_PI / ks);
84 
85  }
86 
88 }
89 
90 //-----------------------------------------------------------------------------
91 // This is private and therefor disabled
94  :
95  Filter(copy.sample_rate_),
96  b_(NULL),
97  window_(NULL),
98  x_history_(NULL),
99  x_ptr_(NULL),
100  x_end_ptr_(NULL),
101  frequency_1_Hz_(0.0),
102  lp_cache_()
103 {}
104 
105 //-----------------------------------------------------------------------------
108 {
109  delete [] window_;
110  delete [] x_history_;
111 
112  KernelCache::iterator itor = lp_cache_.begin();
113  KernelCache::iterator end = lp_cache_.end();
114 
115  while(itor != end)
116  {
117  delete [] itor->b_;
118  ++itor;
119  }
120 }
121 
124 filter(const AudioStream & x)
125 {
126  return Filter::filter(x);
127 }
128 
131 filter(const AudioStream & x, const float64 & f)
132 {
133  return Filter::filter(x, f);
134 }
135 
138 filter(const AudioStream & x, const Buffer & frequencies)
139 {
140  return Filter::filter(x, frequencies);
141 }
142 
143 Buffer
145 filter(const Buffer & x)
146 {
147  return Filter::filter(x);
148 }
149 
150 Buffer
152 filter(const Buffer & x, const float64 & f)
153 {
154  // Instead of calling Filter::filter(x,f), we implement this function here
155  // to avoid calling makeKernel(f) for each sample.
156 
157  reset();
158  makeKernel(f);
159 
160  Buffer::const_iterator itor = x.begin();
161  Buffer::const_iterator end = x.end();
162 
163  Buffer y(x.getLength());
164  while(itor != end)
165  {
166  y << FilterLowPassFIR::filter(*itor);
167  ++itor;
168  }
169 
170  return y;
171 }
172 
173 Buffer
175 filter(const Buffer & x, const Buffer & frequencies)
176 {
177  return Filter::filter(x, frequencies);
178 }
179 
180 float64
182 filter(const float64 & x)
183 {
184  // Write x to history
185  *x_ptr_ = x;
186 
187  // Increment history pointer
188  ++x_ptr_;
189 
190  // Bounds check.
191  if(x_ptr_ >= x_end_ptr_)
192  {
193  x_ptr_ = x_history_;
194  }
195 
196  // Perform the convolution.
197 
198  // y[n] = kernel_[0] * x[n]
199  // + kernel_[1] * x[n - 1]
200  // + kernel_[2] * x[n - 2]
201  // ...
202  // + kernel_[N] * x[n - N]
203 
204  float64 y = 0.0;
205  float64 * history = x_ptr_;
206  for(float64 * b = b_; b != b_ + kernel_size_; ++b)
207  {
208  // When we enter this loop, history is pointing at x[n + 1].
209  --history;
210 
211  // Bounds check
212  if(history < x_history_)
213  {
214  history = x_end_ptr_ - 1;
215  }
216 
217  y += *b * *history;
218  }
219 
220  return y;
221 }
222 
223 float64
225 filter(const float64 & x, const float64 & frequency_Hz)
226 {
227  FilterLowPassFIR::makeKernel(frequency_Hz);
228  return filter(x);
229 }
230 
231 Buffer
234 {
236 }
237 
238 void
240 plot(boolean show_fc, boolean show_phase)
241 {
242  char title[256];
243  sprintf(title,
244  "Low Pass FIR Frequency Response\n"
245  "order = %d, fc = %0.1f Hz, sr = %0.1f Hz",
247 
248  Filter::plot(show_phase);
249 
250  Plotter pylab;
251 
252  uint32 n_rows = 1;
253 
254  if(show_phase)
255  {
256  n_rows = 2;
257  }
258 
259  if(show_fc)
260  {
261  pylab.subplot(n_rows, 1, 1);
262 
263  pylab.axvline(frequency_1_Hz_,"color='red'");
264 
265  pylab.title(title);
266  }
267 }
268 
269 void
271 makeKernel(const float64 & cutoff_frequency_Hz)
272 {
273  // Here we create a key into the kernel cache. I'm only storing kernels
274  // with freqs chopped off to the 10th Hz. So filtering with a kernel
275  // designed at 440.1567 Hz will get stored as 440.1. Any other frequency
276  // that starts at 440.1 will not get a kernel calculated and just use the
277  // one in the cache. So the max diff between what's passed in and what's
278  // stored in the cache is 0.09999 Hz.
279 
281  Kernel new_kernel(
282  static_cast<uint32>(cutoff_frequency_Hz));
283 
284  // See if the kernel is in the cache.
285  KernelCache::const_iterator itor = lp_cache_.find(new_kernel);
286 
287  if(itor != lp_cache_.end())
288  {
289  // The kernel was found in the cache.
290  b_ = itor->b_;
291  return;
292  }
293 
294  // The filter wasn't in the cache, need to make it.
295  new_kernel.b_ = new float64[kernel_size_];
296 
297  // Makes a Windowed-sinc filter with cutoff frequency specified.
298  //
299  // DSP: A Practical Guide for Engineers and Scientists
300  //
301  // ISBN-13: 978-0-7506-7444-7
302  //
303  // Equation 16-4
304 
305  const float64 ks = static_cast<float64>(kernel_size_);
306 
307  const float64 omega = two_pi_over_sample_rate_ * cutoff_frequency_Hz;
308 
309  // cut off of zero Hz!
310  if(cutoff_frequency_Hz < 0.10)
311  {
312  // Create no pass filter.
313  memset(new_kernel.b_, 0, kernel_size_ * sizeof(float64));
314  }
315  else
316  {
317  float64 sum = 0.0;
318 
319  float64 ks_2 = ks / 2.0;
320 
321  float64 * b = new_kernel.b_;
322  float64 * w = window_;
323  for(float64 i = 0.0; i < ks; i += 1.0)
324  {
325  float64 x = i - ks_2 + 1e-16;
326 
327  *b = std::sin(omega * x) / x * *w;
328 
329  sum += *b;
330 
331  ++b;
332  ++w;
333  }
334 
335  // Normalize kernel.
336  for(float64 * b = new_kernel.b_; b < new_kernel.b_ + kernel_size_; ++b)
337  {
338  *b /= sum;
339  }
340  }
341 
342  // Point to the new kernel.
343  b_ = new_kernel.b_;
344 
345  // Add the new kernel to the cache.
346  lp_cache_.insert(new_kernel);
347 
348 //~ for(int i = 0; i < kernel_size_; ++i)
349 //~ {
350 //~ cerr << "b_[" << i + 1 << "] = " << b_[i] << endl;
351 //~ }
352 }
353 
354 void
357 {
358  memset(x_history_, 0, sizeof(float64) * (kernel_size_ + 1));
359  x_ptr_ = x_history_;
360 
362 }
363 
364 void
366 setCutoff(const float64 & fc)
367 {
368  frequency_1_Hz_ = fc;
369  reset();
370 }
371 
372 //-----------------------------------------------------------------------------
374 Kernel::
375 Kernel(const uint32 & frequency)
376  :
377  b_(NULL),
378  frequency_(frequency)
379 {
380 }
381 
382 bool
384 Kernel::
385 operator<(const Kernel & rhs) const
386 {
387  return frequency_ < rhs.frequency_;
388 }
389 
void plot(boolean show_fc=true, boolean show_phase=false)
FilterLowPassFIR(const float64 &sample_rate, uint32 kernel_size, const float64 &cutoff_frequency_Hz)
unsigned int uint32
Definition: Nsound.h:153
void axvline(const float64 &x_pos=0.0, const std::string &kwargs="")
Draws a vertical line at x and spans ymin to ymax (ralitive).
Definition: Plotter.cc:358
#define M_PI
Definition: Nsound.h:121
void title(const std::string &title, const std::string &kwargs="")
Add a title to the plot at the top and centered.
Definition: Plotter.cc:1127
double float64
Definition: Nsound.h:146
bool operator<(const Kernel &rhs) const
uint32 getLength() const
Returns the number of samples in the Buffer.
Definition: Buffer.h:587
Buffer getImpulseResponse(const uint32 n_samples=8192)
Definition: Filter.cc:225
void makeKernel(const float64 &frequency1)
Base class for IIR Filters, defines the interface.
Definition: Filter.h:49
Axes subplot(const uint32 n_rows, const uint32 n_cols, const uint32 n, const std::string &kwargs="", Axes *sharex=NULL, Axes *sharey=NULL)
Creates a figure in a subplot, subplot(A, B, C, **kwargs)
Definition: Plotter.cc:1031
iterator end()
Retruns the itreator at the end of the Buffer.
Definition: Buffer.h:348
FloatVector::const_iterator const_iterator
Definition: Buffer.h:70
uint32 kernel_size_
Definition: Filter.h:116
iterator begin()
Retruns the itreator at the start of the Buffer.
Definition: Buffer.h:285
void plot(boolean show_phase=false)
Definition: Filter.cc:262
AudioStream filter(const AudioStream &x)
Definition: Filter.cc:53
A Buffer for storing audio samples.
Definition: Buffer.h:60
Kernel(const uint32 &frequency)
void setCutoff(const float64 &fc)
Sets the cut off frequency (Hz).
A class to store calculated kernels.
float64 two_pi_over_sample_rate_
Definition: Filter.h:114
AudioStream filter(const AudioStream &x)
float64 sample_rate_
Definition: Filter.h:113