icub-basic-demos
Functions | Variables
IIRGausDeriv.cpp File Reference

Coefficients and poles of IIR Gaussian Derivative Filters (order 0, 1 and 2). More...

#include <complex>
#include <iCub/IIRFilt.h>

Go to the source code of this file.

Functions

void calc_poles (int taps, const double scale, const complex< double > oldpoles[], complex< double > newpoles[])
 5 tap second derivative filter with Linf norm approximation More...
 
void calc_coeffs (int taps, const complex< double > poles[], float *coeffs)
 Compute the coefficients of the filter, given its poles. More...
 
void calc_coeffs (int taps, const complex< double > poles[], const double s, float *coeffs)
 Compute the coefficients of the filter for scale s, given the poles at scale 2. More...
 

Variables

const complex< double > d0_N3_Linf []
 

Detailed Description

Coefficients and poles of IIR Gaussian Derivative Filters (order 0, 1 and 2).

See also
"Recursive Gaussian Derivative Filters", L.van.Vliet,I.T.Young and P.W.Verbeek, 1998
Author
Alex Bernardino, ISR-IST
Date
2006-2007
Note
Released under GNU GPL v2.0

Definition in file IIRGausDeriv.cpp.

Function Documentation

◆ calc_coeffs() [1/2]

void calc_coeffs ( int  taps,
const complex< double >  poles[],
const double  s,
float *  coeffs 
)

Compute the coefficients of the filter for scale s, given the poles at scale 2.

Parameters
tapsNumber of taps (3, 4, or 5)
polesPoles of the scale 2 filter.
scaleRequired filter scale (values between 1 and 100 are OK).
coeffsComputed coefficients $(b_0, a_1, a_2, a_3)$:
  • $ b_0 $ is the gain
  • $ (a_1, a_2, a_3)$ are the autoregressive coefficients

Definition at line 251 of file IIRGausDeriv.cpp.

252 {
253 
254  if((taps < 3)||(taps>5))
255  throw "Invalid number of taps in calc_coeffs";
256 
257  if(coeffs == NULL)
258  throw "NULL Pointer argument in calc_coeffs";
259 
260  complex <double> d1_2, d2_2, d3_2, d4_2, d5_2;
261  d1_2 = poles[0];
262  d2_2 = poles[1];
263  d3_2 = poles[2];
264  if(taps > 3)
265  d4_2 = poles[3];
266  else
267  d4_2 = 0;
268  if(taps > 4)
269  d5_2 = poles[4];
270  else
271  d5_2 = 0;
272 
273  double q, std, lambda;
274  double tol = 0.01;
275  complex <double> j(0,1), var;
276  complex <double> d1_s, d2_s, d3_s, d4_s, d5_s;
277  // computing new values for the poles
278  q = s/2;
279  d1_s = exp(log(abs(d1_2))/q)*exp(j*arg(d1_2)/q);
280  d2_s = exp(log(abs(d2_2))/q)*exp(j*arg(d2_2)/q);
281  d3_s = exp(log(abs(d3_2))/q)*exp(j*arg(d3_2)/q);
282  if( abs(d4_2) != 0 )
283  d4_s = exp(log(abs(d4_2))/q)*exp(j*arg(d4_2)/q);
284  else
285  d4_s = 0;
286  if( abs(d5_2) != 0 )
287  d5_s = exp(log(abs(d5_2))/q)*exp(j*arg(d5_2)/q);
288  else
289  d5_s = 0;
290  // computing the variance of the new filter
291  var = d1_s*2.0/(d1_s-1.0)/(d1_s-1.0) + d2_s*2.0/(d2_s-1.0)/(d2_s-1.0) + d3_s*2.0/(d3_s-1.0)/(d3_s-1.0)+d4_s*2.0/(d4_s-1.0)/(d4_s-1.0)+d5_s*2.0/(d5_s-1.0)/(d5_s-1.0);
292  std = sqrt(var.real());
293  while( fabs(s-std) > tol )
294  {
295  lambda = s/std;
296  q = q*lambda;
297  // computing new values for the poles
298  d1_s = exp(log(abs(d1_2))/q)*exp(j*arg(d1_2)/q);
299  d2_s = exp(log(abs(d2_2))/q)*exp(j*arg(d2_2)/q);
300  d3_s = exp(log(abs(d3_2))/q)*exp(j*arg(d3_2)/q);
301  if( abs(d4_2) != 0)
302  d4_s = exp(log(abs(d4_2))/q)*exp(j*arg(d4_2)/q);
303  else
304  d4_s = 0;
305  if( abs(d5_2) != 0)
306  d5_s = exp(log(abs(d5_2))/q)*exp(j*arg(d5_2)/q);
307  else
308  d5_s = 0;
309  // computing the variance of the new filter
310  var = d1_s*2.0/(d1_s-1.0)/(d1_s-1.0) + d2_s*2.0/(d2_s-1.0)/(d2_s-1.0) + d3_s*2.0/(d3_s-1.0)/(d3_s-1.0)+d4_s*2.0/(d4_s-1.0)/(d4_s-1.0)+d5_s*2.0/(d5_s-1.0)/(d5_s-1.0);
311  std = sqrt(var.real());
312  }
313 
314  //computing the filter coeffs
315  if( taps == 3 )
316  {
317  complex<double> b = complex<double>(1.0,0.0)/d1_s/d2_s/d3_s;
318  coeffs[1] = (float)real(-b*(d2_s*d1_s + d3_s*d1_s + d3_s*d2_s));
319  coeffs[2] = (float)real(b*(d1_s + d2_s + d3_s));
320  coeffs[3] = (float)real(-b);
321  coeffs[4] = 0.0f;
322  coeffs[5] = 0.0f;
323  coeffs[0] = 1.0f + coeffs[1] + coeffs[2] + coeffs[3];
324  }
325  else if(taps == 4)
326  {
327  complex<double> b = complex<double>(1.0,0.0)/d1_s/d2_s/d3_s/d4_s;
328  coeffs[1] = (float)real(-b*(d3_s*d2_s*d1_s + d4_s*d2_s*d1_s + d4_s*d3_s*d1_s + d4_s*d3_s*d2_s));
329  coeffs[2] = (float)real(b*(d2_s*d1_s + d3_s*d1_s + d3_s*d2_s + d4_s*d1_s + d4_s*d2_s + d4_s*d3_s));
330  coeffs[3] = (float)real(-b*(d1_s + d2_s + d3_s + d4_s));
331  coeffs[4] = (float)real(b);
332  coeffs[5] = 0.0f;
333  coeffs[0] = 1.0f + coeffs[1] + coeffs[2] + coeffs[3] + coeffs[4];
334  }
335  else if(taps == 5)
336  {
337  complex <double> b = complex<double>(1.0,0.0)/d1_s/d2_s/d3_s/d4_s/d5_s;
338  coeffs[1] = (float)real(-b*(d4_s*d3_s*d2_s*d1_s + d5_s*d3_s*d2_s*d1_s + d5_s*d4_s*d2_s*d1_s + d5_s*d4_s*d3_s*d1_s + d5_s*d4_s*d3_s*d2_s));
339  coeffs[2] = (float)real(b*(d3_s*d2_s*d1_s + d4_s*d2_s*d1_s + d4_s*d3_s*d1_s + d4_s*d3_s*d2_s + d5_s*d2_s*d1_s + d5_s*d3_s*d1_s + d5_s*d3_s*d2_s + d5_s*d4_s*d1_s + d5_s*d4_s*d2_s + d5_s*d4_s*d3_s));
340  coeffs[3] = (float)real(-b*(d2_s*d1_s + d3_s*d1_s + d3_s*d2_s + d4_s*d1_s + d4_s*d2_s + d4_s*d3_s + d5_s*d1_s + d5_s*d2_s + d5_s*d3_s + d5_s*d4_s));
341  coeffs[4] = (float)real(b*(d1_s + d2_s + d3_s + d4_s + d5_s));
342  coeffs[5] = (float)real(-b);
343  coeffs[0] = 1.0f + coeffs[1] + coeffs[2] + coeffs[3] + coeffs[4] + coeffs[5];
344  }
345 }

◆ calc_coeffs() [2/2]

void calc_coeffs ( int  taps,
const complex< double >  poles[],
float *  coeffs 
)

Compute the coefficients of the filter, given its poles.

Parameters
tapsNumber of taps (3, 4, or 5)
polesPoles of the filter.
coeffsComputed coefficients $(b_0, a_1, a_2, a_3)$:
  • $ b_0 $ is the gain
  • $ (a_1, a_2, a_3)$ are the autoregressive coefficients

Definition at line 192 of file IIRGausDeriv.cpp.

193 {
194  if((taps < 3)||(taps>5))
195  throw "Invalid number of taps in calc_coeffs";
196 
197  if(coeffs == NULL)
198  throw "NULL Pointer argument in calc_coeffs";
199 
200  complex <double> d1_s, d2_s, d3_s, d4_s, d5_s;
201  d1_s = poles[0];
202  d2_s = poles[1];
203  d3_s = poles[2];
204  if(taps > 3)
205  d4_s = poles[3];
206  else
207  d4_s = 0;
208  if(taps > 4)
209  d5_s = poles[4];
210  else
211  d5_s = 0;
212 
213  //computing the filter coeffs
214  if( taps == 3 )
215  {
216  complex<double> b = complex<double>(1.0,0.0)/d1_s/d2_s/d3_s;
217  coeffs[1] = (float)real(-b*(d2_s*d1_s + d3_s*d1_s + d3_s*d2_s));
218  coeffs[2] = (float)real(b*(d1_s + d2_s + d3_s));
219  coeffs[3] = (float)real(-b);
220  coeffs[4] = 0.0f;
221  coeffs[5] = 0.0f;
222  coeffs[0] = 1.0f + coeffs[1] + coeffs[2] + coeffs[3];
223  }
224  else if(taps == 4)
225  {
226  complex<double> b = complex<double>(1.0,0.0)/d1_s/d2_s/d3_s/d4_s;
227  coeffs[1] = (float)real(-b*(d3_s*d2_s*d1_s + d4_s*d2_s*d1_s + d4_s*d3_s*d1_s + d4_s*d3_s*d2_s));
228  coeffs[2] = (float)real(b*(d2_s*d1_s + d3_s*d1_s + d3_s*d2_s + d4_s*d1_s + d4_s*d2_s + d4_s*d3_s));
229  coeffs[3] = (float)real(-b*(d1_s + d2_s + d3_s + d4_s));
230  coeffs[4] = (float)real(b);
231  coeffs[5] = 0.0f;
232  coeffs[0] = 1.0f + coeffs[1] + coeffs[2] + coeffs[3] + coeffs[4];
233  }
234  else if(taps == 5)
235  {
236  complex <double> b = complex<double>(1.0,0.0)/d1_s/d2_s/d3_s/d4_s/d5_s;
237  coeffs[1] = (float)real(-b*(d4_s*d3_s*d2_s*d1_s + d5_s*d3_s*d2_s*d1_s + d5_s*d4_s*d2_s*d1_s + d5_s*d4_s*d3_s*d1_s + d5_s*d4_s*d3_s*d2_s));
238  coeffs[2] = (float)real(b*(d3_s*d2_s*d1_s + d4_s*d2_s*d1_s + d4_s*d3_s*d1_s + d4_s*d3_s*d2_s + d5_s*d2_s*d1_s + d5_s*d3_s*d1_s + d5_s*d3_s*d2_s + d5_s*d4_s*d1_s + d5_s*d4_s*d2_s + d5_s*d4_s*d3_s));
239  coeffs[3] = (float)real(-b*(d2_s*d1_s + d3_s*d1_s + d3_s*d2_s + d4_s*d1_s + d4_s*d2_s + d4_s*d3_s + d5_s*d1_s + d5_s*d2_s + d5_s*d3_s + d5_s*d4_s));
240  coeffs[4] = (float)real(b*(d1_s + d2_s + d3_s + d4_s + d5_s));
241  coeffs[5] = (float)real(-b);
242  coeffs[0] = 1.0f + coeffs[1] + coeffs[2] + coeffs[3] + coeffs[4] + coeffs[5];
243  }
244 }

◆ calc_poles()

void calc_poles ( int  taps,
const double  scale,
const complex< double >  oldpoles[],
complex< double >  newpoles[] 
)

5 tap second derivative filter with Linf norm approximation

Compute the poles of the filter for scale s, given the poles at scale 2

Parameters
tapsNumber of taps (3, 4, or 5)
scaleRequired filter scale (values between 1 and 100 are OK).
oldpolesPoles of the scale 2 filter.
newpolesPoles of the computed filter

Definition at line 121 of file IIRGausDeriv.cpp.

122 {
123  if((taps < 3)||(taps>5))
124  throw "Invalid number of taps in calc_poles";
125 
126  if(newpoles == NULL)
127  throw "NULL Pointer argument in calc_poles";
128 
129  complex <double> d1_2, d2_2, d3_2, d4_2, d5_2;
130  d1_2 = oldpoles[0];
131  d2_2 = oldpoles[1];
132  d3_2 = oldpoles[2];
133  if(taps > 3)
134  d4_2 = oldpoles[3];
135  else
136  d4_2 = 0;
137  if(taps > 4)
138  d5_2 = oldpoles[4];
139  else
140  d5_2 = 0;
141 
142  double q, std, lambda;
143  double tol = 0.01;
144  complex <double> j(0,1), var;
145  complex <double> d1_s, d2_s, d3_s, d4_s, d5_s;
146  // computing new values for the poles
147  q = scale/2;
148  d1_s = exp(log(abs(d1_2))/q)*exp(j*arg(d1_2)/q);
149  d2_s = exp(log(abs(d2_2))/q)*exp(j*arg(d2_2)/q);
150  d3_s = exp(log(abs(d3_2))/q)*exp(j*arg(d3_2)/q);
151  if( abs(d4_2) != 0 )
152  d4_s = exp(log(abs(d4_2))/q)*exp(j*arg(d4_2)/q);
153  else
154  d4_s = 0;
155  if( abs(d5_2) != 0 )
156  d5_s = exp(log(abs(d5_2))/q)*exp(j*arg(d5_2)/q);
157  else
158  d5_s = 0;
159  // computing the variance of the new filter
160  var = d1_s*2.0/(d1_s-1.0)/(d1_s-1.0) + d2_s*2.0/(d2_s-1.0)/(d2_s-1.0) + d3_s*2.0/(d3_s-1.0)/(d3_s-1.0)+d4_s*2.0/(d4_s-1.0)/(d4_s-1.0)+d5_s*2.0/(d5_s-1.0)/(d5_s-1.0);
161  std = sqrt(var.real());
162  while( fabs(scale-std) > tol )
163  {
164  lambda = scale/std;
165  q = q*lambda;
166  // computing new values for the poles
167  d1_s = exp(log(abs(d1_2))/q)*exp(j*arg(d1_2)/q);
168  d2_s = exp(log(abs(d2_2))/q)*exp(j*arg(d2_2)/q);
169  d3_s = exp(log(abs(d3_2))/q)*exp(j*arg(d3_2)/q);
170  if( abs(d4_2) != 0)
171  d4_s = exp(log(abs(d4_2))/q)*exp(j*arg(d4_2)/q);
172  else
173  d4_s = 0;
174  if( abs(d5_2) != 0)
175  d5_s = exp(log(abs(d5_2))/q)*exp(j*arg(d5_2)/q);
176  else
177  d5_s = 0;
178  // computing the variance of the new filter
179  var = d1_s*2.0/(d1_s-1.0)/(d1_s-1.0) + d2_s*2.0/(d2_s-1.0)/(d2_s-1.0) + d3_s*2.0/(d3_s-1.0)/(d3_s-1.0)+d4_s*2.0/(d4_s-1.0)/(d4_s-1.0)+d5_s*2.0/(d5_s-1.0)/(d5_s-1.0);
180  std = sqrt(var.real());
181  }
182  newpoles[0] = d1_s;
183  newpoles[1] = d2_s;
184  newpoles[2] = d3_s;
185  newpoles[3] = d4_s;
186  newpoles[4] = d5_s;
187 }

Variable Documentation

◆ d0_N3_Linf

const complex<double> d0_N3_Linf[]
extern
Initial value:
= {
complex<double>(1.40098,1.00236),
complex<double>(1.40098,-1.00236),
complex<double>(1.85132,0),
}
Generated by doxygen 1.9.1