| File: | /tmp/pyrefcon/scipy/scipy/signal/sigtoolsmodule.c |
| Warning: | line 1425, column 9 PyObject ownership leak with reference count of 1 |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | /* SIGTOOLS module by Travis Oliphant | |||
| 2 | ||||
| 3 | Copyright 2005 Travis Oliphant | |||
| 4 | Permission to use, copy, modify, and distribute this software without fee | |||
| 5 | is granted under the SciPy License. | |||
| 6 | */ | |||
| 7 | #include <Python.h> | |||
| 8 | #define PY_ARRAY_UNIQUE_SYMBOL_scipy_signal_ARRAY_API _scipy_signal_ARRAY_API | |||
| 9 | #include "numpy/ndarrayobject.h" | |||
| 10 | ||||
| 11 | #include "sigtools.h" | |||
| 12 | #include <setjmp.h> | |||
| 13 | #include <stdlib.h> | |||
| 14 | ||||
| 15 | #define PYERR(message){PyErr_SetString(PyExc_ValueError, message); goto fail;} {PyErr_SetString(PyExc_ValueError, message); goto fail;} | |||
| 16 | ||||
| 17 | ||||
| 18 | jmp_buf MALLOC_FAIL; | |||
| 19 | ||||
| 20 | char *check_malloc(size_t size) | |||
| 21 | { | |||
| 22 | char *the_block = malloc(size); | |||
| 23 | if (the_block == NULL((void*)0)) { | |||
| 24 | printf("\nERROR: unable to allocate %zu bytes!\n", size)__printf_chk (2 - 1, "\nERROR: unable to allocate %zu bytes!\n" , size); | |||
| 25 | longjmp(MALLOC_FAIL,-1); | |||
| 26 | } | |||
| 27 | return the_block; | |||
| 28 | } | |||
| 29 | ||||
| 30 | ||||
| 31 | /************************************************************************ | |||
| 32 | * Start of portable, non-python specific routines. * | |||
| 33 | ************************************************************************/ | |||
| 34 | ||||
| 35 | /* Some core routines are written | |||
| 36 | in a portable way so that they could be used in other applications. The | |||
| 37 | order filtering, however uses python-specific constructs in its guts | |||
| 38 | and is therefore Python dependent. This could be changed in a | |||
| 39 | straightforward way but I haven't done it for lack of time.*/ | |||
| 40 | ||||
| 41 | static int index_out_of_bounds(npy_intp *indices, npy_intp *max_indices, int ndims) { | |||
| 42 | int bad_index = 0, k = 0; | |||
| 43 | ||||
| 44 | while (!bad_index && (k++ < ndims)) { | |||
| 45 | bad_index = ((*(indices) >= *(max_indices++)) || (*(indices) < 0)); | |||
| 46 | indices++; | |||
| 47 | } | |||
| 48 | return bad_index; | |||
| 49 | } | |||
| 50 | ||||
| 51 | /* This maybe could be redone with stride information so it could be | |||
| 52 | * called with non-contiguous arrays: I think offsets is related to | |||
| 53 | * the difference between the strides. I'm not sure about init_offset | |||
| 54 | * just yet. I think it needs to be calculated because of mode_dep | |||
| 55 | * but probably with dim1 being the size of the "original, unsliced" array | |||
| 56 | */ | |||
| 57 | ||||
| 58 | static npy_intp compute_offsets (npy_uintp *offsets, npy_intp *offsets2, npy_intp *dim1, | |||
| 59 | npy_intp *dim2, npy_intp *dim3, npy_intp *mode_dep, | |||
| 60 | int nd) { | |||
| 61 | int k,i; | |||
| 62 | npy_intp init_offset = 0; | |||
| 63 | ||||
| 64 | for (k = 0; k < nd - 1; k++) | |||
| 65 | { | |||
| 66 | init_offset += mode_dep[k]; | |||
| 67 | init_offset *= dim1[k+1]; | |||
| 68 | } | |||
| 69 | init_offset += mode_dep[k] - 2; | |||
| 70 | ||||
| 71 | k = nd; | |||
| 72 | while(k--) { | |||
| 73 | offsets[k] = 0; | |||
| 74 | offsets2[k] = 0; | |||
| 75 | for (i = k + 1; i < nd - 1; i++) { | |||
| 76 | offsets[k] += dim1[i] - dim2[i]; | |||
| 77 | offsets[k] *= dim1[i+1]; | |||
| 78 | ||||
| 79 | offsets2[k] += dim1[i] - dim3[i]; | |||
| 80 | offsets2[k] *= dim1[i+1]; | |||
| 81 | } | |||
| 82 | ||||
| 83 | if (k < nd - 1) { | |||
| 84 | offsets[k] += dim1[i] - dim2[i]; | |||
| 85 | offsets2[k] += dim1[i] - dim3[i]; | |||
| 86 | } | |||
| 87 | offsets[k] += 1; | |||
| 88 | offsets2[k] += 1; | |||
| 89 | } | |||
| 90 | return init_offset; | |||
| 91 | } | |||
| 92 | ||||
| 93 | /* increment by 1 the index into an N-D array, doing the necessary | |||
| 94 | carrying when the index reaches the dimension along that axis */ | |||
| 95 | static int increment(npy_intp *ret_ind, int nd, npy_intp *max_ind) { | |||
| 96 | int k, incr = 1; | |||
| 97 | ||||
| 98 | k = nd - 1; | |||
| 99 | if (++ret_ind[k] >= max_ind[k]) { | |||
| 100 | while (k >= 0 && (ret_ind[k] >= max_ind[k]-1)) { | |||
| 101 | incr++; | |||
| 102 | ret_ind[k--] = 0; | |||
| 103 | } | |||
| 104 | if (k >= 0) ret_ind[k]++; | |||
| 105 | } | |||
| 106 | return incr; | |||
| 107 | } | |||
| 108 | ||||
| 109 | /******************************************************** | |||
| 110 | * | |||
| 111 | * Code taken from remez.c by Erik Kvaleberg which was | |||
| 112 | * converted from an original FORTRAN by | |||
| 113 | * | |||
| 114 | * AUTHORS: JAMES H. MCCLELLAN | |||
| 115 | * | |||
| 116 | * DEPARTMENT OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE | |||
| 117 | * MASSACHUSETTS INSTITUTE OF TECHNOLOGY | |||
| 118 | * CAMBRIDGE, MASS. 02139 | |||
| 119 | * | |||
| 120 | * THOMAS W. PARKS | |||
| 121 | * DEPARTMENT OF ELECTRICAL ENGINEERING | |||
| 122 | * RICE UNIVERSITY | |||
| 123 | * HOUSTON, TEXAS 77001 | |||
| 124 | * | |||
| 125 | * LAWRENCE R. RABINER | |||
| 126 | * BELL LABORATORIES | |||
| 127 | * MURRAY HILL, NEW JERSEY 07974 | |||
| 128 | * | |||
| 129 | * | |||
| 130 | * Adaptation to C by | |||
| 131 | * egil kvaleberg | |||
| 132 | * husebybakken 14a | |||
| 133 | * 0379 oslo, norway | |||
| 134 | * Email: | |||
| 135 | * egil@kvaleberg.no | |||
| 136 | * Web: | |||
| 137 | * http://www.kvaleberg.com/ | |||
| 138 | * | |||
| 139 | * | |||
| 140 | *********************************************************/ | |||
| 141 | ||||
| 142 | ||||
| 143 | #define BANDPASS1 1 | |||
| 144 | #define DIFFERENTIATOR2 2 | |||
| 145 | #define HILBERT3 3 | |||
| 146 | ||||
| 147 | #define GOBACKgoto goto | |||
| 148 | #define DOloop(a,from,to)for ( (a) = (from); (a) <= (to); ++(a)) for ( (a) = (from); (a) <= (to); ++(a)) | |||
| 149 | #define PI3.14159265358979323846 3.14159265358979323846 | |||
| 150 | #define TWOPI(3.14159265358979323846 +3.14159265358979323846) (PI3.14159265358979323846+PI3.14159265358979323846) | |||
| 151 | ||||
| 152 | /* | |||
| 153 | *----------------------------------------------------------------------- | |||
| 154 | * FUNCTION: lagrange_interp (d) | |||
| 155 | * FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION | |||
| 156 | * COEFFICIENTS FOR USE IN THE FUNCTION gee. | |||
| 157 | *----------------------------------------------------------------------- | |||
| 158 | */ | |||
| 159 | static double lagrange_interp(int k, int n, int m, double *x) | |||
| 160 | { | |||
| 161 | int j, l; | |||
| 162 | double q, retval; | |||
| 163 | ||||
| 164 | retval = 1.0; | |||
| 165 | q = x[k]; | |||
| 166 | DOloop(l,1,m)for ( (l) = (1); (l) <= (m); ++(l)) { | |||
| 167 | for (j = l; j <= n; j += m) { | |||
| 168 | if (j != k) | |||
| 169 | retval *= 2.0 * (q - x[j]); | |||
| 170 | } | |||
| 171 | } | |||
| 172 | return 1.0 / retval; | |||
| 173 | } | |||
| 174 | ||||
| 175 | /* | |||
| 176 | *----------------------------------------------------------------------- | |||
| 177 | * FUNCTION: freq_eval (gee) | |||
| 178 | * FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE | |||
| 179 | * LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM | |||
| 180 | *----------------------------------------------------------------------- | |||
| 181 | */ | |||
| 182 | static double freq_eval(int k, int n, double *grid, double *x, double *y, double *ad) | |||
| 183 | { | |||
| 184 | int j; | |||
| 185 | double p,c,d,xf; | |||
| 186 | ||||
| 187 | d = 0.0; | |||
| 188 | p = 0.0; | |||
| 189 | xf = cos(TWOPI(3.14159265358979323846 +3.14159265358979323846) * grid[k]); | |||
| 190 | ||||
| 191 | DOloop(j,1,n)for ( (j) = (1); (j) <= (n); ++(j)) { | |||
| 192 | c = ad[j] / (xf - x[j]); | |||
| 193 | d += c; | |||
| 194 | p += c * y[j]; | |||
| 195 | } | |||
| 196 | ||||
| 197 | return p/d; | |||
| 198 | } | |||
| 199 | ||||
| 200 | ||||
| 201 | /* | |||
| 202 | *----------------------------------------------------------------------- | |||
| 203 | * SUBROUTINE: remez | |||
| 204 | * THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM | |||
| 205 | * FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS | |||
| 206 | * FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE | |||
| 207 | * ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE | |||
| 208 | * DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE | |||
| 209 | * GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE | |||
| 210 | * EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV | |||
| 211 | * ERROR BY DETERMINING THE BSMINEST LOCATION OF THE EXTREMAL | |||
| 212 | * FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES | |||
| 213 | * THE COEFFICIENTS OF THE BEST APPROXIMATION. | |||
| 214 | *----------------------------------------------------------------------- | |||
| 215 | */ | |||
| 216 | static int remez(double *dev, double des[], double grid[], double edge[], | |||
| 217 | double wt[], int ngrid, int nbands, int iext[], double alpha[], | |||
| 218 | int nfcns, int itrmax, double *work, int dimsize, int *niter_out) | |||
| 219 | /* dev, iext, alpha are output types */ | |||
| 220 | /* des, grid, edge, wt, ngrid, nbands, nfcns are input types */ | |||
| 221 | { | |||
| 222 | int k, k1, kkk, kn, knz, klow, kup, nz, nzz, nm1; | |||
| 223 | int cn; | |||
| 224 | int j, jchnge, jet, jm1, jp1; | |||
| 225 | int l, luck=0, nu, nut, nut1=0, niter; | |||
| 226 | ||||
| 227 | double ynz=0.0, comp=0.0, devl, gtemp, fsh, y1=0.0, err, dtemp, delf, dnum, dden; | |||
| 228 | double aa=0.0, bb=0.0, ft, xe, xt; | |||
| 229 | ||||
| 230 | static double *a, *p, *q; | |||
| 231 | static double *ad, *x, *y; | |||
| 232 | ||||
| 233 | a = work; p = a + dimsize+1; q = p + dimsize+1; | |||
| 234 | ad = q + dimsize+1; x = ad + dimsize+1; y = x + dimsize+1; | |||
| 235 | devl = -1.0; | |||
| 236 | nz = nfcns+1; | |||
| 237 | nzz = nfcns+2; | |||
| 238 | niter = 0; | |||
| 239 | ||||
| 240 | do { | |||
| 241 | L100: | |||
| 242 | iext[nzz] = ngrid + 1; | |||
| 243 | ++niter; | |||
| 244 | ||||
| 245 | if (niter > itrmax) break; | |||
| 246 | ||||
| 247 | /* printf("ITERATION %2d: ",niter); */ | |||
| 248 | ||||
| 249 | DOloop(j,1,nz)for ( (j) = (1); (j) <= (nz); ++(j)) { | |||
| 250 | x[j] = cos(grid[iext[j]]*TWOPI(3.14159265358979323846 +3.14159265358979323846)); | |||
| 251 | } | |||
| 252 | jet = (nfcns-1) / 15 + 1; | |||
| 253 | ||||
| 254 | DOloop(j,1,nz)for ( (j) = (1); (j) <= (nz); ++(j)) { | |||
| 255 | ad[j] = lagrange_interp(j,nz,jet,x); | |||
| 256 | } | |||
| 257 | ||||
| 258 | dnum = 0.0; | |||
| 259 | dden = 0.0; | |||
| 260 | k = 1; | |||
| 261 | ||||
| 262 | DOloop(j,1,nz)for ( (j) = (1); (j) <= (nz); ++(j)) { | |||
| 263 | l = iext[j]; | |||
| 264 | dnum += ad[j] * des[l]; | |||
| 265 | dden += (double)k * ad[j] / wt[l]; | |||
| 266 | k = -k; | |||
| 267 | } | |||
| 268 | *dev = dnum / dden; | |||
| 269 | ||||
| 270 | /* printf("DEVIATION = %lg\n",*dev); */ | |||
| 271 | ||||
| 272 | nu = 1; | |||
| 273 | if ( (*dev) > 0.0 ) nu = -1; | |||
| 274 | (*dev) = -(double)nu * (*dev); | |||
| 275 | k = nu; | |||
| 276 | DOloop(j,1,nz)for ( (j) = (1); (j) <= (nz); ++(j)) { | |||
| 277 | l = iext[j]; | |||
| 278 | y[j] = des[l] + (double)k * (*dev) / wt[l]; | |||
| 279 | k = -k; | |||
| 280 | } | |||
| 281 | if ( (*dev) <= devl ) { | |||
| 282 | /* finished */ | |||
| 283 | *niter_out = niter; | |||
| 284 | return -1; | |||
| 285 | } | |||
| 286 | devl = (*dev); | |||
| 287 | jchnge = 0; | |||
| 288 | k1 = iext[1]; | |||
| 289 | knz = iext[nz]; | |||
| 290 | klow = 0; | |||
| 291 | nut = -nu; | |||
| 292 | j = 1; | |||
| 293 | ||||
| 294 | /* | |||
| 295 | * SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST APPROXIMATION | |||
| 296 | */ | |||
| 297 | ||||
| 298 | L200: | |||
| 299 | if (j == nzz) ynz = comp; | |||
| 300 | if (j >= nzz) goto L300; | |||
| 301 | kup = iext[j+1]; | |||
| 302 | l = iext[j]+1; | |||
| 303 | nut = -nut; | |||
| 304 | if (j == 2) y1 = comp; | |||
| 305 | comp = (*dev); | |||
| 306 | if (l >= kup) goto L220; | |||
| 307 | err = (freq_eval(l,nz,grid,x,y,ad)-des[l]) * wt[l]; | |||
| 308 | if (((double)nut*err-comp) <= 0.0) goto L220; | |||
| 309 | comp = (double)nut * err; | |||
| 310 | L210: | |||
| 311 | if (++l >= kup) goto L215; | |||
| 312 | err = (freq_eval(l,nz,grid,x,y,ad)-des[l]) * wt[l]; | |||
| 313 | if (((double)nut*err-comp) <= 0.0) goto L215; | |||
| 314 | comp = (double)nut * err; | |||
| 315 | GOBACKgoto L210; | |||
| 316 | ||||
| 317 | L215: | |||
| 318 | iext[j++] = l - 1; | |||
| 319 | klow = l - 1; | |||
| 320 | ++jchnge; | |||
| 321 | GOBACKgoto L200; | |||
| 322 | ||||
| 323 | L220: | |||
| 324 | --l; | |||
| 325 | L225: | |||
| 326 | if (--l <= klow) goto L250; | |||
| 327 | err = (freq_eval(l,nz,grid,x,y,ad)-des[l]) * wt[l]; | |||
| 328 | if (((double)nut*err-comp) > 0.0) goto L230; | |||
| 329 | if (jchnge <= 0) goto L225; | |||
| 330 | goto L260; | |||
| 331 | ||||
| 332 | L230: | |||
| 333 | comp = (double)nut * err; | |||
| 334 | L235: | |||
| 335 | if (--l <= klow) goto L240; | |||
| 336 | err = (freq_eval(l,nz,grid,x,y,ad)-des[l]) * wt[l]; | |||
| 337 | if (((double)nut*err-comp) <= 0.0) goto L240; | |||
| 338 | comp = (double)nut * err; | |||
| 339 | GOBACKgoto L235; | |||
| 340 | L240: | |||
| 341 | klow = iext[j]; | |||
| 342 | iext[j] = l+1; | |||
| 343 | ++j; | |||
| 344 | ++jchnge; | |||
| 345 | GOBACKgoto L200; | |||
| 346 | ||||
| 347 | L250: | |||
| 348 | l = iext[j]+1; | |||
| 349 | if (jchnge > 0) GOBACKgoto L215; | |||
| 350 | ||||
| 351 | L255: | |||
| 352 | if (++l >= kup) goto L260; | |||
| 353 | err = (freq_eval(l,nz,grid,x,y,ad)-des[l]) * wt[l]; | |||
| 354 | if (((double)nut*err-comp) <= 0.0) GOBACKgoto L255; | |||
| 355 | comp = (double)nut * err; | |||
| 356 | ||||
| 357 | GOBACKgoto L210; | |||
| 358 | L260: | |||
| 359 | klow = iext[j++]; | |||
| 360 | GOBACKgoto L200; | |||
| 361 | ||||
| 362 | L300: | |||
| 363 | if (j > nzz) goto L320; | |||
| 364 | if (k1 > iext[1] ) k1 = iext[1]; | |||
| 365 | if (knz < iext[nz]) knz = iext[nz]; | |||
| 366 | nut1 = nut; | |||
| 367 | nut = -nu; | |||
| 368 | l = 0; | |||
| 369 | kup = k1; | |||
| 370 | comp = ynz*(1.00001); | |||
| 371 | luck = 1; | |||
| 372 | L310: | |||
| 373 | if (++l >= kup) goto L315; | |||
| 374 | err = (freq_eval(l,nz,grid,x,y,ad)-des[l]) * wt[l]; | |||
| 375 | if (((double)nut*err-comp) <= 0.0) GOBACKgoto L310; | |||
| 376 | comp = (double) nut * err; | |||
| 377 | j = nzz; | |||
| 378 | GOBACKgoto L210; | |||
| 379 | ||||
| 380 | L315: | |||
| 381 | luck = 6; | |||
| 382 | goto L325; | |||
| 383 | ||||
| 384 | L320: | |||
| 385 | if (luck > 9) goto L350; | |||
| 386 | if (comp > y1) y1 = comp; | |||
| 387 | k1 = iext[nzz]; | |||
| 388 | L325: | |||
| 389 | l = ngrid+1; | |||
| 390 | klow = knz; | |||
| 391 | nut = -nut1; | |||
| 392 | comp = y1*(1.00001); | |||
| 393 | L330: | |||
| 394 | if (--l <= klow) goto L340; | |||
| 395 | err = (freq_eval(l,nz,grid,x,y,ad)-des[l]) * wt[l]; | |||
| 396 | if (((double)nut*err-comp) <= 0.0) GOBACKgoto L330; | |||
| 397 | j = nzz; | |||
| 398 | comp = (double) nut * err; | |||
| 399 | luck = luck + 10; | |||
| 400 | GOBACKgoto L235; | |||
| 401 | L340: | |||
| 402 | if (luck == 6) goto L370; | |||
| 403 | DOloop(j,1,nfcns)for ( (j) = (1); (j) <= (nfcns); ++(j)) { | |||
| 404 | iext[nzz-j] = iext[nz-j]; | |||
| 405 | } | |||
| 406 | iext[1] = k1; | |||
| 407 | GOBACKgoto L100; | |||
| 408 | L350: | |||
| 409 | kn = iext[nzz]; | |||
| 410 | DOloop(j,1,nfcns)for ( (j) = (1); (j) <= (nfcns); ++(j)) iext[j] = iext[j+1]; | |||
| 411 | iext[nz] = kn; | |||
| 412 | ||||
| 413 | GOBACKgoto L100; | |||
| 414 | L370: | |||
| 415 | ; | |||
| 416 | } while (jchnge > 0); | |||
| 417 | ||||
| 418 | /* | |||
| 419 | * CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION | |||
| 420 | * USING THE INVERSE DISCRETE FOURIER TRANSFORM | |||
| 421 | */ | |||
| 422 | nm1 = nfcns - 1; | |||
| 423 | fsh = 1.0e-06; | |||
| 424 | gtemp = grid[1]; | |||
| 425 | x[nzz] = -2.0; | |||
| 426 | cn = 2*nfcns - 1; | |||
| 427 | delf = 1.0/cn; | |||
| 428 | l = 1; | |||
| 429 | kkk = 0; | |||
| 430 | ||||
| 431 | if (edge[1] == 0.0 && edge[2*nbands] == 0.5) kkk = 1; | |||
| 432 | ||||
| 433 | if (nfcns <= 3) kkk = 1; | |||
| 434 | if (kkk != 1) { | |||
| 435 | dtemp = cos(TWOPI(3.14159265358979323846 +3.14159265358979323846)*grid[1]); | |||
| 436 | dnum = cos(TWOPI(3.14159265358979323846 +3.14159265358979323846)*grid[ngrid]); | |||
| 437 | aa = 2.0/(dtemp-dnum); | |||
| 438 | bb = -(dtemp+dnum)/(dtemp-dnum); | |||
| 439 | } | |||
| 440 | ||||
| 441 | DOloop(j,1,nfcns)for ( (j) = (1); (j) <= (nfcns); ++(j)) { | |||
| 442 | ft = (j - 1) * delf; | |||
| 443 | xt = cos(TWOPI(3.14159265358979323846 +3.14159265358979323846)*ft); | |||
| 444 | if (kkk != 1) { | |||
| 445 | xt = (xt-bb)/aa; | |||
| 446 | #if 0 | |||
| 447 | /*XX* ckeck up !! */ | |||
| 448 | xt1 = sqrt(1.0-xt*xt); | |||
| 449 | ft = atan2(xt1,xt)/TWOPI(3.14159265358979323846 +3.14159265358979323846); | |||
| 450 | #else | |||
| 451 | ft = acos(xt)/TWOPI(3.14159265358979323846 +3.14159265358979323846); | |||
| 452 | #endif | |||
| 453 | } | |||
| 454 | L410: | |||
| 455 | xe = x[l]; | |||
| 456 | if (xt > xe) goto L420; | |||
| 457 | if ((xe-xt) < fsh) goto L415; | |||
| 458 | ++l; | |||
| 459 | GOBACKgoto L410; | |||
| 460 | L415: | |||
| 461 | a[j] = y[l]; | |||
| 462 | goto L425; | |||
| 463 | L420: | |||
| 464 | if ((xt-xe) < fsh) GOBACKgoto L415; | |||
| 465 | grid[1] = ft; | |||
| 466 | a[j] = freq_eval(1,nz,grid,x,y,ad); | |||
| 467 | L425: | |||
| 468 | if (l > 1) l = l-1; | |||
| 469 | } | |||
| 470 | ||||
| 471 | grid[1] = gtemp; | |||
| 472 | dden = TWOPI(3.14159265358979323846 +3.14159265358979323846) / cn; | |||
| 473 | DOloop (j,1,nfcns)for ( (j) = (1); (j) <= (nfcns); ++(j)) { | |||
| 474 | dtemp = 0.0; | |||
| 475 | dnum = (j-1) * dden; | |||
| 476 | if (nm1 >= 1) { | |||
| 477 | DOloop(k,1,nm1)for ( (k) = (1); (k) <= (nm1); ++(k)) { | |||
| 478 | dtemp += a[k+1] * cos(dnum*k); | |||
| 479 | } | |||
| 480 | } | |||
| 481 | alpha[j] = 2.0 * dtemp + a[1]; | |||
| 482 | } | |||
| 483 | ||||
| 484 | DOloop(j,2,nfcns)for ( (j) = (2); (j) <= (nfcns); ++(j)) alpha[j] *= 2.0 / cn; | |||
| 485 | alpha[1] /= cn; | |||
| 486 | ||||
| 487 | if (kkk != 1) { | |||
| 488 | p[1] = 2.0*alpha[nfcns]*bb+alpha[nm1]; | |||
| 489 | p[2] = 2.0*aa*alpha[nfcns]; | |||
| 490 | q[1] = alpha[nfcns-2]-alpha[nfcns]; | |||
| 491 | DOloop(j,2,nm1)for ( (j) = (2); (j) <= (nm1); ++(j)) { | |||
| 492 | if (j >= nm1) { | |||
| 493 | aa *= 0.5; | |||
| 494 | bb *= 0.5; | |||
| 495 | } | |||
| 496 | p[j+1] = 0.0; | |||
| 497 | DOloop(k,1,j)for ( (k) = (1); (k) <= (j); ++(k)) { | |||
| 498 | a[k] = p[k]; | |||
| 499 | p[k] = 2.0 * bb * a[k]; | |||
| 500 | } | |||
| 501 | p[2] += a[1] * 2.0 *aa; | |||
| 502 | jm1 = j - 1; | |||
| 503 | DOloop(k,1,jm1)for ( (k) = (1); (k) <= (jm1); ++(k)) p[k] += q[k] + aa * a[k+1]; | |||
| 504 | jp1 = j + 1; | |||
| 505 | DOloop(k,3,jp1)for ( (k) = (3); (k) <= (jp1); ++(k)) p[k] += aa * a[k-1]; | |||
| 506 | ||||
| 507 | if (j != nm1) { | |||
| 508 | DOloop(k,1,j)for ( (k) = (1); (k) <= (j); ++(k)) q[k] = -a[k]; | |||
| 509 | q[1] += alpha[nfcns - 1 - j]; | |||
| 510 | } | |||
| 511 | } | |||
| 512 | DOloop(j,1,nfcns)for ( (j) = (1); (j) <= (nfcns); ++(j)) alpha[j] = p[j]; | |||
| 513 | } | |||
| 514 | ||||
| 515 | if (nfcns <= 3) { | |||
| 516 | alpha[nfcns+1] = alpha[nfcns+2] = 0.0; | |||
| 517 | } | |||
| 518 | return 0; | |||
| 519 | } | |||
| 520 | ||||
| 521 | ||||
| 522 | /* | |||
| 523 | *----------------------------------------------------------------------- | |||
| 524 | * FUNCTION: eff | |||
| 525 | * FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE | |||
| 526 | * AS A FUNCTION OF FREQUENCY. | |||
| 527 | * AN ARBITRARY FUNCTION OF FREQUENCY CAN BE | |||
| 528 | * APPROXIMATED IF THE USER REPLACES THIS FUNCTION | |||
| 529 | * WITH THE APPROPRIATE CODE TO EVALUATE THE IDEAL | |||
| 530 | * MAGNITUDE. NOTE THAT THE PARAMETER FREQ IS THE | |||
| 531 | * VALUE OF NORMALIZED FREQUENCY NEEDED FOR EVALUATION. | |||
| 532 | *----------------------------------------------------------------------- | |||
| 533 | */ | |||
| 534 | static double eff(double freq, double *fx, int lband, int jtype) | |||
| 535 | { | |||
| 536 | if (jtype != 2) return fx[lband]; | |||
| 537 | else return fx[lband] * freq; | |||
| 538 | } | |||
| 539 | ||||
| 540 | /* | |||
| 541 | *----------------------------------------------------------------------- | |||
| 542 | * FUNCTION: wate | |||
| 543 | * FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION | |||
| 544 | * OF FREQUENCY. SIMILAR TO THE FUNCTION eff, THIS FUNCTION CAN | |||
| 545 | * BE REPLACED BY A USER-WRITTEN ROUTINE TO CALCULATE ANY | |||
| 546 | * DESIRED WEIGHTING FUNCTION. | |||
| 547 | *----------------------------------------------------------------------- | |||
| 548 | */ | |||
| 549 | static double wate(double freq, double *fx, double *wtx, int lband, int jtype) | |||
| 550 | { | |||
| 551 | if (jtype != 2) return wtx[lband]; | |||
| 552 | if (fx[lband] >= 0.0001) return wtx[lband] / freq; | |||
| 553 | return wtx[lband]; | |||
| 554 | } | |||
| 555 | ||||
| 556 | /*********************************************************/ | |||
| 557 | ||||
| 558 | /* This routine accepts basic input information and puts it in | |||
| 559 | * the form expected by remez. | |||
| 560 | ||||
| 561 | * Adpated from main() by Travis Oliphant | |||
| 562 | */ | |||
| 563 | ||||
| 564 | static int pre_remez(double *h2, int numtaps, int numbands, double *bands, | |||
| 565 | double *response, double *weight, int type, int maxiter, | |||
| 566 | int grid_density, int *niter_out) { | |||
| 567 | ||||
| 568 | int jtype, nbands, nfilt, lgrid, nz; | |||
| 569 | int neg, nodd, nm1; | |||
| 570 | int j, k, l, lband, dimsize; | |||
| 571 | double delf, change, fup, temp; | |||
| 572 | double *tempstor, *edge, *h, *fx, *wtx; | |||
| 573 | double *des, *grid, *wt, *alpha, *work; | |||
| 574 | double dev; | |||
| 575 | int ngrid; | |||
| 576 | int *iext; | |||
| 577 | int nfcns, wrksize, total_dsize, total_isize; | |||
| 578 | ||||
| 579 | lgrid = grid_density; | |||
| 580 | dimsize = (int) ceil(numtaps/2.0 + 2); | |||
| 581 | wrksize = grid_density * dimsize; | |||
| 582 | nfilt = numtaps; | |||
| 583 | jtype = type; nbands = numbands; | |||
| 584 | /* Note: code assumes these arrays start at 1 */ | |||
| 585 | edge = bands-1; | |||
| 586 | h = h2 - 1; | |||
| 587 | fx = response - 1; | |||
| 588 | wtx = weight - 1; | |||
| 589 | ||||
| 590 | total_dsize = (dimsize+1)*7 + 3*(wrksize+1); | |||
| 591 | total_isize = (dimsize+1); | |||
| 592 | /* Need space for: (all arrays ignore the first element). | |||
| 593 | ||||
| 594 | des (wrksize+1) | |||
| 595 | grid (wrksize+1) | |||
| 596 | wt (wrksize+1) | |||
| 597 | iext (dimsize+1) (integer) | |||
| 598 | alpha (dimsize+1) | |||
| 599 | work (dimsize+1)*6 | |||
| 600 | ||||
| 601 | */ | |||
| 602 | tempstor = malloc((total_dsize)*sizeof(double)+(total_isize)*sizeof(int)); | |||
| 603 | if (tempstor == NULL((void*)0)) return -2; | |||
| 604 | ||||
| 605 | des = tempstor; grid = des + wrksize+1; | |||
| 606 | wt = grid + wrksize+1; alpha = wt + wrksize+1; | |||
| 607 | work = alpha + dimsize+1; iext = (int *)(work + (dimsize+1)*6); | |||
| 608 | ||||
| 609 | /* Set up problem on dense_grid */ | |||
| 610 | ||||
| 611 | neg = 1; | |||
| 612 | if (jtype == 1) neg = 0; | |||
| 613 | nodd = nfilt % 2; | |||
| 614 | nfcns = nfilt / 2; | |||
| 615 | if (nodd == 1 && neg == 0) nfcns = nfcns + 1; | |||
| 616 | ||||
| 617 | /* | |||
| 618 | * SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID | |||
| 619 | * IS (FILTER LENGTH + 1)*GRID DENSITY/2 | |||
| 620 | */ | |||
| 621 | grid[1] = edge[1]; | |||
| 622 | delf = lgrid * nfcns; | |||
| 623 | delf = 0.5 / delf; | |||
| 624 | if (neg != 0) { | |||
| 625 | if (edge[1] < delf) grid[1] = delf; | |||
| 626 | } | |||
| 627 | j = 1; | |||
| 628 | l = 1; | |||
| 629 | lband = 1; | |||
| 630 | ||||
| 631 | /* | |||
| 632 | * CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT | |||
| 633 | * FUNCTION ON THE GRID | |||
| 634 | */ | |||
| 635 | for (;;) { | |||
| 636 | fup = edge[l + 1]; | |||
| 637 | do { | |||
| 638 | temp = grid[j]; | |||
| 639 | des[j] = eff(temp,fx,lband,jtype); | |||
| 640 | wt[j] = wate(temp,fx,wtx,lband,jtype); | |||
| 641 | if (++j > wrksize) { | |||
| 642 | /* too many points, or too dense grid */ | |||
| 643 | free(tempstor); | |||
| 644 | return -1; | |||
| 645 | } | |||
| 646 | grid[j] = temp + delf; | |||
| 647 | } while (grid[j] <= fup); | |||
| 648 | ||||
| 649 | grid[j-1] = fup; | |||
| 650 | des[j-1] = eff(fup,fx,lband,jtype); | |||
| 651 | wt[j-1] = wate(fup,fx,wtx,lband,jtype); | |||
| 652 | ++lband; | |||
| 653 | l += 2; | |||
| 654 | if (lband > nbands) break; | |||
| 655 | grid[j] = edge[l]; | |||
| 656 | } | |||
| 657 | ||||
| 658 | ngrid = j - 1; | |||
| 659 | if (neg == nodd) { | |||
| 660 | if (grid[ngrid] > (0.5-delf)) --ngrid; | |||
| 661 | } | |||
| 662 | ||||
| 663 | /* | |||
| 664 | * SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT | |||
| 665 | * TO THE ORIGINAL PROBLEM | |||
| 666 | */ | |||
| 667 | if (neg <= 0) { | |||
| 668 | if (nodd != 1) { | |||
| 669 | DOloop(j,1,ngrid)for ( (j) = (1); (j) <= (ngrid); ++(j)) { | |||
| 670 | change = cos(PI3.14159265358979323846*grid[j]); | |||
| 671 | des[j] = des[j] / change; | |||
| 672 | wt[j] = wt[j] * change; | |||
| 673 | } | |||
| 674 | } | |||
| 675 | } else { | |||
| 676 | if (nodd != 1) { | |||
| 677 | DOloop(j,1,ngrid)for ( (j) = (1); (j) <= (ngrid); ++(j)) { | |||
| 678 | change = sin(PI3.14159265358979323846*grid[j]); | |||
| 679 | des[j] = des[j] / change; | |||
| 680 | wt[j] = wt[j] * change; | |||
| 681 | } | |||
| 682 | } else { | |||
| 683 | DOloop(j,1,ngrid)for ( (j) = (1); (j) <= (ngrid); ++(j)) { | |||
| 684 | change = sin(TWOPI(3.14159265358979323846 +3.14159265358979323846) * grid[j]); | |||
| 685 | des[j] = des[j] / change; | |||
| 686 | wt[j] = wt[j] * change; | |||
| 687 | } | |||
| 688 | } | |||
| 689 | } | |||
| 690 | ||||
| 691 | /*XX*/ | |||
| 692 | temp = (double)(ngrid-1) / (double)nfcns; | |||
| 693 | DOloop(j,1,nfcns)for ( (j) = (1); (j) <= (nfcns); ++(j)) { | |||
| 694 | iext[j] = (int)((j-1)*temp) + 1; /* round? !! */ | |||
| 695 | } | |||
| 696 | iext[nfcns+1] = ngrid; | |||
| 697 | nm1 = nfcns - 1; | |||
| 698 | nz = nfcns + 1; | |||
| 699 | ||||
| 700 | if (remez(&dev, des, grid, edge, wt, ngrid, numbands, iext, alpha, nfcns, | |||
| 701 | maxiter, work, dimsize, niter_out) < 0) { | |||
| 702 | free(tempstor); | |||
| 703 | return -1; | |||
| 704 | } | |||
| 705 | ||||
| 706 | /* | |||
| 707 | * CALCULATE THE IMPULSE RESPONSE. | |||
| 708 | */ | |||
| 709 | if (neg <= 0) { | |||
| 710 | ||||
| 711 | if (nodd != 0) { | |||
| 712 | DOloop(j,1,nm1)for ( (j) = (1); (j) <= (nm1); ++(j)) { | |||
| 713 | h[j] = 0.5 * alpha[nz-j]; | |||
| 714 | } | |||
| 715 | h[nfcns] = alpha[1]; | |||
| 716 | } else { | |||
| 717 | h[1] = 0.25 * alpha[nfcns]; | |||
| 718 | DOloop(j,2,nm1)for ( (j) = (2); (j) <= (nm1); ++(j)) { | |||
| 719 | h[j] = 0.25 * (alpha[nz-j] + alpha[nfcns+2-j]); | |||
| 720 | } | |||
| 721 | h[nfcns] = 0.5*alpha[1] + 0.25*alpha[2]; | |||
| 722 | } | |||
| 723 | } else { | |||
| 724 | if (nodd != 0) { | |||
| 725 | h[1] = 0.25 * alpha[nfcns]; | |||
| 726 | h[2] = 0.25 * alpha[nm1]; | |||
| 727 | DOloop(j,3,nm1)for ( (j) = (3); (j) <= (nm1); ++(j)) { | |||
| 728 | h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns+3-j]); | |||
| 729 | } | |||
| 730 | h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[3]; | |||
| 731 | h[nz] = 0.0; | |||
| 732 | } else { | |||
| 733 | h[1] = 0.25 * alpha[nfcns]; | |||
| 734 | DOloop(j,2,nm1)for ( (j) = (2); (j) <= (nm1); ++(j)) { | |||
| 735 | h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns+2-j]); | |||
| 736 | } | |||
| 737 | h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[2]; | |||
| 738 | } | |||
| 739 | } | |||
| 740 | ||||
| 741 | DOloop(j,1,nfcns)for ( (j) = (1); (j) <= (nfcns); ++(j)){ | |||
| 742 | k = nfilt + 1 - j; | |||
| 743 | if (neg == 0) | |||
| 744 | h[k] = h[j]; | |||
| 745 | else | |||
| 746 | h[k] = -h[j]; | |||
| 747 | } | |||
| 748 | if (neg == 1 && nodd == 1) h[nz] = 0.0; | |||
| 749 | ||||
| 750 | free(tempstor); | |||
| 751 | return 0; | |||
| 752 | ||||
| 753 | } | |||
| 754 | ||||
| 755 | /************************************************************** | |||
| 756 | * End of remez routines | |||
| 757 | **************************************************************/ | |||
| 758 | ||||
| 759 | ||||
| 760 | /****************************************************/ | |||
| 761 | /* End of python-independent routines */ | |||
| 762 | /****************************************************/ | |||
| 763 | ||||
| 764 | /************************/ | |||
| 765 | /* N-D Order Filtering. */ | |||
| 766 | ||||
| 767 | ||||
| 768 | static void fill_buffer(char *ip1, PyArrayObject *ap1, PyArrayObject *ap2, | |||
| 769 | char *sort_buffer, int nels2, int check, | |||
| 770 | npy_intp *loop_ind, npy_intp *temp_ind, npy_uintp *offset){ | |||
| 771 | int i, k, incr = 1; | |||
| 772 | int ndims = PyArray_NDIM(ap1); | |||
| 773 | npy_intp *dims2 = PyArray_DIMS(ap2); | |||
| 774 | npy_intp *dims1 = PyArray_DIMS(ap1); | |||
| 775 | npy_intp is1 = PyArray_ITEMSIZE(ap1); | |||
| 776 | npy_intp is2 = PyArray_ITEMSIZE(ap2); | |||
| 777 | char *ip2 = PyArray_DATA(ap2); | |||
| 778 | int elsize = PyArray_ITEMSIZE(ap1); | |||
| 779 | char *ptr; | |||
| 780 | ||||
| 781 | i = nels2; | |||
| 782 | ptr = PyArray_Zero(*(char * (*)(PyArrayObject *)) _scipy_signal_ARRAY_API[47])(ap2); | |||
| 783 | temp_ind[ndims-1]--; | |||
| 784 | while (i--) { | |||
| 785 | /* Adjust index array and move ptr1 to right place */ | |||
| 786 | k = ndims - 1; | |||
| 787 | while(--incr) { | |||
| 788 | temp_ind[k] -= dims2[k] - 1; /* Return to start for these dimensions */ | |||
| 789 | k--; | |||
| 790 | } | |||
| 791 | ip1 += offset[k]*is1; /* Precomputed offset array */ | |||
| 792 | temp_ind[k]++; | |||
| 793 | ||||
| 794 | if (!(check && index_out_of_bounds(temp_ind,dims1,ndims)) && \ | |||
| 795 | memcmp(ip2, ptr, PyArray_ITEMSIZE(ap2))) { | |||
| 796 | memcpy(sort_buffer, ip1, elsize); | |||
| 797 | sort_buffer += elsize; | |||
| 798 | } | |||
| 799 | /* Returns number of N-D indices incremented. */ | |||
| 800 | incr = increment(loop_ind, ndims, dims2); | |||
| 801 | ip2 += is2; | |||
| 802 | ||||
| 803 | } | |||
| 804 | PyDataMem_FREE(*(void (*)(void *)) _scipy_signal_ARRAY_API[289])(ptr); | |||
| 805 | return; | |||
| 806 | } | |||
| 807 | ||||
| 808 | #define COMPARE(fname, type)int fname(type *ip1, type *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } \ | |||
| 809 | int fname(type *ip1, type *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 810 | ||||
| 811 | COMPARE(DOUBLE_compare, double)int DOUBLE_compare(double *ip1, double *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 812 | COMPARE(FLOAT_compare, float)int FLOAT_compare(float *ip1, float *ip2) { return *ip1 < * ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 813 | COMPARE(LONGDOUBLE_compare, npy_longdouble)int LONGDOUBLE_compare(npy_longdouble *ip1, npy_longdouble *ip2 ) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 814 | COMPARE(BYTE_compare, npy_byte)int BYTE_compare(npy_byte *ip1, npy_byte *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 815 | COMPARE(SHORT_compare, short)int SHORT_compare(short *ip1, short *ip2) { return *ip1 < * ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 816 | COMPARE(INT_compare, int)int INT_compare(int *ip1, int *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 817 | COMPARE(LONG_compare, long)int LONG_compare(long *ip1, long *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 818 | COMPARE(LONGLONG_compare, npy_longlong)int LONGLONG_compare(npy_longlong *ip1, npy_longlong *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 819 | COMPARE(UBYTE_compare, npy_ubyte)int UBYTE_compare(npy_ubyte *ip1, npy_ubyte *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 820 | COMPARE(USHORT_compare, npy_ushort)int USHORT_compare(npy_ushort *ip1, npy_ushort *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 821 | COMPARE(UINT_compare, npy_uint)int UINT_compare(npy_uint *ip1, npy_uint *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 822 | COMPARE(ULONG_compare, npy_ulong)int ULONG_compare(npy_ulong *ip1, npy_ulong *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 823 | COMPARE(ULONGLONG_compare, npy_ulonglong)int ULONGLONG_compare(npy_ulonglong *ip1, npy_ulonglong *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; } | |||
| 824 | ||||
| 825 | ||||
| 826 | int OBJECT_compare(PyObject **ip1, PyObject **ip2) { | |||
| 827 | /* PyObject_RichCompareBool returns -1 on error; not handled here */ | |||
| 828 | if(PyObject_RichCompareBool(*ip1, *ip2, Py_LT0) == 1) | |||
| 829 | return -1; | |||
| 830 | else if(PyObject_RichCompareBool(*ip1, *ip2, Py_EQ2) == 1) | |||
| 831 | return 0; | |||
| 832 | else | |||
| 833 | return 1; | |||
| 834 | } | |||
| 835 | ||||
| 836 | typedef int (*CompareFunction)(const void *, const void *); | |||
| 837 | ||||
| 838 | CompareFunction compare_functions[] = \ | |||
| 839 | {NULL((void*)0), (CompareFunction)BYTE_compare,(CompareFunction)UBYTE_compare,\ | |||
| 840 | (CompareFunction)SHORT_compare,(CompareFunction)USHORT_compare, \ | |||
| 841 | (CompareFunction)INT_compare,(CompareFunction)UINT_compare, \ | |||
| 842 | (CompareFunction)LONG_compare,(CompareFunction)ULONG_compare, \ | |||
| 843 | (CompareFunction)LONGLONG_compare,(CompareFunction)ULONGLONG_compare, | |||
| 844 | (CompareFunction)FLOAT_compare,(CompareFunction)DOUBLE_compare, | |||
| 845 | (CompareFunction)LONGDOUBLE_compare, NULL((void*)0), NULL((void*)0), NULL((void*)0), | |||
| 846 | (CompareFunction)OBJECT_compare, NULL((void*)0), NULL((void*)0), NULL((void*)0)}; | |||
| 847 | ||||
| 848 | PyObject *PyArray_OrderFilterND(PyObject *op1, PyObject *op2, int order) { | |||
| 849 | PyArrayObject *ap1=NULL((void*)0), *ap2=NULL((void*)0), *ret=NULL((void*)0); | |||
| 850 | npy_intp *a_ind=NULL((void*)0), *b_ind=NULL((void*)0), *temp_ind=NULL((void*)0), *mode_dep=NULL((void*)0), *check_ind=NULL((void*)0); | |||
| 851 | npy_uintp *offsets=NULL((void*)0); | |||
| 852 | npy_intp *offsets2=NULL((void*)0); | |||
| 853 | npy_uintp offset1; | |||
| 854 | int i, n2, n2_nonzero, k, check, incr = 1; | |||
| 855 | int typenum, bytes_in_array; | |||
| 856 | int is1, os; | |||
| 857 | char *op, *ap1_ptr, *ap2_ptr, *sort_buffer=NULL((void*)0); | |||
| 858 | npy_intp *ret_ind=NULL((void*)0); | |||
| 859 | CompareFunction compare_func=NULL((void*)0); | |||
| 860 | char *zptr=NULL((void*)0); | |||
| 861 | PyArray_CopySwapFunc *copyswap; | |||
| 862 | ||||
| 863 | /* Get Array objects from input */ | |||
| 864 | typenum = PyArray_ObjectType(*(int (*)(PyObject *, int)) _scipy_signal_ARRAY_API[54])(op1, 0); | |||
| 865 | typenum = PyArray_ObjectType(*(int (*)(PyObject *, int)) _scipy_signal_ARRAY_API[54])(op2, typenum); | |||
| 866 | ||||
| 867 | ap1 = (PyArrayObject *)PyArray_ContiguousFromObject(op1, typenum, 0, 0)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(op1, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(typenum), 0, 0, ((0x0001 | (0x0100 | 0x0400))) | 0x0040, ((void*)0)); | |||
| 868 | if (ap1 == NULL((void*)0)) return NULL((void*)0); | |||
| 869 | ap2 = (PyArrayObject *)PyArray_ContiguousFromObject(op2, typenum, 0, 0)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(op2, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(typenum), 0, 0, ((0x0001 | (0x0100 | 0x0400))) | 0x0040, ((void*)0)); | |||
| 870 | if (ap2 == NULL((void*)0)) goto fail; | |||
| 871 | ||||
| 872 | if (PyArray_NDIM(ap1) != PyArray_NDIM(ap2)) { | |||
| 873 | PyErr_SetString(PyExc_ValueError, | |||
| 874 | "All input arrays must have the same number of dimensions."); | |||
| 875 | goto fail; | |||
| 876 | } | |||
| 877 | ||||
| 878 | n2 = PyArray_Size(*(npy_intp (*)(PyObject *)) _scipy_signal_ARRAY_API[59])((PyObject *)ap2); | |||
| 879 | n2_nonzero = 0; | |||
| 880 | ap2_ptr = PyArray_DATA(ap2); | |||
| 881 | /* | |||
| 882 | * Find out the number of non-zero entries in domain (allows for | |||
| 883 | * different shapped rank-filters to be used besides just rectangles) | |||
| 884 | */ | |||
| 885 | zptr = PyArray_Zero(*(char * (*)(PyArrayObject *)) _scipy_signal_ARRAY_API[47])(ap2); | |||
| 886 | if (zptr == NULL((void*)0)) goto fail; | |||
| 887 | for (k=0; k < n2; k++) { | |||
| 888 | n2_nonzero += (memcmp(ap2_ptr,zptr,PyArray_ITEMSIZE(ap2)) != 0); | |||
| 889 | ap2_ptr += PyArray_ITEMSIZE(ap2); | |||
| 890 | } | |||
| 891 | ||||
| 892 | if ((order >= n2_nonzero) || (order < 0)) { | |||
| 893 | PyErr_SetString(PyExc_ValueError, | |||
| 894 | "Order must be non-negative and less than number of nonzero elements in domain."); | |||
| 895 | goto fail; | |||
| 896 | } | |||
| 897 | ||||
| 898 | ret = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ap1),(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int, npy_intp const *, void *, int, int, PyObject *)) _scipy_signal_ARRAY_API [93])(&(*(PyTypeObject *)_scipy_signal_ARRAY_API[2]), PyArray_NDIM (ap1), PyArray_DIMS(ap1), typenum, ((void*)0), ((void*)0), 0, 0, ((void*)0)) | |||
| 899 | PyArray_DIMS(ap1),(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int, npy_intp const *, void *, int, int, PyObject *)) _scipy_signal_ARRAY_API [93])(&(*(PyTypeObject *)_scipy_signal_ARRAY_API[2]), PyArray_NDIM (ap1), PyArray_DIMS(ap1), typenum, ((void*)0), ((void*)0), 0, 0, ((void*)0)) | |||
| 900 | typenum)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int, npy_intp const *, void *, int, int, PyObject *)) _scipy_signal_ARRAY_API [93])(&(*(PyTypeObject *)_scipy_signal_ARRAY_API[2]), PyArray_NDIM (ap1), PyArray_DIMS(ap1), typenum, ((void*)0), ((void*)0), 0, 0, ((void*)0)); | |||
| 901 | if (ret == NULL((void*)0)) goto fail; | |||
| 902 | ||||
| 903 | if (PyArray_TYPE(ap1) < sizeof(compare_functions) / sizeof(compare_functions[0])) { | |||
| 904 | compare_func = compare_functions[PyArray_TYPE(ap1)]; | |||
| 905 | } | |||
| 906 | if (compare_func == NULL((void*)0)) { | |||
| 907 | PyErr_SetString(PyExc_ValueError, | |||
| 908 | "order_filterND not available for this type"); | |||
| 909 | goto fail; | |||
| 910 | } | |||
| 911 | ||||
| 912 | is1 = PyArray_ITEMSIZE(ap1); | |||
| 913 | ||||
| 914 | if (!(sort_buffer = malloc(n2_nonzero*is1))) goto fail; | |||
| 915 | ||||
| 916 | os = PyArray_ITEMSIZE(ret); | |||
| 917 | op = PyArray_DATA(ret); | |||
| 918 | ||||
| 919 | copyswap = PyArray_DESCR(ret)->f->copyswap; | |||
| 920 | ||||
| 921 | bytes_in_array = PyArray_NDIM(ap1)*sizeof(npy_intp); | |||
| 922 | mode_dep = malloc(bytes_in_array); | |||
| 923 | if (mode_dep == NULL((void*)0)) goto fail; | |||
| 924 | for (k = 0; k < PyArray_NDIM(ap1); k++) { | |||
| 925 | mode_dep[k] = -((PyArray_DIMS(ap2)[k]-1) >> 1); | |||
| 926 | } | |||
| 927 | ||||
| 928 | b_ind = (npy_intp *)malloc(bytes_in_array); /* loop variables */ | |||
| 929 | if (b_ind == NULL((void*)0)) goto fail; | |||
| 930 | memset(b_ind,0,bytes_in_array); | |||
| 931 | a_ind = (npy_intp *)malloc(bytes_in_array); | |||
| 932 | ret_ind = (npy_intp *)malloc(bytes_in_array); | |||
| 933 | if (a_ind == NULL((void*)0) || ret_ind == NULL((void*)0)) goto fail; | |||
| 934 | memset(ret_ind,0,bytes_in_array); | |||
| 935 | temp_ind = (npy_intp *)malloc(bytes_in_array); | |||
| 936 | check_ind = (npy_intp*)malloc(bytes_in_array); | |||
| 937 | offsets = (npy_uintp *)malloc(PyArray_NDIM(ap1)*sizeof(npy_uintp)); | |||
| 938 | offsets2 = (npy_intp *)malloc(PyArray_NDIM(ap1)*sizeof(npy_intp)); | |||
| 939 | if (temp_ind == NULL((void*)0) || check_ind == NULL((void*)0) || offsets == NULL((void*)0) || offsets2 == NULL((void*)0)) goto fail; | |||
| 940 | offset1 = compute_offsets(offsets, offsets2, PyArray_DIMS(ap1), | |||
| 941 | PyArray_DIMS(ap2), PyArray_DIMS(ret), | |||
| 942 | mode_dep, PyArray_NDIM(ap1)); | |||
| 943 | /* The filtering proceeds by looping through the output array | |||
| 944 | and for each value filling a buffer from the | |||
| 945 | element-by-element product of the two input arrays. The buffer | |||
| 946 | is then sorted and the order_th element is kept as output. Index | |||
| 947 | counters are used for book-keeping in the area so that we | |||
| 948 | can tell where we are in all of the arrays and be sure that | |||
| 949 | we are not trying to access areas outside the arrays definition. | |||
| 950 | ||||
| 951 | The inner loop is implemented separately but equivalently for each | |||
| 952 | datatype. The outer loop is similar in structure and form to | |||
| 953 | to the inner loop. | |||
| 954 | */ | |||
| 955 | /* Need to keep track of a ptr to place in big (first) input | |||
| 956 | array where we start the multiplication (we pass over it in the | |||
| 957 | inner loop (and not dereferenced) | |||
| 958 | if it is pointing outside dataspace) | |||
| 959 | */ | |||
| 960 | /* Calculate it once and the just move it around appropriately */ | |||
| 961 | PyDataMem_FREE(*(void (*)(void *)) _scipy_signal_ARRAY_API[289])(zptr); | |||
| 962 | zptr = PyArray_Zero(*(char * (*)(PyArrayObject *)) _scipy_signal_ARRAY_API[47])(ap1); | |||
| 963 | if (zptr == NULL((void*)0)) goto fail; | |||
| 964 | ap1_ptr = (char *)PyArray_DATA(ap1) + offset1*is1; | |||
| 965 | for (k=0; k < PyArray_NDIM(ap1); k++) { | |||
| 966 | a_ind[k] = mode_dep[k]; | |||
| 967 | check_ind[k] = PyArray_DIMS(ap1)[k] - PyArray_DIMS(ap2)[k] - mode_dep[k] - 1; | |||
| 968 | } | |||
| 969 | a_ind[PyArray_NDIM(ap1)-1]--; | |||
| 970 | i = PyArray_Size(*(npy_intp (*)(PyObject *)) _scipy_signal_ARRAY_API[59])((PyObject *)ret); | |||
| 971 | while (i--) { | |||
| 972 | /* | |||
| 973 | * Zero out the sort_buffer (has effect of zero-padding | |||
| 974 | * on boundaries). Treat object arrays right. | |||
| 975 | */ | |||
| 976 | ap2_ptr = sort_buffer; | |||
| 977 | for (k=0; k < n2_nonzero; k++) { | |||
| 978 | memcpy(ap2_ptr,zptr,is1); | |||
| 979 | ap2_ptr += is1; | |||
| 980 | } | |||
| 981 | ||||
| 982 | k = PyArray_NDIM(ap1) - 1; | |||
| 983 | while(--incr) { | |||
| 984 | a_ind[k] -= PyArray_DIMS(ret)[k] - 1; /* Return to start */ | |||
| 985 | k--; | |||
| 986 | } | |||
| 987 | ap1_ptr += offsets2[k]*is1; | |||
| 988 | a_ind[k]++; | |||
| 989 | memcpy(temp_ind, a_ind, bytes_in_array); | |||
| 990 | ||||
| 991 | check = 0; k = -1; | |||
| 992 | while(!check && (++k < PyArray_NDIM(ap1))) | |||
| 993 | check = (check || (ret_ind[k] < -mode_dep[k]) || | |||
| 994 | (ret_ind[k] > check_ind[k])); | |||
| 995 | ||||
| 996 | fill_buffer(ap1_ptr,ap1,ap2,sort_buffer,n2,check,b_ind,temp_ind,offsets); | |||
| 997 | qsort(sort_buffer, n2_nonzero, is1, compare_func); | |||
| 998 | ||||
| 999 | /* | |||
| 1000 | * Use copyswap for correct refcounting with object arrays | |||
| 1001 | * (sort_buffer has borrowed references, op owns references). Note | |||
| 1002 | * also that os == PyArray_ITEMSIZE(ret) and we are copying a single | |||
| 1003 | * scalar here. | |||
| 1004 | */ | |||
| 1005 | copyswap(op, sort_buffer + order*is1, 0, NULL((void*)0)); | |||
| 1006 | ||||
| 1007 | /* increment index counter */ | |||
| 1008 | incr = increment(ret_ind,PyArray_NDIM(ret),PyArray_DIMS(ret)); | |||
| 1009 | /* increment to next output index */ | |||
| 1010 | op += os; | |||
| 1011 | ||||
| 1012 | } | |||
| 1013 | free(b_ind); free(a_ind); free(ret_ind); | |||
| 1014 | free(offsets); free(offsets2); free(temp_ind); | |||
| 1015 | free(check_ind); free(mode_dep); | |||
| 1016 | free(sort_buffer); | |||
| 1017 | ||||
| 1018 | PyDataMem_FREE(*(void (*)(void *)) _scipy_signal_ARRAY_API[289])(zptr); | |||
| 1019 | Py_DECREF(ap1)_Py_DECREF(((PyObject*)(ap1))); | |||
| 1020 | Py_DECREF(ap2)_Py_DECREF(((PyObject*)(ap2))); | |||
| 1021 | ||||
| 1022 | return PyArray_Return(*(PyObject * (*)(PyArrayObject *)) _scipy_signal_ARRAY_API[76 ])(ret); | |||
| 1023 | ||||
| 1024 | fail: | |||
| 1025 | if (zptr) PyDataMem_FREE(*(void (*)(void *)) _scipy_signal_ARRAY_API[289])(zptr); | |||
| 1026 | free(sort_buffer); | |||
| 1027 | free(mode_dep); | |||
| 1028 | free(b_ind); | |||
| 1029 | free(a_ind); | |||
| 1030 | free(ret_ind); | |||
| 1031 | free(temp_ind); | |||
| 1032 | free(check_ind); | |||
| 1033 | free(offsets); | |||
| 1034 | free(offsets2); | |||
| 1035 | Py_XDECREF(ap1)_Py_XDECREF(((PyObject*)(ap1))); | |||
| 1036 | Py_XDECREF(ap2)_Py_XDECREF(((PyObject*)(ap2))); | |||
| 1037 | Py_XDECREF(ret)_Py_XDECREF(((PyObject*)(ret))); | |||
| 1038 | return NULL((void*)0); | |||
| 1039 | } | |||
| 1040 | ||||
| 1041 | ||||
| 1042 | /******************************************/ | |||
| 1043 | ||||
| 1044 | static char doc_correlateND[] = "out = _correlateND(a,kernel,mode) \n\n mode = 0 - 'valid', 1 - 'same', \n 2 - 'full' (default)"; | |||
| 1045 | ||||
| 1046 | /*******************************************************************/ | |||
| 1047 | ||||
| 1048 | static char doc_convolve2d[] = "out = _convolve2d(in1, in2, flip, mode, boundary, fillvalue)"; | |||
| 1049 | ||||
| 1050 | extern int pylab_convolve_2d(char*, npy_intp*, char*, npy_intp*, char*, | |||
| 1051 | npy_intp*, npy_intp*, npy_intp*, int, char*); | |||
| 1052 | ||||
| 1053 | static PyObject *sigtools_convolve2d(PyObject *NPY_UNUSED(dummy)(__NPY_UNUSED_TAGGEDdummy) __attribute__ ((__unused__)), PyObject *args) { | |||
| 1054 | ||||
| 1055 | PyObject *in1=NULL((void*)0), *in2=NULL((void*)0), *fill_value=NULL((void*)0); | |||
| 1056 | int mode=2, boundary=0, typenum, flag, flip=1, ret; | |||
| 1057 | npy_intp *aout_dimens=NULL((void*)0); | |||
| 1058 | int i; | |||
| 1059 | PyArrayObject *ain1=NULL((void*)0), *ain2=NULL((void*)0), *aout=NULL((void*)0); | |||
| 1060 | PyArrayObject *afill=NULL((void*)0); | |||
| 1061 | ||||
| 1062 | if (!PyArg_ParseTuple(args, "OO|iiiO", &in1, &in2, &flip, &mode, &boundary, | |||
| 1063 | &fill_value)) { | |||
| 1064 | return NULL((void*)0); | |||
| 1065 | } | |||
| 1066 | ||||
| 1067 | typenum = PyArray_ObjectType(*(int (*)(PyObject *, int)) _scipy_signal_ARRAY_API[54])(in1, 0); | |||
| 1068 | typenum = PyArray_ObjectType(*(int (*)(PyObject *, int)) _scipy_signal_ARRAY_API[54])(in2, typenum); | |||
| 1069 | ain1 = (PyArrayObject *)PyArray_FromObject(in1, typenum, 2, 2)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(in1, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(typenum), 2, 2, (0x0100 | 0x0400) | 0x0040, ((void*)0)); | |||
| 1070 | if (ain1 == NULL((void*)0)) goto fail; | |||
| 1071 | ain2 = (PyArrayObject *)PyArray_FromObject(in2, typenum, 2, 2)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(in2, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(typenum), 2, 2, (0x0100 | 0x0400) | 0x0040, ((void*)0)); | |||
| 1072 | if (ain2 == NULL((void*)0)) goto fail; | |||
| 1073 | ||||
| 1074 | if ((boundary != PAD0) && (boundary != REFLECT4) && (boundary != CIRCULAR8)) | |||
| 1075 | PYERR("Incorrect boundary value."){PyErr_SetString(PyExc_ValueError, "Incorrect boundary value." ); goto fail;}; | |||
| 1076 | ||||
| 1077 | if ((boundary == PAD0) & (fill_value != NULL((void*)0))) { | |||
| 1078 | afill = (PyArrayObject *)PyArray_FromObject(fill_value, typenum, 0, 0)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(fill_value, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(typenum), 0, 0, (0x0100 | 0x0400) | 0x0040, ((void*)0)); | |||
| 1079 | if (afill == NULL((void*)0)) { | |||
| 1080 | /* For backwards compatibility try go via complex */ | |||
| 1081 | PyArrayObject *tmp; | |||
| 1082 | PyErr_Clear(); | |||
| 1083 | tmp = (PyArrayObject *)PyArray_FromObject(fill_value,(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(fill_value, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(NPY_CDOUBLE), 0, 0, (0x0100 | 0x0400) | 0x0040, ((void*)0)) | |||
| 1084 | NPY_CDOUBLE, 0, 0)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(fill_value, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(NPY_CDOUBLE), 0, 0, (0x0100 | 0x0400) | 0x0040, ((void*)0)); | |||
| 1085 | if (tmp == NULL((void*)0)) goto fail; | |||
| 1086 | afill = (PyArrayObject *)PyArray_Cast(tmp, typenum)(*(PyObject * (*)(PyArrayObject *, PyArray_Descr *, int)) _scipy_signal_ARRAY_API [49])(tmp, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API [45])(typenum), 0); | |||
| 1087 | Py_DECREF(tmp)_Py_DECREF(((PyObject*)(tmp))); | |||
| 1088 | if (afill == NULL((void*)0)) goto fail; | |||
| 1089 | /* Deprecated 2017-07, scipy version 1.0.0 */ | |||
| 1090 | if (DEPRECATE("could not cast `fillvalue` directly to the output "PyErr_WarnEx(PyExc_DeprecationWarning,"could not cast `fillvalue` directly to the output " "type (it was first converted to complex). " "This is deprecated and will raise an error in the " "future.",1) | |||
| 1091 | "type (it was first converted to complex). "PyErr_WarnEx(PyExc_DeprecationWarning,"could not cast `fillvalue` directly to the output " "type (it was first converted to complex). " "This is deprecated and will raise an error in the " "future.",1) | |||
| 1092 | "This is deprecated and will raise an error in the "PyErr_WarnEx(PyExc_DeprecationWarning,"could not cast `fillvalue` directly to the output " "type (it was first converted to complex). " "This is deprecated and will raise an error in the " "future.",1) | |||
| 1093 | "future.")PyErr_WarnEx(PyExc_DeprecationWarning,"could not cast `fillvalue` directly to the output " "type (it was first converted to complex). " "This is deprecated and will raise an error in the " "future.",1) < 0) { | |||
| 1094 | goto fail; | |||
| 1095 | } | |||
| 1096 | } | |||
| 1097 | if (PyArray_SIZE(afill)(*(npy_intp (*)(npy_intp const *, int)) _scipy_signal_ARRAY_API [158])(PyArray_DIMS(afill), PyArray_NDIM(afill)) != 1) { | |||
| 1098 | if (PyArray_SIZE(afill)(*(npy_intp (*)(npy_intp const *, int)) _scipy_signal_ARRAY_API [158])(PyArray_DIMS(afill), PyArray_NDIM(afill)) == 0) { | |||
| 1099 | PyErr_SetString(PyExc_ValueError, | |||
| 1100 | "`fillvalue` cannot be an empty array."); | |||
| 1101 | goto fail; | |||
| 1102 | } | |||
| 1103 | /* Deprecated 2017-07, scipy version 1.0.0 */ | |||
| 1104 | if (DEPRECATE("`fillvalue` must be scalar or an array with "PyErr_WarnEx(PyExc_DeprecationWarning,"`fillvalue` must be scalar or an array with " "one element. " "This will raise an error in the future.",1) | |||
| 1105 | "one element. "PyErr_WarnEx(PyExc_DeprecationWarning,"`fillvalue` must be scalar or an array with " "one element. " "This will raise an error in the future.",1) | |||
| 1106 | "This will raise an error in the future.")PyErr_WarnEx(PyExc_DeprecationWarning,"`fillvalue` must be scalar or an array with " "one element. " "This will raise an error in the future.",1) < 0) { | |||
| 1107 | goto fail; | |||
| 1108 | } | |||
| 1109 | } | |||
| 1110 | } | |||
| 1111 | else { | |||
| 1112 | /* Create a zero filled array */ | |||
| 1113 | afill = (PyArrayObject *)PyArray_ZEROS(0, NULL, typenum, 0)(*(PyObject * (*)(int, npy_intp const *, PyArray_Descr *, int )) _scipy_signal_ARRAY_API[183])(0, ((void*)0), (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(typenum), 0); | |||
| 1114 | if (afill == NULL((void*)0)) goto fail; | |||
| 1115 | } | |||
| 1116 | ||||
| 1117 | aout_dimens = malloc(PyArray_NDIM(ain1)*sizeof(npy_intp)); | |||
| 1118 | if (aout_dimens == NULL((void*)0)) goto fail; | |||
| 1119 | switch(mode & OUTSIZE_MASK3) { | |||
| 1120 | case VALID0: | |||
| 1121 | for (i = 0; i < PyArray_NDIM(ain1); i++) { | |||
| 1122 | aout_dimens[i] = PyArray_DIMS(ain1)[i] - PyArray_DIMS(ain2)[i] + 1; | |||
| 1123 | if (aout_dimens[i] < 0) { | |||
| 1124 | PyErr_SetString(PyExc_ValueError, | |||
| 1125 | "no part of the output is valid, use option 1 (same) or 2 " | |||
| 1126 | "(full) for third argument"); | |||
| 1127 | goto fail; | |||
| 1128 | } | |||
| 1129 | } | |||
| 1130 | break; | |||
| 1131 | case SAME1: | |||
| 1132 | for (i = 0; i < PyArray_NDIM(ain1); i++) { | |||
| 1133 | aout_dimens[i] = PyArray_DIMS(ain1)[i]; | |||
| 1134 | } | |||
| 1135 | break; | |||
| 1136 | case FULL2: | |||
| 1137 | for (i = 0; i < PyArray_NDIM(ain1); i++) { | |||
| 1138 | aout_dimens[i] = PyArray_DIMS(ain1)[i] + PyArray_DIMS(ain2)[i] - 1; | |||
| 1139 | } | |||
| 1140 | break; | |||
| 1141 | default: | |||
| 1142 | PyErr_SetString(PyExc_ValueError, | |||
| 1143 | "mode must be 0 (valid), 1 (same), or 2 (full)"); | |||
| 1144 | goto fail; | |||
| 1145 | } | |||
| 1146 | ||||
| 1147 | aout = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ain1), aout_dimens,(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int, npy_intp const *, void *, int, int, PyObject *)) _scipy_signal_ARRAY_API [93])(&(*(PyTypeObject *)_scipy_signal_ARRAY_API[2]), PyArray_NDIM (ain1), aout_dimens, typenum, ((void*)0), ((void*)0), 0, 0, ( (void*)0)) | |||
| 1148 | typenum)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int, npy_intp const *, void *, int, int, PyObject *)) _scipy_signal_ARRAY_API [93])(&(*(PyTypeObject *)_scipy_signal_ARRAY_API[2]), PyArray_NDIM (ain1), aout_dimens, typenum, ((void*)0), ((void*)0), 0, 0, ( (void*)0)); | |||
| 1149 | if (aout == NULL((void*)0)) goto fail; | |||
| 1150 | ||||
| 1151 | flag = mode + boundary + (typenum << TYPE_SHIFT5) + \ | |||
| 1152 | (flip != 0) * FLIP_MASK16; | |||
| 1153 | ||||
| 1154 | ret = pylab_convolve_2d (PyArray_DATA(ain1), /* Input data Ns[0] x Ns[1] */ | |||
| 1155 | PyArray_STRIDES(ain1), /* Input strides */ | |||
| 1156 | PyArray_DATA(aout), /* Output data */ | |||
| 1157 | PyArray_STRIDES(aout), /* Output strides */ | |||
| 1158 | PyArray_DATA(ain2), /* coefficients in filter */ | |||
| 1159 | PyArray_STRIDES(ain2), /* coefficients strides */ | |||
| 1160 | PyArray_DIMS(ain2), /* Size of kernel Nwin[2] */ | |||
| 1161 | PyArray_DIMS(ain1), /* Size of image Ns[0] x Ns[1] */ | |||
| 1162 | flag, /* convolution parameters */ | |||
| 1163 | PyArray_DATA(afill)); /* fill value */ | |||
| 1164 | ||||
| 1165 | ||||
| 1166 | switch (ret) { | |||
| 1167 | case 0: | |||
| 1168 | free(aout_dimens); | |||
| 1169 | Py_DECREF(ain1)_Py_DECREF(((PyObject*)(ain1))); | |||
| 1170 | Py_DECREF(ain2)_Py_DECREF(((PyObject*)(ain2))); | |||
| 1171 | Py_XDECREF(afill)_Py_XDECREF(((PyObject*)(afill))); | |||
| 1172 | return (PyObject *)aout; | |||
| 1173 | break; | |||
| 1174 | case -5: | |||
| 1175 | case -4: | |||
| 1176 | PyErr_SetString(PyExc_ValueError, | |||
| 1177 | "convolve2d not available for this type."); | |||
| 1178 | goto fail; | |||
| 1179 | case -3: | |||
| 1180 | PyErr_NoMemory(); | |||
| 1181 | goto fail; | |||
| 1182 | case -2: | |||
| 1183 | PyErr_SetString(PyExc_ValueError, | |||
| 1184 | "Invalid boundary type."); | |||
| 1185 | goto fail; | |||
| 1186 | case -1: | |||
| 1187 | PyErr_SetString(PyExc_ValueError, | |||
| 1188 | "Invalid output flag."); | |||
| 1189 | goto fail; | |||
| 1190 | } | |||
| 1191 | ||||
| 1192 | fail: | |||
| 1193 | free(aout_dimens); | |||
| 1194 | Py_XDECREF(ain1)_Py_XDECREF(((PyObject*)(ain1))); | |||
| 1195 | Py_XDECREF(ain2)_Py_XDECREF(((PyObject*)(ain2))); | |||
| 1196 | Py_XDECREF(aout)_Py_XDECREF(((PyObject*)(aout))); | |||
| 1197 | Py_XDECREF(afill)_Py_XDECREF(((PyObject*)(afill))); | |||
| 1198 | return NULL((void*)0); | |||
| 1199 | } | |||
| 1200 | ||||
| 1201 | /*******************************************************************/ | |||
| 1202 | ||||
| 1203 | static char doc_order_filterND[] = "out = _order_filterND(a,domain,order)"; | |||
| 1204 | ||||
| 1205 | static PyObject *sigtools_order_filterND(PyObject *NPY_UNUSED(dummy)(__NPY_UNUSED_TAGGEDdummy) __attribute__ ((__unused__)), | |||
| 1206 | PyObject *args) { | |||
| 1207 | PyObject *domain, *a0; | |||
| 1208 | int order=0; | |||
| 1209 | ||||
| 1210 | if (!PyArg_ParseTuple(args, "OO|i", &a0, &domain, &order)) return NULL((void*)0); | |||
| 1211 | ||||
| 1212 | return PyArray_OrderFilterND(a0, domain, order); | |||
| 1213 | } | |||
| 1214 | ||||
| 1215 | ||||
| 1216 | static char doc_remez[] = | |||
| 1217 | "h = _remez(numtaps, bands, des, weight, type, fs, maxiter, grid_density)\n" | |||
| 1218 | " returns the optimal (in the Chebyshev/minimax sense) FIR filter impulse\n" | |||
| 1219 | " response given a set of band edges, the desired response on those bands,\n" | |||
| 1220 | " and the weight given to the error in those bands. Bands is a monotonic\n" | |||
| 1221 | " vector with band edges given in frequency domain where fs is the sampling\n" | |||
| 1222 | " frequency."; | |||
| 1223 | ||||
| 1224 | static PyObject *sigtools_remez(PyObject *NPY_UNUSED(dummy)(__NPY_UNUSED_TAGGEDdummy) __attribute__ ((__unused__)), PyObject *args) | |||
| 1225 | { | |||
| 1226 | PyObject *bands, *des, *weight; | |||
| 1227 | int k, numtaps, numbands, type = BANDPASS1, err; | |||
| 1228 | PyArrayObject *a_bands=NULL((void*)0), *a_des=NULL((void*)0), *a_weight=NULL((void*)0); | |||
| 1229 | PyArrayObject *h=NULL((void*)0); | |||
| 1230 | npy_intp ret_dimens; int maxiter = 25, grid_density = 16; | |||
| 1231 | double oldvalue, *dptr, fs = 1.0; | |||
| 1232 | char mystr[255]; | |||
| 1233 | int niter = -1; | |||
| 1234 | ||||
| 1235 | if (!PyArg_ParseTuple(args, "iOOO|idii", &numtaps, &bands, &des, &weight, | |||
| 1236 | &type, &fs, &maxiter, &grid_density)) { | |||
| 1237 | return NULL((void*)0); | |||
| 1238 | } | |||
| 1239 | ||||
| 1240 | if (type != BANDPASS1 && type != DIFFERENTIATOR2 && type != HILBERT3) { | |||
| 1241 | PyErr_SetString(PyExc_ValueError, | |||
| 1242 | "The type must be BANDPASS, DIFFERENTIATOR, or HILBERT."); | |||
| 1243 | return NULL((void*)0); | |||
| 1244 | } | |||
| 1245 | ||||
| 1246 | if (numtaps < 2) { | |||
| 1247 | PyErr_SetString(PyExc_ValueError, | |||
| 1248 | "The number of taps must be greater than 1."); | |||
| 1249 | return NULL((void*)0); | |||
| 1250 | } | |||
| 1251 | ||||
| 1252 | a_bands = (PyArrayObject *)PyArray_ContiguousFromObject(bands, NPY_DOUBLE,1,1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(bands, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(NPY_DOUBLE), 1, 1, ( (0x0001 | (0x0100 | 0x0400))) | 0x0040, ((void*)0)); | |||
| 1253 | if (a_bands == NULL((void*)0)) goto fail; | |||
| 1254 | a_des = (PyArrayObject *)PyArray_ContiguousFromObject(des, NPY_DOUBLE,1,1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(des, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(NPY_DOUBLE), 1, 1, ( (0x0001 | (0x0100 | 0x0400))) | 0x0040, ((void*)0)); | |||
| 1255 | if (a_des == NULL((void*)0)) goto fail; | |||
| 1256 | a_weight = (PyArrayObject *)PyArray_ContiguousFromObject(weight, NPY_DOUBLE,1,1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(weight, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(NPY_DOUBLE), 1, 1, ( (0x0001 | (0x0100 | 0x0400))) | 0x0040, ((void*)0)); | |||
| 1257 | if (a_weight == NULL((void*)0)) goto fail; | |||
| 1258 | ||||
| 1259 | numbands = PyArray_DIMS(a_des)[0]; | |||
| 1260 | if ((PyArray_DIMS(a_bands)[0] != 2*numbands) || | |||
| 1261 | (PyArray_DIMS(a_weight)[0] != numbands)) { | |||
| 1262 | PyErr_SetString(PyExc_ValueError, | |||
| 1263 | "The inputs desired and weight must have same length.\n " | |||
| 1264 | "The input bands must have twice this length."); | |||
| 1265 | goto fail; | |||
| 1266 | } | |||
| 1267 | ||||
| 1268 | /* | |||
| 1269 | * Check the bands input to see if it is monotonic, divide by | |||
| 1270 | * fs to take from range 0 to 0.5 and check to see if in that range | |||
| 1271 | */ | |||
| 1272 | dptr = (double *)PyArray_DATA(a_bands); | |||
| 1273 | oldvalue = 0; | |||
| 1274 | for (k=0; k < 2*numbands; k++) { | |||
| 1275 | if (*dptr < oldvalue) { | |||
| 1276 | PyErr_SetString(PyExc_ValueError, | |||
| 1277 | "Bands must be monotonic starting at zero."); | |||
| 1278 | goto fail; | |||
| 1279 | } | |||
| 1280 | if (*dptr * 2 > fs) { | |||
| 1281 | PyErr_SetString(PyExc_ValueError, | |||
| 1282 | "Band edges should be less than 1/2 the sampling frequency"); | |||
| 1283 | goto fail; | |||
| 1284 | } | |||
| 1285 | oldvalue = *dptr; | |||
| 1286 | *dptr = oldvalue / fs; /* Change so that sampling frequency is 1.0 */ | |||
| 1287 | dptr++; | |||
| 1288 | } | |||
| 1289 | ||||
| 1290 | ret_dimens = numtaps; | |||
| 1291 | h = (PyArrayObject *)PyArray_SimpleNew(1, &ret_dimens, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int, npy_intp const *, void *, int, int, PyObject *)) _scipy_signal_ARRAY_API [93])(&(*(PyTypeObject *)_scipy_signal_ARRAY_API[2]), 1, & ret_dimens, NPY_DOUBLE, ((void*)0), ((void*)0), 0, 0, ((void* )0)); | |||
| 1292 | if (h == NULL((void*)0)) goto fail; | |||
| 1293 | ||||
| 1294 | err = pre_remez((double *)PyArray_DATA(h), numtaps, numbands, | |||
| 1295 | (double *)PyArray_DATA(a_bands), | |||
| 1296 | (double *)PyArray_DATA(a_des), | |||
| 1297 | (double *)PyArray_DATA(a_weight), | |||
| 1298 | type, maxiter, grid_density, &niter); | |||
| 1299 | if (err < 0) { | |||
| 1300 | if (err == -1) { | |||
| 1301 | sprintf(mystr, "Failure to converge at iteration %d, try reducing transition band width.\n", niter)__builtin___sprintf_chk (mystr, 2 - 1, __builtin_object_size ( mystr, 2 > 1), "Failure to converge at iteration %d, try reducing transition band width.\n" , niter); | |||
| 1302 | PyErr_SetString(PyExc_ValueError, mystr); | |||
| 1303 | goto fail; | |||
| 1304 | } | |||
| 1305 | else if (err == -2) { | |||
| 1306 | PyErr_NoMemory(); | |||
| 1307 | goto fail; | |||
| 1308 | } | |||
| 1309 | } | |||
| 1310 | ||||
| 1311 | Py_DECREF(a_bands)_Py_DECREF(((PyObject*)(a_bands))); | |||
| 1312 | Py_DECREF(a_des)_Py_DECREF(((PyObject*)(a_des))); | |||
| 1313 | Py_DECREF(a_weight)_Py_DECREF(((PyObject*)(a_weight))); | |||
| 1314 | ||||
| 1315 | return PyArray_Return(*(PyObject * (*)(PyArrayObject *)) _scipy_signal_ARRAY_API[76 ])(h); | |||
| 1316 | ||||
| 1317 | fail: | |||
| 1318 | Py_XDECREF(a_bands)_Py_XDECREF(((PyObject*)(a_bands))); | |||
| 1319 | Py_XDECREF(a_des)_Py_XDECREF(((PyObject*)(a_des))); | |||
| 1320 | Py_XDECREF(a_weight)_Py_XDECREF(((PyObject*)(a_weight))); | |||
| 1321 | Py_XDECREF(h)_Py_XDECREF(((PyObject*)(h))); | |||
| 1322 | return NULL((void*)0); | |||
| 1323 | } | |||
| 1324 | ||||
| 1325 | static char doc_median2d[] = "filt = _median2d(data, size)"; | |||
| 1326 | ||||
| 1327 | extern void f_medfilt2(float*,float*,npy_intp*,npy_intp*); | |||
| 1328 | extern void d_medfilt2(double*,double*,npy_intp*,npy_intp*); | |||
| 1329 | extern void b_medfilt2(unsigned char*,unsigned char*,npy_intp*,npy_intp*); | |||
| 1330 | ||||
| 1331 | static PyObject *sigtools_median2d(PyObject *NPY_UNUSED(dummy)(__NPY_UNUSED_TAGGEDdummy) __attribute__ ((__unused__)), PyObject *args) | |||
| 1332 | { | |||
| 1333 | PyObject *image=NULL((void*)0), *size=NULL((void*)0); | |||
| 1334 | int typenum; | |||
| 1335 | PyArrayObject *a_image=NULL((void*)0), *a_size=NULL((void*)0); | |||
| 1336 | PyArrayObject *a_out=NULL((void*)0); | |||
| 1337 | npy_intp Nwin[2] = {3,3}; | |||
| 1338 | ||||
| 1339 | if (!PyArg_ParseTuple(args, "O|O", &image, &size)) return NULL((void*)0); | |||
| 1340 | ||||
| 1341 | typenum = PyArray_ObjectType(*(int (*)(PyObject *, int)) _scipy_signal_ARRAY_API[54])(image, 0); | |||
| 1342 | a_image = (PyArrayObject *)PyArray_ContiguousFromObject(image, typenum, 2, 2)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(image, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(typenum), 2, 2, ((0x0001 | (0x0100 | 0x0400))) | 0x0040, ((void*)0)); | |||
| 1343 | if (a_image == NULL((void*)0)) goto fail; | |||
| 1344 | ||||
| 1345 | if (size != NULL((void*)0)) { | |||
| 1346 | a_size = (PyArrayObject *)PyArray_ContiguousFromObject(size, NPY_INTP, 1, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int, PyObject *)) _scipy_signal_ARRAY_API[69])(size, (*(PyArray_Descr * (*)(int)) _scipy_signal_ARRAY_API[45])(NPY_LONG), 1, 1, (( 0x0001 | (0x0100 | 0x0400))) | 0x0040, ((void*)0)); | |||
| 1347 | if (a_size == NULL((void*)0)) goto fail; | |||
| 1348 | if ((PyArray_NDIM(a_size) != 1) || (PyArray_DIMS(a_size)[0] < 2)) | |||
| 1349 | PYERR("Size must be a length two sequence"){PyErr_SetString(PyExc_ValueError, "Size must be a length two sequence" ); goto fail;}; | |||
| 1350 | Nwin[0] = ((npy_intp *)PyArray_DATA(a_size))[0]; | |||
| 1351 | Nwin[1] = ((npy_intp *)PyArray_DATA(a_size))[1]; | |||
| 1352 | } | |||
| 1353 | ||||
| 1354 | a_out = (PyArrayObject *)PyArray_SimpleNew(2, PyArray_DIMS(a_image), typenum)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int, npy_intp const *, void *, int, int, PyObject *)) _scipy_signal_ARRAY_API [93])(&(*(PyTypeObject *)_scipy_signal_ARRAY_API[2]), 2, PyArray_DIMS (a_image), typenum, ((void*)0), ((void*)0), 0, 0, ((void*)0)); | |||
| 1355 | if (a_out == NULL((void*)0)) goto fail; | |||
| 1356 | ||||
| 1357 | if (setjmp(MALLOC_FAIL)_setjmp (MALLOC_FAIL)) { | |||
| 1358 | PYERR("Memory allocation error."){PyErr_SetString(PyExc_ValueError, "Memory allocation error." ); goto fail;}; | |||
| 1359 | } | |||
| 1360 | else { | |||
| 1361 | switch (typenum) { | |||
| 1362 | case NPY_UBYTE: | |||
| 1363 | b_medfilt2((unsigned char *)PyArray_DATA(a_image), | |||
| 1364 | (unsigned char *)PyArray_DATA(a_out), | |||
| 1365 | Nwin, PyArray_DIMS(a_image)); | |||
| 1366 | break; | |||
| 1367 | case NPY_FLOAT: | |||
| 1368 | f_medfilt2((float *)PyArray_DATA(a_image), | |||
| 1369 | (float *)PyArray_DATA(a_out), Nwin, | |||
| 1370 | PyArray_DIMS(a_image)); | |||
| 1371 | break; | |||
| 1372 | case NPY_DOUBLE: | |||
| 1373 | d_medfilt2((double *)PyArray_DATA(a_image), | |||
| 1374 | (double *)PyArray_DATA(a_out), Nwin, | |||
| 1375 | PyArray_DIMS(a_image)); | |||
| 1376 | break; | |||
| 1377 | default: | |||
| 1378 | PYERR("2D median filter only supports uint8, float32, and float64."){PyErr_SetString(PyExc_ValueError, "2D median filter only supports uint8, float32, and float64." ); goto fail;}; | |||
| 1379 | } | |||
| 1380 | } | |||
| 1381 | ||||
| 1382 | Py_DECREF(a_image)_Py_DECREF(((PyObject*)(a_image))); | |||
| 1383 | Py_XDECREF(a_size)_Py_XDECREF(((PyObject*)(a_size))); | |||
| 1384 | ||||
| 1385 | return PyArray_Return(*(PyObject * (*)(PyArrayObject *)) _scipy_signal_ARRAY_API[76 ])(a_out); | |||
| 1386 | ||||
| 1387 | fail: | |||
| 1388 | Py_XDECREF(a_image)_Py_XDECREF(((PyObject*)(a_image))); | |||
| 1389 | Py_XDECREF(a_size)_Py_XDECREF(((PyObject*)(a_size))); | |||
| 1390 | Py_XDECREF(a_out)_Py_XDECREF(((PyObject*)(a_out))); | |||
| 1391 | return NULL((void*)0); | |||
| 1392 | ||||
| 1393 | } | |||
| 1394 | ||||
| 1395 | static char doc_linear_filter[] = | |||
| 1396 | "(y,Vf) = _linear_filter(b,a,X,Dim=-1,Vi=None) " \ | |||
| 1397 | "implemented using Direct Form II transposed flow " \ | |||
| 1398 | "diagram. If Vi is not given, Vf is not returned."; | |||
| 1399 | ||||
| 1400 | static struct PyMethodDef toolbox_module_methods[] = { | |||
| 1401 | {"_correlateND", scipy_signal_sigtools_correlateND, METH_VARARGS0x0001, doc_correlateND}, | |||
| 1402 | {"_convolve2d", sigtools_convolve2d, METH_VARARGS0x0001, doc_convolve2d}, | |||
| 1403 | {"_order_filterND", sigtools_order_filterND, METH_VARARGS0x0001, doc_order_filterND}, | |||
| 1404 | {"_linear_filter", scipy_signal_sigtools_linear_filter, METH_VARARGS0x0001, doc_linear_filter}, | |||
| 1405 | {"_remez",sigtools_remez, METH_VARARGS0x0001, doc_remez}, | |||
| 1406 | {"_medfilt2d", sigtools_median2d, METH_VARARGS0x0001, doc_median2d}, | |||
| 1407 | {NULL((void*)0), NULL((void*)0), 0, NULL((void*)0)} /* sentinel */ | |||
| 1408 | }; | |||
| 1409 | ||||
| 1410 | static struct PyModuleDef moduledef = { | |||
| 1411 | PyModuleDef_HEAD_INIT{ { 1, ((void*)0) }, ((void*)0), 0, ((void*)0), }, | |||
| 1412 | "sigtools", | |||
| 1413 | NULL((void*)0), | |||
| 1414 | -1, | |||
| 1415 | toolbox_module_methods, | |||
| 1416 | NULL((void*)0), | |||
| 1417 | NULL((void*)0), | |||
| 1418 | NULL((void*)0), | |||
| 1419 | NULL((void*)0) | |||
| 1420 | }; | |||
| 1421 | PyObject *PyInit_sigtools(void) | |||
| 1422 | { | |||
| 1423 | PyObject *m; | |||
| 1424 | ||||
| 1425 | m = PyModule_Create(&moduledef)PyModule_Create2(&moduledef, 1013); | |||
| ||||
| ||||
| 1426 | import_array(){if (_import_array() < 0) {PyErr_Print(); PyErr_SetString( PyExc_ImportError, "numpy.core.multiarray failed to import"); return ((void*)0); } }; | |||
| 1427 | ||||
| 1428 | scipy_signal_sigtools_linear_filter_module_init(); | |||
| 1429 | ||||
| 1430 | return m; | |||
| 1431 | } |
| 1 | #ifndef PyModule_Create2 |
| 2 | struct _object; |
| 3 | typedef struct _object PyObject; |
| 4 | PyObject* clang_analyzer_PyObject_New_Reference(); |
| 5 | PyObject* PyModule_Create2(PyModuleDef *def, int module_api_version) { |
| 6 | return clang_analyzer_PyObject_New_Reference(); |
| 7 | } |
| 8 | #else |
| 9 | #warning "API PyModule_Create2 is defined as a macro." |
| 10 | #endif |