Nsound  0.9.4
FilterParametricEqualizer.cc
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------
2 //
3 // $Id: FilterParametricEqualizer.cc 874 2014-09-08 02:21:29Z weegreenblobbie $
4 //
5 // Copyright (c) 2007 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 <cstdio>
36 #include <cstring>
37 #include <iostream>
38 
39 using namespace Nsound;
40 
41 static const uint64 N_POLES = 3;
42 
43 //-----------------------------------------------------------------------------
46  const Type & type,
47  const float64 & sample_rate,
48  const float64 & frequency,
49  const float64 & resonance,
50  const float64 & boost)
51  :
52  Filter(sample_rate),
53  type_(type),
54  frequency_(frequency),
55  resonance_(resonance),
56  boost_(boost),
57  a_(new float64 [N_POLES]),
58  b_(new float64 [N_POLES]),
59  x_history_(NULL),
60  x_ptr_(NULL),
61  x_end_ptr_(NULL),
62  y_history_(NULL),
63  y_ptr_(NULL),
64  y_end_ptr_(NULL),
65  kernel_cache_()
66 {
67  x_history_ = new float64 [N_POLES + 1];
70 
71  y_history_ = new float64 [N_POLES + 1];
73  y_end_ptr_ = y_history_ + N_POLES + 1;
74 
75  reset();
76 }
77 
78 //-----------------------------------------------------------------------------
81 {
82  delete [] a_;
83  delete [] b_;
84  delete [] x_history_;
85  delete [] y_history_;
86 }
87 
90 filter(const AudioStream & x)
91 {
92  uint32 n_channels = x.getNChannels();
93 
94  AudioStream y(x.getSampleRate(), n_channels);
95 
96  for(uint32 channel = 0; channel < n_channels; ++channel)
97  {
98  y[channel] = filter(x[channel]);
99  }
100 
101  return y;
102 }
103 
106 filter(const AudioStream & x, const Buffer & frequencies)
107 {
108  uint32 n_channels = x.getNChannels();
109 
110  AudioStream y(x.getSampleRate(), n_channels);
111 
112  for(uint32 channel = 0; channel < n_channels; ++channel)
113  {
114  y[channel] = filter(x[channel], frequencies);
115  }
116 
117  return y;
118 }
119 
123  const AudioStream & x,
124  const Buffer & frequencies,
125  const Buffer & resonance)
126 {
127  uint32 n_channels = x.getNChannels();
128 
129  AudioStream y(x.getSampleRate(), n_channels);
130 
131  for(uint32 channel = 0; channel < n_channels; ++channel)
132  {
133  y[channel] = filter(x[channel], frequencies, resonance);
134  }
135 
136  return y;
137 }
138 
142  const AudioStream & x,
143  const Buffer & frequencies,
144  const Buffer & resonance,
145  const Buffer & boost)
146 {
147  uint32 n_channels = x.getNChannels();
148 
149  AudioStream y(x.getSampleRate(), n_channels);
150 
151  for(uint32 channel = 0; channel < n_channels; ++channel)
152  {
153  y[channel] = filter(x[channel], frequencies, resonance, boost);
154  }
155 
156  return y;
157 }
158 
159 
160 Buffer
162 filter(const Buffer & x)
163 {
164  reset();
165 
166  uint32 n_samples = x.getLength();
167 
168  Buffer y;
169 
170  for(uint32 n = 0; n < n_samples; ++n)
171  {
172  y << filter(x[n]);
173  }
174 
175  return y;
176 }
177 
178 Buffer
180 filter(const Buffer & x, const Buffer & frequencies)
181 {
182  reset();
183 
184  uint32 n_samples = x.getLength();
185  uint32 n_freqs = frequencies.getLength();
186 
187  Buffer y;
188 
189  for(uint32 n = 0; n < n_samples; ++n)
190  {
191  y << filter(x[n], frequencies[n % n_freqs]);
192  }
193 
194  return y;
195 }
196 
197 Buffer
200  const Buffer & x,
201  const Buffer & frequencies,
202  const Buffer & resonance)
203 {
204  reset();
205 
206  uint32 n_samples = x.getLength();
207  uint32 n_freqs = frequencies.getLength();
208  uint32 n_rezzies = resonance.getLength();
209 
210  Buffer y;
211 
212  for(uint32 n = 0; n < n_samples; ++n)
213  {
214  y << filter(x[n], frequencies[n % n_freqs], resonance[n % n_rezzies]);
215  }
216 
217  return y;
218 }
219 
220 Buffer
223  const Buffer & x,
224  const Buffer & frequencies,
225  const Buffer & resonance,
226  const Buffer & boost)
227 {
228  reset();
229 
230  uint32 n_samples = x.getLength();
231  uint32 n_freqs = frequencies.getLength();
232  uint32 n_rezzies = resonance.getLength();
233  uint32 n_boosts = boost.getLength();
234 
235  Buffer y;
236 
237  for(uint32 n = 0; n < n_samples; ++n)
238  {
239  y << filter(x[n],
240  frequencies[n % n_freqs],
241  resonance[n % n_rezzies],
242  boost[n % n_boosts]);
243  }
244 
245  return y;
246 }
247 
248 float64
250 filter(const float64 & x)
251 {
252  // Write x to history
253  *x_ptr_ = x;
254 
255  // Increment history pointer
256  ++x_ptr_;
257 
258  // Bounds check
259  if(x_ptr_ >= x_end_ptr_)
260  {
261  x_ptr_ = x_history_;
262  }
263 
264  float64 y = 0.0;
265  float64 * x_hist = x_ptr_;
266  for(float64 * b = b_; b < b_ + N_POLES; ++b)
267  {
268  // When we enter this loop, x_hist is pointing at x[n + 1]
269  --x_hist;
270 
271  // Bounds check
272  if(x_hist < x_history_)
273  {
274  x_hist = x_end_ptr_ - 1;
275  }
276 
277  y += *b * *x_hist;
278  }
279 
280  float64 * y_hist = y_ptr_;
281  for(float64 * a = a_ + 1; a < a_ + N_POLES; ++a)
282  {
283  // When we enter this loop, y_hist is pointing at y[n + 1]
284  --y_hist;
285 
286  // Bounds check
287  if(y_hist < y_history_)
288  {
289  y_hist = y_end_ptr_ - 1;
290  }
291 
292  y -= *a * *y_hist;
293 
294  }
295 
296  // insert output into history buffer
297  *y_ptr_ = y;
298 
299  // Increment history pointer
300  ++y_ptr_;
301 
302  // Bounds check
303  if(y_ptr_ >= y_end_ptr_)
304  {
305  y_ptr_ = y_history_;
306  }
307 
308  return y;
309 }
310 
311 float64
313 filter(const float64 & x, const float64 & frequency)
314 {
315  makeKernel(frequency, resonance_, boost_);
316 
317  return filter(x);
318 }
319 
320 float64
323  const float64 & x,
324  const float64 & frequency,
325  const float64 & resonance)
326 {
327  makeKernel(frequency, resonance, boost_);
328 
329  return filter(x);
330 }
331 
332 float64
335  const float64 & x,
336  const float64 & frequency,
337  const float64 & resonance,
338  const float64 & boost)
339 {
340  makeKernel(frequency, resonance, boost);
341 
342  return filter(x);
343 }
344 
345 void
348  const float64 & frequency,
349  const float64 & resonance,
350  const float64 & boost)
351 {
352  Kernel new_kernel(N_POLES, N_POLES,
353  static_cast<int32>(std::fabs(frequency)),
354  static_cast<int32>(resonance * 1e3),
355  static_cast<int32>(boost * 1e3));
356 
357  // See if the kernel is in the cache.
358  KernelCache::iterator itor = kernel_cache_.find(new_kernel);
359 
360  // Kernel was found in the cache, set pointers and return.
361  if(itor != kernel_cache_.end())
362  {
363  memcpy(a_, itor->getA(), sizeof(float64) * N_POLES);
364  memcpy(b_, itor->getB(), sizeof(float64) * N_POLES);
365 
366  return;
367  }
368 
369  memset(a_, 0, sizeof(float64) * N_POLES);
370  memset(b_, 0, sizeof(float64) * N_POLES);
371 
372  float64 omega = two_pi_over_sample_rate_ * frequency;
373 
374  switch(type_)
375  {
376  case PEAKING:
377  {
378 
379 //~ K = tan(omega/2)
380 //~
381 //~ b0 = 1 + V*K + K^2
382 //~ b1 = 2*(K^2 - 1)
383 //~ b2 = 1 - V*K + K^2
384 //~
385 //~ a0 = 1 + K/Q + K^2
386 //~ a1 = 2*(K^2 - 1)
387 //~ a2 = 1 - K/Q + K^2
388 
389  float64 K = tan(omega / 2.0);
390  float64 KK = K * K;
391  float64 vkdq = boost * K / resonance;
392 
393  b_[0] = 1.0 + vkdq + KK;
394  b_[1] = 2.0 * (KK - 1.0);
395  b_[2] = 1.0 - vkdq + KK;
396 
397  a_[0] = 1.0 + K / resonance + KK;
398  a_[1] = b_[1];
399  a_[2] = 1.0 - K / resonance + KK;
400 
401  break;
402  }
403 
404  case LOW_SHELF:
405  {
406 
407 //~ K = tan(omega/2)
408 //~
409 //~ b0 = 1 + sqrt(2*V)*K + V*K^2
410 //~ b1 = 2*(V*K^2 - 1)
411 //~ b2 = 1 - sqrt(2*V)*K + V*K^2
412 //~
413 //~ a0 = 1 + K/Q + K^2
414 //~ a1 = 2*(K^2 - 1)
415 //~ a2 = 1 - K/Q + K^2
416 
417  float64 omega = 2.0 * M_PI * frequency / sample_rate_;
418  float64 K = tan(omega / 2.0);
419  float64 KK = K * K;
420 
421  float64 sqvk = sqrt(2.0 * boost) * K;
422  float64 vkk = boost * KK;
423 
424  b_[0] = 1.0 + sqvk + vkk;
425  b_[1] = 2.0 * (vkk - 1.0);
426  b_[2] = 1.0 - sqvk + vkk;
427 
428  a_[0] = 1.0 + K / resonance + KK;
429  a_[1] = 2.0 * (KK - 1.0);
430  a_[2] = 1.0 - K / resonance + KK;
431 
432  break;
433  }
434 
435  case HIGH_SHELF:
436  {
437 
438 //~ K = tan((pi-omega)/2)
439 //~
440 //~ b0 = 1 + sqrt(2*V)*K + V*K^2
441 //~ b1 = -2*(V*K^2 - 1)
442 //~ b1 = 1 - sqrt(2*V)*K + V*K^2
443 //~
444 //~ a0 = 1 + K/Q + K^2
445 //~ a1 = -2*(K^2 - 1)
446 //~ a2 = 1 - K/Q + K^2
447 
448  float64 K = tan((M_PI - omega) * 0.5);
449  float64 KK = K * K;
450  float64 vkk = boost * KK;
451 
452  float64 sq = sqrt(2.0 * boost);
453 
454  b_[0] = 1.0 + sq * K + vkk;
455  b_[1] = -2.0 * (vkk - 1.0);
456  b_[2] = 1.0 - sq * K + vkk;
457 
458  a_[0] = 1.0 + K / resonance + KK;
459  a_[1] = -2.0 * (KK - 1.0);
460  a_[2] = 1.0 - K / resonance + KK;
461 
462  break;
463  }
464  }
465 
466  a_[0] = 1.0 / a_[0];
467 
468  a_[1] *= a_[0];
469  a_[2] *= a_[0];
470 
471  b_[0] *= a_[0];
472  b_[1] *= a_[0];
473  b_[2] *= a_[0];
474 
475  new_kernel.setB(b_);
476  new_kernel.setA(a_);
477 
478  // Add the new kernel to the cache.
479  kernel_cache_.insert(new_kernel);
480 }
481 
482 
483 void
485 plot(boolean show_fc, boolean show_phase)
486 {
487  Filter::plot(show_phase);
488 
489  if(show_fc)
490  {
491  Plotter pylab;
492 
493  pylab.axvline(frequency_,"color='red'");
494 
495  char title[128];
496  sprintf(title,
497  "Parametric Equalizing Frequency Response\n"
498  "order = %.0f, f = %0.1f Hz, sr = %0.1f Hz",
499  static_cast<float64>(N_POLES - 1),
500  frequency_,
501  static_cast<float64>(sample_rate_));
502  pylab.title(title);
503  }
504 }
505 
506 void
509 {
510  memset(x_history_, 0, sizeof(float64) * (N_POLES + 1));
511  memset(y_history_, 0, sizeof(float64) * (N_POLES + 1));
512 
513  x_ptr_ = x_history_;
514  y_ptr_ = y_history_;
515 
517 }
518 
unsigned int uint32
Definition: Nsound.h:153
FilterParametricEqualizer(const Type &type, const float64 &sample_rate, const float64 &frequency, const float64 &resonance=0.707106781187, const float64 &boost_dB=0.0)
boost is in dB
float64 getSampleRate() const
Returns the sample rate of the stream.
Definition: AudioStream.h:217
void plot(boolean show_fc=true, boolean show_phase=false)
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
void setB(const float64 *b)
Copy getBLength() values from the array a into the Kernel. Use setLength() to set the number of value...
Definition: Kernel.cc:233
static const uint64 N_POLES
double float64
Definition: Nsound.h:146
uint32 getLength() const
Returns the number of samples in the Buffer.
Definition: Buffer.h:587
Base class for IIR Filters, defines the interface.
Definition: Filter.h:49
AudioStream filter(const AudioStream &x)
uint32 getNChannels(void) const
Returns the number of audio channels in the stream.
Definition: AudioStream.h:212
unsigned long long uint64
Definition: Nsound.h:154
void plot(boolean show_phase=false)
Definition: Filter.cc:262
void makeKernel(const float64 &frequency, const float64 &resonance, const float64 &boost)
A Buffer for storing audio samples.
Definition: Buffer.h:60
float64 two_pi_over_sample_rate_
Definition: Filter.h:114
float64 sample_rate_
Definition: Filter.h:113
void setA(const float64 *a)
Copy getALength() values from the array a into the Kernel. Use setLength() to set the number of value...
Definition: Kernel.cc:211