| Array2dScalarReal64* Impala::Core::Array::MakeGaussIIR1d | ( | double | sigma, | |
| int | derivativeOrder, | |||
| int | filterOrder | |||
| ) | [inline] |
Definition at line 249 of file MakeGaussIIR1d.h.
References ComputeStDev(), FillPoleCoef(), FilterCoef(), Impala::Core::Array::complex::im, MAXDERI, MAXRECUR, and Impala::Core::Array::complex::re.
Referenced by Impala::Visualization::Background::Background(), Impala::Samples::Timing::DoOneRecConvSep(), Impala::Application::Src::mainSrc(), RecGauss(), and Impala::Application::WindowBackground::WindowBackground().
00250 { 00251 int size = 2*filterOrder+1; 00252 double* filter = new double[size]; 00253 00254 double q, cc; 00255 double phase, modulus; 00256 int jj; 00257 00258 int nn, mm; 00259 complex bb[MAXRECUR+1]; 00260 complex pp[MAXRECUR+1]; 00261 complex tmp; 00262 double q0, qs, dq; 00263 00264 00265 if (filterOrder > MAXRECUR) 00266 filterOrder = MAXRECUR; 00267 00268 if (derivativeOrder > MAXDERI) { 00269 std::cout << "MakeGaussIIR1d : derivative order too large" << std::endl; 00270 derivativeOrder = MAXDERI; 00271 } 00272 00273 for (jj = 0; jj <= 2*filterOrder; jj++) 00274 filter[jj] = 0; 00275 00276 /* 00277 * Fetch the desired poles, 00278 * depending on the filter order and derivative order 00279 */ 00280 FillPoleCoef( filterOrder, pp, derivativeOrder ); 00281 nn = pp[0].re; /* filter order */ 00282 00283 /* 00284 * Compute the correct value for q (based on the poles in the z-domain) 00285 * that yields a recursive filter with exactly the right 00286 * standard deviation. 00287 */ 00288 q = sigma / 2.0; /* fit was done for sigma = 2 */ 00289 for ( q0 = q, qs = ComputeStDev(nn, pp, q); fabs(sigma - qs)>0.000001; ) { 00290 dq = qs/2.0 - q; 00291 q = q0 - dq; 00292 qs = ComputeStDev(nn, pp, q); 00293 } 00294 00295 /* 00296 * Scale the poles 00297 */ 00298 for (mm = 1; mm <= nn; mm++) 00299 { 00300 modulus = sqrt(pp[mm].re * pp[mm].re + pp[mm].im * pp[mm].im); 00301 phase = atan(pp[mm].im / pp[mm].re); 00302 modulus = exp( log( modulus ) / q ); 00303 phase = phase / q; 00304 pp[mm].re = modulus * cos(phase); 00305 pp[mm].im = modulus * sin(phase); 00306 } 00307 00308 /* 00309 * Compute the actual filter coefficients 00310 */ 00311 for (mm = 0; mm <= nn; mm++) 00312 { 00313 FilterCoef( mm, nn, pp, mm, nn, tmp, bb ); 00314 } 00315 for (mm = 0; mm <= nn; mm++) 00316 { 00317 bb[nn - mm].re *= pow(-1.0, (double) (nn-mm)) / bb[0].re; 00318 bb[nn - mm].im *= pow(-1.0, (double) (nn-mm)) / bb[0].re; 00319 00320 //filter[nn - mm+filterOrder] = HxVec3Double( 00321 // -bb[nn - mm].re, -bb[nn - mm].re, -bb[nn - mm].re); 00322 filter[nn - mm+filterOrder] = -bb[nn - mm].re; 00323 } 00324 00325 /* 00326 * Normalization 00327 */ 00328 cc = 0.0; 00329 for (mm = 0; mm <= nn; mm++) 00330 cc += filter[mm+filterOrder]; 00331 if (cc < 0.0) 00332 cc = -cc; 00333 filter[filterOrder] = cc; 00334 for (mm = 0; mm < filterOrder; mm++) 00335 filter[mm] = filter[2*filterOrder-mm]; 00336 00337 return ArrayCreate<Array2dScalarReal64>(size, 1, 0, 0, filter); 00338 }
Here is the call graph for this function:
