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
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;}
16
17
18 jmp_buf MALLOC_FAIL;
19
check_malloc(size_t size)20 char *check_malloc(size_t size)
21 {
22 char *the_block = malloc(size);
23 if (the_block == NULL) {
24 printf("\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
index_out_of_bounds(npy_intp * indices,npy_intp * max_indices,int ndims)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
compute_offsets(npy_uintp * offsets,npy_intp * offsets2,npy_intp * dim1,npy_intp * dim2,npy_intp * dim3,npy_intp * mode_dep,int nd)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 */
increment(npy_intp * ret_ind,int nd,npy_intp * max_ind)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 BANDPASS 1
144 #define DIFFERENTIATOR 2
145 #define HILBERT 3
146
147 #define GOBACK goto
148 #define DOloop(a,from,to) for ( (a) = (from); (a) <= (to); ++(a))
149 #define PI 3.14159265358979323846
150 #define TWOPI (PI+PI)
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 */
lagrange_interp(int k,int n,int m,double * x)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) {
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 */
freq_eval(int k,int n,double * grid,double * x,double * y,double * ad)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 * grid[k]);
190
191 DOloop(j,1,n) {
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 */
remez(double * dev,double des[],double grid[],double edge[],double wt[],int ngrid,int nbands,int iext[],double alpha[],int nfcns,int itrmax,double * work,int dimsize,int * niter_out)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) {
250 x[j] = cos(grid[iext[j]]*TWOPI);
251 }
252 jet = (nfcns-1) / 15 + 1;
253
254 DOloop(j,1,nz) {
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) {
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) {
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 GOBACK L210;
316
317 L215:
318 iext[j++] = l - 1;
319 klow = l - 1;
320 ++jchnge;
321 GOBACK 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 GOBACK L235;
340 L240:
341 klow = iext[j];
342 iext[j] = l+1;
343 ++j;
344 ++jchnge;
345 GOBACK L200;
346
347 L250:
348 l = iext[j]+1;
349 if (jchnge > 0) GOBACK 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) GOBACK L255;
355 comp = (double)nut * err;
356
357 GOBACK L210;
358 L260:
359 klow = iext[j++];
360 GOBACK 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) GOBACK L310;
376 comp = (double) nut * err;
377 j = nzz;
378 GOBACK 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) GOBACK L330;
397 j = nzz;
398 comp = (double) nut * err;
399 luck = luck + 10;
400 GOBACK L235;
401 L340:
402 if (luck == 6) goto L370;
403 DOloop(j,1,nfcns) {
404 iext[nzz-j] = iext[nz-j];
405 }
406 iext[1] = k1;
407 GOBACK L100;
408 L350:
409 kn = iext[nzz];
410 DOloop(j,1,nfcns) iext[j] = iext[j+1];
411 iext[nz] = kn;
412
413 GOBACK 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*grid[1]);
436 dnum = cos(TWOPI*grid[ngrid]);
437 aa = 2.0/(dtemp-dnum);
438 bb = -(dtemp+dnum)/(dtemp-dnum);
439 }
440
441 DOloop(j,1,nfcns) {
442 ft = (j - 1) * delf;
443 xt = cos(TWOPI*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;
450 #else
451 ft = acos(xt)/TWOPI;
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 GOBACK L410;
460 L415:
461 a[j] = y[l];
462 goto L425;
463 L420:
464 if ((xt-xe) < fsh) GOBACK 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 / cn;
473 DOloop (j,1,nfcns) {
474 dtemp = 0.0;
475 dnum = (j-1) * dden;
476 if (nm1 >= 1) {
477 DOloop(k,1,nm1) {
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) 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) {
492 if (j >= nm1) {
493 aa *= 0.5;
494 bb *= 0.5;
495 }
496 p[j+1] = 0.0;
497 DOloop(k,1,j) {
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) p[k] += q[k] + aa * a[k+1];
504 jp1 = j + 1;
505 DOloop(k,3,jp1) p[k] += aa * a[k-1];
506
507 if (j != nm1) {
508 DOloop(k,1,j) q[k] = -a[k];
509 q[1] += alpha[nfcns - 1 - j];
510 }
511 }
512 DOloop(j,1,nfcns) 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 */
eff(double freq,double * fx,int lband,int jtype)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 */
wate(double freq,double * fx,double * wtx,int lband,int jtype)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
pre_remez(double * h2,int numtaps,int numbands,double * bands,double * response,double * weight,int type,int maxiter,int grid_density,int * niter_out)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) 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) {
670 change = cos(PI*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) {
678 change = sin(PI*grid[j]);
679 des[j] = des[j] / change;
680 wt[j] = wt[j] * change;
681 }
682 } else {
683 DOloop(j,1,ngrid) {
684 change = sin(TWOPI * 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) {
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) {
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) {
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) {
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) {
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){
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
fill_buffer(char * ip1,PyArrayObject * ap1,PyArrayObject * ap2,char * sort_buffer,int nels2,int check,npy_intp * loop_ind,npy_intp * temp_ind,npy_uintp * offset)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(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(ptr);
805 return;
806 }
807
808 #define COMPARE(fname, type) \
809 int fname(type *ip1, type *ip2) { return *ip1 < *ip2 ? -1 : *ip1 == *ip2 ? 0 : 1; }
810
COMPARE(DOUBLE_compare,double)811 COMPARE(DOUBLE_compare, double)
812 COMPARE(FLOAT_compare, float)
813 COMPARE(LONGDOUBLE_compare, npy_longdouble)
814 COMPARE(BYTE_compare, npy_byte)
815 COMPARE(SHORT_compare, short)
816 COMPARE(INT_compare, int)
817 COMPARE(LONG_compare, long)
818 COMPARE(LONGLONG_compare, npy_longlong)
819 COMPARE(UBYTE_compare, npy_ubyte)
820 COMPARE(USHORT_compare, npy_ushort)
821 COMPARE(UINT_compare, npy_uint)
822 COMPARE(ULONG_compare, npy_ulong)
823 COMPARE(ULONGLONG_compare, npy_ulonglong)
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_LT) == 1)
829 return -1;
830 else if(PyObject_RichCompareBool(*ip1, *ip2, Py_EQ) == 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, (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, NULL, NULL,
846 (CompareFunction)OBJECT_compare, NULL, NULL, NULL};
847
PyArray_OrderFilterND(PyObject * op1,PyObject * op2,int order)848 PyObject *PyArray_OrderFilterND(PyObject *op1, PyObject *op2, int order) {
849 PyArrayObject *ap1=NULL, *ap2=NULL, *ret=NULL;
850 npy_intp *a_ind=NULL, *b_ind=NULL, *temp_ind=NULL, *mode_dep=NULL, *check_ind=NULL;
851 npy_uintp *offsets=NULL;
852 npy_intp *offsets2=NULL;
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;
858 npy_intp *ret_ind=NULL;
859 CompareFunction compare_func=NULL;
860 char *zptr=NULL;
861 PyArray_CopySwapFunc *copyswap;
862
863 /* Get Array objects from input */
864 typenum = PyArray_ObjectType(op1, 0);
865 typenum = PyArray_ObjectType(op2, typenum);
866
867 ap1 = (PyArrayObject *)PyArray_ContiguousFromObject(op1, typenum, 0, 0);
868 if (ap1 == NULL) return NULL;
869 ap2 = (PyArrayObject *)PyArray_ContiguousFromObject(op2, typenum, 0, 0);
870 if (ap2 == NULL) 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((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(ap2);
886 if (zptr == NULL) 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),
899 PyArray_DIMS(ap1),
900 typenum);
901 if (ret == NULL) 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) {
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) 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) 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 || ret_ind == NULL) 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 || check_ind == NULL || offsets == NULL || offsets2 == NULL) 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(zptr);
962 zptr = PyArray_Zero(ap1);
963 if (zptr == NULL) 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((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);
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(zptr);
1019 Py_DECREF(ap1);
1020 Py_DECREF(ap2);
1021
1022 return PyArray_Return(ret);
1023
1024 fail:
1025 if (zptr) PyDataMem_FREE(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);
1036 Py_XDECREF(ap2);
1037 Py_XDECREF(ret);
1038 return NULL;
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
sigtools_convolve2d(PyObject * NPY_UNUSED (dummy),PyObject * args)1053 static PyObject *sigtools_convolve2d(PyObject *NPY_UNUSED(dummy), PyObject *args) {
1054
1055 PyObject *in1=NULL, *in2=NULL, *fill_value=NULL;
1056 int mode=2, boundary=0, typenum, flag, flip=1, ret;
1057 npy_intp *aout_dimens=NULL;
1058 int i;
1059 PyArrayObject *ain1=NULL, *ain2=NULL, *aout=NULL;
1060 PyArrayObject *afill=NULL;
1061
1062 if (!PyArg_ParseTuple(args, "OO|iiiO", &in1, &in2, &flip, &mode, &boundary,
1063 &fill_value)) {
1064 return NULL;
1065 }
1066
1067 typenum = PyArray_ObjectType(in1, 0);
1068 typenum = PyArray_ObjectType(in2, typenum);
1069 ain1 = (PyArrayObject *)PyArray_FromObject(in1, typenum, 2, 2);
1070 if (ain1 == NULL) goto fail;
1071 ain2 = (PyArrayObject *)PyArray_FromObject(in2, typenum, 2, 2);
1072 if (ain2 == NULL) goto fail;
1073
1074 if ((boundary != PAD) && (boundary != REFLECT) && (boundary != CIRCULAR))
1075 PYERR("Incorrect boundary value.");
1076
1077 if ((boundary == PAD) & (fill_value != NULL)) {
1078 afill = (PyArrayObject *)PyArray_FromObject(fill_value, typenum, 0, 0);
1079 if (afill == NULL) {
1080 /* For backwards compatibility try go via complex */
1081 PyArrayObject *tmp;
1082 PyErr_Clear();
1083 tmp = (PyArrayObject *)PyArray_FromObject(fill_value,
1084 NPY_CDOUBLE, 0, 0);
1085 if (tmp == NULL) goto fail;
1086 afill = (PyArrayObject *)PyArray_Cast(tmp, typenum);
1087 Py_DECREF(tmp);
1088 if (afill == NULL) goto fail;
1089 /* Deprecated 2017-07, scipy version 1.0.0 */
1090 if (DEPRECATE("could not cast `fillvalue` directly to the output "
1091 "type (it was first converted to complex). "
1092 "This is deprecated and will raise an error in the "
1093 "future.") < 0) {
1094 goto fail;
1095 }
1096 }
1097 if (PyArray_SIZE(afill) != 1) {
1098 if (PyArray_SIZE(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 "
1105 "one element. "
1106 "This will raise an error in the future.") < 0) {
1107 goto fail;
1108 }
1109 }
1110 }
1111 else {
1112 /* Create a zero filled array */
1113 afill = (PyArrayObject *)PyArray_ZEROS(0, NULL, typenum, 0);
1114 if (afill == NULL) goto fail;
1115 }
1116
1117 aout_dimens = malloc(PyArray_NDIM(ain1)*sizeof(npy_intp));
1118 if (aout_dimens == NULL) goto fail;
1119 switch(mode & OUTSIZE_MASK) {
1120 case VALID:
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 SAME:
1132 for (i = 0; i < PyArray_NDIM(ain1); i++) {
1133 aout_dimens[i] = PyArray_DIMS(ain1)[i];
1134 }
1135 break;
1136 case FULL:
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,
1148 typenum);
1149 if (aout == NULL) goto fail;
1150
1151 flag = mode + boundary + (typenum << TYPE_SHIFT) + \
1152 (flip != 0) * FLIP_MASK;
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);
1170 Py_DECREF(ain2);
1171 Py_XDECREF(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);
1195 Py_XDECREF(ain2);
1196 Py_XDECREF(aout);
1197 Py_XDECREF(afill);
1198 return NULL;
1199 }
1200
1201 /*******************************************************************/
1202
1203 static char doc_order_filterND[] = "out = _order_filterND(a,domain,order)";
1204
sigtools_order_filterND(PyObject * NPY_UNUSED (dummy),PyObject * args)1205 static PyObject *sigtools_order_filterND(PyObject *NPY_UNUSED(dummy),
1206 PyObject *args) {
1207 PyObject *domain, *a0;
1208 int order=0;
1209
1210 if (!PyArg_ParseTuple(args, "OO|i", &a0, &domain, &order)) return NULL;
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
sigtools_remez(PyObject * NPY_UNUSED (dummy),PyObject * args)1224 static PyObject *sigtools_remez(PyObject *NPY_UNUSED(dummy), PyObject *args)
1225 {
1226 PyObject *bands, *des, *weight;
1227 int k, numtaps, numbands, type = BANDPASS, err;
1228 PyArrayObject *a_bands=NULL, *a_des=NULL, *a_weight=NULL;
1229 PyArrayObject *h=NULL;
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;
1238 }
1239
1240 if (type != BANDPASS && type != DIFFERENTIATOR && type != HILBERT) {
1241 PyErr_SetString(PyExc_ValueError,
1242 "The type must be BANDPASS, DIFFERENTIATOR, or HILBERT.");
1243 return NULL;
1244 }
1245
1246 if (numtaps < 2) {
1247 PyErr_SetString(PyExc_ValueError,
1248 "The number of taps must be greater than 1.");
1249 return NULL;
1250 }
1251
1252 a_bands = (PyArrayObject *)PyArray_ContiguousFromObject(bands, NPY_DOUBLE,1,1);
1253 if (a_bands == NULL) goto fail;
1254 a_des = (PyArrayObject *)PyArray_ContiguousFromObject(des, NPY_DOUBLE,1,1);
1255 if (a_des == NULL) goto fail;
1256 a_weight = (PyArrayObject *)PyArray_ContiguousFromObject(weight, NPY_DOUBLE,1,1);
1257 if (a_weight == NULL) 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);
1292 if (h == NULL) 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);
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);
1312 Py_DECREF(a_des);
1313 Py_DECREF(a_weight);
1314
1315 return PyArray_Return(h);
1316
1317 fail:
1318 Py_XDECREF(a_bands);
1319 Py_XDECREF(a_des);
1320 Py_XDECREF(a_weight);
1321 Py_XDECREF(h);
1322 return NULL;
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
sigtools_median2d(PyObject * NPY_UNUSED (dummy),PyObject * args)1331 static PyObject *sigtools_median2d(PyObject *NPY_UNUSED(dummy), PyObject *args)
1332 {
1333 PyObject *image=NULL, *size=NULL;
1334 int typenum;
1335 PyArrayObject *a_image=NULL, *a_size=NULL;
1336 PyArrayObject *a_out=NULL;
1337 npy_intp Nwin[2] = {3,3};
1338
1339 if (!PyArg_ParseTuple(args, "O|O", &image, &size)) return NULL;
1340
1341 typenum = PyArray_ObjectType(image, 0);
1342 a_image = (PyArrayObject *)PyArray_ContiguousFromObject(image, typenum, 2, 2);
1343 if (a_image == NULL) goto fail;
1344
1345 if (size != NULL) {
1346 a_size = (PyArrayObject *)PyArray_ContiguousFromObject(size, NPY_INTP, 1, 1);
1347 if (a_size == NULL) goto fail;
1348 if ((PyArray_NDIM(a_size) != 1) || (PyArray_DIMS(a_size)[0] < 2))
1349 PYERR("Size must be a length two sequence");
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);
1355 if (a_out == NULL) goto fail;
1356
1357 if (setjmp(MALLOC_FAIL)) {
1358 PYERR("Memory allocation error.");
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.");
1379 }
1380 }
1381
1382 Py_DECREF(a_image);
1383 Py_XDECREF(a_size);
1384
1385 return PyArray_Return(a_out);
1386
1387 fail:
1388 Py_XDECREF(a_image);
1389 Py_XDECREF(a_size);
1390 Py_XDECREF(a_out);
1391 return NULL;
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_VARARGS, doc_correlateND},
1402 {"_convolve2d", sigtools_convolve2d, METH_VARARGS, doc_convolve2d},
1403 {"_order_filterND", sigtools_order_filterND, METH_VARARGS, doc_order_filterND},
1404 {"_linear_filter", scipy_signal_sigtools_linear_filter, METH_VARARGS, doc_linear_filter},
1405 {"_remez",sigtools_remez, METH_VARARGS, doc_remez},
1406 {"_medfilt2d", sigtools_median2d, METH_VARARGS, doc_median2d},
1407 {NULL, NULL, 0, NULL} /* sentinel */
1408 };
1409
1410 static struct PyModuleDef moduledef = {
1411 PyModuleDef_HEAD_INIT,
1412 "sigtools",
1413 NULL,
1414 -1,
1415 toolbox_module_methods,
1416 NULL,
1417 NULL,
1418 NULL,
1419 NULL
1420 };
PyInit_sigtools(void)1421 PyObject *PyInit_sigtools(void)
1422 {
1423 PyObject *m;
1424
1425 m = PyModule_Create(&moduledef);
1426 import_array();
1427
1428 scipy_signal_sigtools_linear_filter_module_init();
1429
1430 return m;
1431 }
1432