1 /*
2 * $Id: filter.i,v 1.1 2005-09-18 22:05:56 dhmunro Exp $
3 * analog signal processing routines
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 func filter(filt, dt, signal, pad=, shift=)
12 /* DOCUMENT filter(filt, dt, signal)
13
14 apply the filter FILT to the input SIGNAL, which is sampled
15 at times spaced by DT. The filter is assumed to be normalized
16 to an angular frequency (e.g.- radians per second), unless
17 DT<0, in which case FILT is assumed to be normalized to a
18 circular frequency (e.g.- Hz or GHz).
19
20 The result will have the same length as SIGNAL; be sure to pad
21 SIGNAL if you need the response to go beyond that time, or
22 you can use the pad=n keyword to force the returned result to
23 have N samples more than SIGNAL.
24
25 If the shift= keyword is non-nil and non-0, then the result
26 is shifted backward in time by the filter group delay at
27 zero frequency.
28
29 The impulse response of the FILT is also assumed to be shorter
30 than the duration of signal, and SIGNAL is assumed to be sampled
31 finely enough to resolve the FILT impulse response.
32
33 FILT is an array of double, which represents a filter with
34 a particular finite list of zeroes and poles. See the specific
35 functions to construct filters from poles and zeroes (fil_make),
36 or classic Bessel, Butterworth, Chebyshev, inverse Chebyshev, or
37 Cauer (elliptic) designs. With fil_analyze, you can find the
38 poles and zeroes of a FILT. The format for FILT is:
39
40 FILT is an array of double with the following meanings:
41 FILT(1) = np = number of poles (integer >= 0)
42 FILT(2) = nz = number of zeroes (integer >= 0)
43 FILT(3) = reserved
44 FILT(4:4+nz) = coefficients for numerator
45 = [a0, a1, a2, a3, ..., anz]
46 FILT(5+nz:4+nz+np) = coefficents for denominator (if np>0)
47 = [b1, b2, b3, ..., bnp]
48
49 The Laplace transform (s-transform) of the filter response is
50
51 L[FILT] = (a0 + a1*s + a2*s^2 + a3*s^3 + ...) /
52 ( 1 + b1*s + b2*s^2 + b3*s^3 + ...)
53
54 SEE ALSO: filter, fil_bessel, fil_butter, fil_cheby1, fil_cheby2,
55 fil_cauer, fil_response, fil_make, fil_analyze,
56 to_db, to_phase
57 */
58 {
59 m= numberof(signal);
60 if (is_void(pad)) pad= 0;
61 n= fft_good(max(2*m,m+pad));
62 signal= fft(grow(double(signal), array(0.,n-m)))/n;
63 if (dt>0.) w= (2*pi/n/dt)*roll(indgen(n/2+1-n:n/2),n/2+1);
64 else w= (-1./n/dt)*roll(indgen(n/2+1-n:n/2),n/2+1);
65 signal= double(fft(signal*fil_response(filt,w),-1));
66 if (shift) {
67 d= long(fil_delay(filt,dt<0.)/abs(dt));
68 signal= signal(d:0);
69 if (numberof(signal)<m+pad)
70 grow, signal, array(0.,m+pad-numberof(signal));
71 }
72 return signal(1:m+pad);
73 }
74
fil_response(filt,w)75 func fil_response(filt, w)
76 /* DOCUMENT fil_response(filt, w)
77
78 return the complex response of FILT at the frequencies W.
79 The frequency scale for W depends on how FILT has been scaled;
80 filters are rational functions in W.
81
82 The to_db and to_phase functions may be useful for extracting
83 the attenuation and phase parts of the complex response.
84
85 SEE ALSO: filter, fil_butter, fil_bessel, fil_cheby1, fil_cheby2,
86 fil_delay, to_db, to_phase
87 */
88 {
89 np= long(filt(1));
90 nz= long(filt(2));
91 w*= 1i;
92 y= np? fil_poly(filt(5+nz:4+nz+np),w) : 0.0;
93 return fil_poly(filt(4:4+nz),w)/(1.+y*w);
94 }
95
fil_delay(filt,hz)96 func fil_delay(filt, hz)
97 /* DOCUMENT fil_delay(filt)
98 or fil_delay(filt, 1)
99
100 return the group delay d(phase)/dw at w=0 (zero frequency) for
101 filter FILT. By default, FILT is assumed to be normalized
102 to an angular frequency (e.g.- radians per second), but if
103 the 2nd parameter is non-nil and non-0 FILT is assumed to be
104 normalized to a circular frequency (e.g.- Hz or GHz).
105
106 SEE ALSO: filter, fil_butter, fil_bessel, fil_cheby1, fil_cheby2,
107 fil_response, to_db, to_phase
108 */
109 {
110 np= long(filt(1));
111 nz= long(filt(2));
112 dpdw= np? filt(5+nz) : 0.0;
113 if (nz) dpdw-= filt(5)/filt(4);
114 if (hz) dpdw/= 2.*pi;
115 return dpdw;
116 }
117
fil_poly(c,x)118 func fil_poly(c, x)
119 /* DOCUMENT fil_poly(c, x)
120 return c(1) + c(2)*x + c(3)*x^2 + c(4)*x^3 + ...
121 */
122 {
123 n= numberof(c);
124 y= array(c(n), dimsof(x));
125 while(--n) y= y*x+c(n);
126 return y;
127 }
128
fil_make(poles,zeroes)129 func fil_make(poles, zeroes)
130 /* DOCUMENT filt= fil_make(poles, zeroes)
131
132 given the complex POLES and ZEROES, return a FILT. The real
133 parts of POLES must all be negative to make a stable FILT.
134 Both POLES and ZEROES must occur in conjugate pairs in order to
135 make a real filter (the returned filter is always real).
136
137 The returned filter always has a0=1 (its DC gain is 1).
138
139 SEE ALSO: filter, fil_analyze
140 */
141 {
142 np= numberof(poles);
143 nz= numberof(zeroes);
144 filt= array(0., 4+nz+np);
145 filt(1:4)= [np,nz,0,1];
146 if (nz) {
147 a= array(0i, nz+1);
148 a(0)= 1.0;
149 for (i=1 ; i<=nz ; i++) a(1:-1)-= a(2:0)*zeroes(i);
150 a/= double(a(1));
151 filt(4:4+nz)= a; /* discards imaginary part */
152 }
153 if (np) {
154 a= array(0i, np+1);
155 a(0)= 1.0;
156 for (i=1 ; i<=np ; i++) a(1:-1)-= a(2:0)*poles(i);
157 a= a(2:0)/double(a(1));
158 filt(5+nz:4+nz+np)= a; /* discards imaginary part */
159 }
160 return filt;
161 }
162
163 func fil_analyze(filt, &poles, &zeroes)
164 /* DOCUMENT fil_analyze, filt, poles, zeroes
165
166 given a FILT, return the complex POLES and ZEROES, sorted in
167 order of increasing imaginary part. The real parts of POLES will
168 all be negative if the FILT is stable.
169
170 SEE ALSO: filter, fil_make
171 */
172 {
173 require, "zroots.i";
174 np= long(filt(1));
175 nz= long(filt(2));
176 poles= zeroes= [];
177 if (nz>0) zeroes= zroots(filt(4:4+nz), imsort=1);
178 if (np>0) poles= zroots(grow([1.],filt(5+nz:4+nz+np)), imsort=1);
179 }
180
181 func fil_bessel(np, wc, db, natural=)
182 /* DOCUMENT filt= fil_bessel(np, wc, db)
183
184 returns the lowpass Bessel filter with NP poles, normalized
185 such that at angular frequency WC, the attenuation is DB decibels.
186 (That is, the attenuation factor is 10^(.05*DB) at WC,
187 so that to_db(response(filt,WC))==DB.)
188
189 A Bessel filter has the most nearly constant group delay time
190 d(phase)/dw of any filter of the same order. It minimizes pulse
191 distortion, but does not cut off very rapidly in frequency.
192
193 If WC is nil or zero, it defaults to 1.0.
194
195 If DB is nil, the filter is normalized such that both the s^0
196 and s^NP terms are 1, unless the natural= keyword is non-zero,
197 in which case the filter is normalized such that the group delay
198 d(phase)/dw is -1 at w=0.
199
200 SEE ALSO: filter, fil_analyze
201 */
202 {
203 filt= array(0., 4+np);
204 filt(1:4)= [np,0,0,1];
205 if (np) {
206 ca= cb= array(0., np+1);
207 ca(1)= 1.;
208 cb(1:2)= c= [1.,1.];
209 for (i=2 ; i<=np ; i++,ca=cb,cb=c) c= (2.*i-1.)*cb + roll(ca,2);
210 /* several different normalizations are used:
211 * (0) note c(1) = (2n-1)!! (product of odd numbers <= 2*np-1)
212 * (1) c/= c(1) normalizes to a group delay -d(phase)/dw = 1.0
213 * for all orders (coefficient c(2)=c(1)=1)
214 * (2) c/= c(1)^span(1.,0.,np+1) normalizes both leading and trailing
215 * polynomial coefficients to 1.0
216 * corresponds to w->w*((2n-1)!!)^(1/n)
217 * limit as n->big is w->w*(2*n/e^2), so
218 * phase delay becomes proportional to order n
219 * (3) scale w so that polynomial value at w=1 is sqrt(0.5) (3 db)
220 * (or some other value) independent of order
221 */
222 if (!natural) c/= c(1)^span(1.,0.,np+1);
223 else c/= c(1);
224 filt(4:4+np)= c;
225 }
226 return fil_normalize(filt, wc, db);
227 }
228
fil_butter(np,wc,db)229 func fil_butter(np, wc, db)
230 /* DOCUMENT filt= fil_butter(np, wc, db)
231
232 returns the lowpass Butterworth filter with NP poles, normalized
233 such that at angular frequency WC, the attenuation is DB decibels.
234 (That is, the attenuation factor is 10^(.05*DB) at WC,
235 so that to_db(response(filt,WC))==DB.)
236
237 A Butterworth filter is the best Taylor series approximation to
238 the ideal lowpass filter (a step in frequency) response at both
239 w=0 and w=infinity.
240
241 For wc=1 and db=10*log10(2), the square of the Butterworth frequency
242 response is 1/(1+w^(2*np)).
243
244 If WC is nil or zero, it defaults to 1.0.
245
246 If DB is nil, the filter is normalized "naturally", which is the
247 same as DB=10*log10(2).
248
249 SEE ALSO: filter, fil_analyze, butter
250 */
251 {
252 filt= array(0., 4+np);
253 filt(1:4)= [np,0,0,1];
254 /* poles: 1i*exp(0.5i*pi/np*(2*indgen(np)-1)) */
255 a= pi/(2*np);
256 b= cos(a*indgen(0:np-1))/sin(a*indgen(np));
257 for (i=1,a=1.0 ; i<=np ; i++) filt(4+i)= a= a*b(i);
258 if (db) {
259 if (!wc) wc= 1.0;
260 wc*= (10.^(0.1*db)-1.)^(-0.5/np);
261 }
262 if (wc) filt(5:4+np)/= wc^indgen(np);
263 return filt;
264 }
265
fil_cheby1(np,ripple,wc,db)266 func fil_cheby1(np, ripple, wc, db)
267 /* DOCUMENT filt= fil_cheby1(np, ripple, wc, db)
268
269 returns the lowpass Chebyshev type I filter with NP poles, and
270 passband ripple RIPPLE decibels, normalized such that at
271 angular frequency WC, the attenuation is DB decibels.
272 (That is, the attenuation factor is 10^(.05*DB) at WC,
273 so that to_db(response(filter,WC))==DB.)
274
275 A Chebyshev type I filter gives the smallest maximum error over the
276 passband for any filter that is a Taylor series approximation to
277 the ideal lowpass filter (a step in frequency) response at
278 w=infinity. It has NP/2 ripples of amplitude RIPPLE in its passband,
279 and a smooth stopband.
280
281 For wc=1 and db=ripple, the square of the Chebyshev frequency
282 response is 1/(1+eps2*Tnp(w)), where eps2 = 10^(ripple/10)-1,
283 and Tnp is the np-th Chebyshev polynomial, cosh(np*acosh(x)) or
284 cos(np*acos(x)).
285
286 If WC is nil or zero, it defaults to 1.0.
287
288 If DB is nil, the filter is normalized "naturally", which is the
289 same as DB=RIPPLE.
290
291 SEE ALSO: filter, fil_analyze, cheby1
292 */
293 { /* passband ripple */
294 eps21= 10.^(0.1*ripple);
295 reps= 1./sqrt(eps21-1.);
296 ripple= asinh(reps)/np;
297 poles= 1i*cos(0.5*pi/np*(2*indgen(np)-1) - 1i*ripple);
298 filt= fil_make(poles);
299 if (!(np%2)) filt(4)/= sqrt(eps21);
300 if (db) {
301 if (!wc) wc= 1.0;
302 wc*= sech(acosh(reps*sqrt(10.^(0.1*db)-1.))/np);
303 }
304 if (wc) filt(5:4+np)/= wc^indgen(np);
305 return filt;
306 }
307
fil_cheby2(np,atten,wc,db)308 func fil_cheby2(np, atten, wc, db)
309 /* DOCUMENT filt= fil_cheby2(np, atten, wc, db)
310
311 returns the lowpass Chebyshev type II filter with NP poles, and
312 stopband attenuation ATTEN decibels, normalized such that at
313 angular frequency WC, the attenuation is DB decibels.
314 (That is, the attenuation factor is 10^(.05*DB) at WC,
315 so that to_db(response(filter,WC))==DB.)
316
317 This is also called an inverse Chebyshev filter, since its poles
318 are the reciprocals of a Chebyshev type I filter. It has NP zeroes
319 as well as NP poles.
320
321 A Chebyshev type II filter gives the smallest maximum leakage over
322 the stopband for any filter that is a Taylor series approximation to
323 the ideal lowpass filter (a step in frequency) response at
324 w=0. It has NP/2 ripples of amplitude ATTEN in its stopband,
325 and a smooth passband.
326
327 For wc=1 and db=ripple, the square of the inverse Chebyshev frequency
328 response is 1 - 1/(1+eps2*Tnp(1/w)), where eps2 = 10^(ripple/10)-1 =
329 1/(10^(atten/10)-1) and Tnp is the np-th Chebyshev polynomial,
330 cosh(np*acosh(x)) or cos(np*acos(x)).
331
332 If WC is nil or zero, it defaults to 1.0.
333
334 If DB is nil, the filter is normalized "naturally", which is the
335 same as DB=ATTEN.
336
337 SEE ALSO: filter, fil_analyze, cheby2
338 */
339 { /* stopband attenuation (ripple) */
340 reps= sqrt(10.^(0.1*atten)-1.);
341 atten= asinh(reps)/np;
342 wk= 0.5*pi/np*(2*indgen(np)-1);
343 poles= -1i/cos(wk - 1i*atten);
344 zeroes= 1i/cos(wk);
345 if (np==1) zeroes= [];
346 else if (np%2) zeroes= grow(zeroes(1:(np-1)/2),zeroes((np+3)/2:0));
347 filt= fil_make(poles, zeroes);
348 if (db) {
349 if (!wc) wc= 1.0;
350 wc*= cosh(acosh(reps/sqrt(10.^(0.1*db)-1.))/np);
351 }
352 if (wc) {
353 nz= numberof(zeroes);
354 if (nz) filt(5:4+nz)/= wc^indgen(nz);
355 filt(5+nz:4+nz+np)/= wc^indgen(np);
356 }
357 return filt;
358 }
359
fil_cauer(np,ripple,atten,wc,db)360 func fil_cauer(np, ripple, atten, wc, db)
361 /* DOCUMENT filt= fil_cauer(np, ripple, atten, wc, db)
362 or filt= fil_cauer(np, ripple, -skirt, wc, db)
363
364 returns the lowpass Cauer (elliptic) filter with NP poles, passband
365 ripple RIPPLE and stopband attenuation ATTEN decibels, normalized
366 such that at angular frequency WC, the attenuation is DB decibels.
367 (That is, the attenuation factor is 10^(.05*DB) at WC,
368 so that to_db(response(filter,WC))==DB.)
369
370 If the third parameter is negative, its absolute value is SKIRT,
371 the ratio of the frequency at which the stopband attenuation is
372 first reached to the frequency at which the passband ends (where
373 the attenuation is RIPPLE). The closer to 1.0 SKIRT is, the
374 smaller the equivalent ATTEN would be. The external variable
375 cauer_other is set to ATTEN if you provide SKIRT, and to SKIRT
376 if you provide ATTEN.
377
378 The Cauer filter has NP zeroes as well as NP poles.
379
380 Consider the four parameters: (1) filter order, (2) transition
381 ("skirt") bandwidth, (3) passband ripple, and (4) stopband ripple.
382 Given any three of these, the Cauer filter minimizes the fourth.
383
384 If WC is nil or zero, it defaults to 1.0.
385
386 If DB is nil, the filter is normalized "naturally", which is the
387 same as DB=RIPPLE.
388
389 SEE ALSO: filter, fil_analyze, cauer
390 */
391 {
392 extern cauer_other;
393 eps2= 10.^(0.1*ripple)-1.;
394 if (atten > 0.) {
395 mb= eps2/(10.^(0.1*atten)-1.);
396 ma= _cauer_parameter(np, mb);
397 cauer_other= 1./sqrt(ma);
398 } else {
399 ma= 1./(atten*atten);
400 mb= _cauer_parameter(1./np, ma);
401 cauer_other= 10.*log10(1.+eps2/mb);
402 }
403 ekb= ellip_k(mb);
404 eka= ellip_k(ma);
405 rat= np*ekb/eka; /* equals ellip_k(1-mb)/ellip_k(1-ma) */
406
407 if (np>1) {
408 zeroes= indgen(np-1:1:-2)(-:1:2,)(*);
409 zeroes(1:0:2)*= -1.;
410 zeroes= 1i/(sqrt(ma)*sin(ell_am(zeroes*eka/np,ma)));
411 }
412 a= ell_f(atan(1./sqrt(eps2)),1.-mb)/rat;
413 b= (2.*indgen(np)-(np+1))*eka/np;
414 phi= ell_am(b,ma);
415 d= dn_(phi);
416 s= sin(phi);
417 c= cos(phi);
418 phi= ell_am(a,1.-ma);
419 dp= dn_(phi);
420 sp= sin(phi);
421 cp= cos(phi);
422 poles= (-c*d*sp*cp + 1i*s*dp)/(cp*cp+ma*(s*sp)^2);
423
424 filt= fil_make(poles, zeroes);
425 if (db) {
426 if (!wc) wc= 1.0;
427 wc*= dn_(ell_am(ell_f(asin(sqrt((1.-eps2/(10.^(0.1*db)-1.))/(1.-mb))),
428 1.-mb) / rat, 1.-ma));
429 }
430 if (wc) {
431 nz= numberof(zeroes);
432 if (nz) filt(5:4+nz)/= wc^indgen(nz);
433 filt(5+nz:4+nz+np)/= wc^indgen(np);
434 }
435 return filt;
436 }
437
_cauer_parameter(n,mb)438 func _cauer_parameter(n, mb)
439 {
440 /* K'(mb)/K'(ma) = n*K(mb)/K(ma)
441 * where the rational Chebyshev polynomial is
442 * sn( n*K(mb)/K(ma) * asn(w|ma) | mb )
443 * <asn means inverse of sn>
444 *
445 * the passband ripple is determined by the additional
446 * parameter e, and the stopband (w>1) ripple is e/mb > e
447 *
448 * now K'(ma)/K(ma) = K'(mb)/K(mb) * (1/n)
449 * and q(m) = exp(-pi*K'(m)/K(m)) is called the "nome" of m,
450 * m ~ 16*q for small q
451 */
452 require, "elliptic.i";
453 if (!mb) return 0.;
454 if (mb==1.) return 1.;
455 if (mb < 1.e-4) u= -log(poly(mb/16., 0., 1., 8., 84., 992.))/pi;
456 else u= ellip_k(1.-mb)/ellip_k(mb);
457 u/= n; /* K'(ma)/K(ma) */
458 if (!u) return 1.0;
459 recip= (u<1.);
460 if (recip) u= 1./u;
461 q= qn= exp(-pi*u); /* < 1/20 */
462 q2n= num= 1.;
463 den= 0.5;
464 for (;;) {
465 den+= q2n*qn;
466 q2n*= qn*qn;
467 numx= num+q2n;
468 if (numx==num) break;
469 num= numx;
470 qn*= q;
471 }
472 ma= q*(num/den)^4;
473 return recip? (1.-ma) : ma;
474 }
475
fil_normalize(filt,wc,db)476 func fil_normalize(filt, wc, db)
477 /* xxDOCUMENT filt= fil_normalize(filt, wc, db)
478 normalize the all pole FILT to give DB decibels attenuation
479 at frequency WC, by rescaling the frequency.
480 Assumes w=1 is a good starting point for search for initial filt;
481 probably should not be used except for cases in this file.
482 */
483 {
484 if (filt(1)>1.5 && !filt(2)) {
485 np= long(filt(1));
486 if (db) {
487 c= filt(4:4+np); /* polynomial (denominator) */
488 c1= c(1);
489 c(1)= 1.0;
490 c/= c1;
491 d= c(2:0)*indgen(np); /* derivative */
492 db= 10.^(0.05*db); /* y value to be found */
493 x= 1.0; /* initial guess */
494 for (i=1 ; i<100 ; i++) {
495 py= fil_poly(c, 1i*x);
496 y= abs(py);
497 dydx= -(fil_poly(d, 1i*x)*conj(py)).im/y;
498 dx= (db-y)/dydx;
499 if (abs(dx)<1.e-12*x) break;
500 x+= dx;
501 }
502 if (i>=100) error, "failed to converge to specified dB level";
503 if (!wc) wc= 1.0;
504 wc/= x;
505 }
506 if (wc) filt(4:4+np)/= wc^indgen(0:np);
507 }
508 return filt;
509 }
510
511 /* ------------------------------------------------------------------------ */
512
to_db(signal,ref)513 func to_db(signal, ref)
514 /* DOCUMENT to_db(signal, ref)
515 or to_db(signal)
516
517 return 20.*log10(abs(SIGNAL)/REF), the number of decibels
518 corresponding to the input SIGNAL. REF defaults to 1.0.
519
520 SEE ALSO: fil_response, to_phase
521 */
522 {
523 if (is_void(ref)) ref= 1.0;
524 return 20.*log10(max(abs(signal)/ref,1.e-37));
525 }
526
to_phase(signal,degrees)527 func to_phase(signal, degrees)
528 /* DOCUMENT to_phase(signal)
529 or to_phase(signal, 1)
530
531 return atan(SIGNAL.im,SIGNAL.re), the phase of the input SIGNAL.
532 If the second argument is present and non-0, the phase will be in
533 degrees; by default the phase is in radians.
534
535 To_phase attempts to unroll any jumps from -180 to +180 degrees
536 or vice-versa; zero phase will be taken somewhere near the middle
537 of the signal. The external variable to_phase_eps controls the
538 details of this unrolling; you can turn off unrolling by setting
539 to_phase_eps=0.0 (initially it is 0.3).
540
541 SEE ALSO: fil_response, to_phase
542 */
543 {
544 p= atan(signal.im, signal.re+!abs(signal));
545 np= numberof(p);
546 if (np>2) {
547 p0= p(1:-1);
548 p1= p(2:0);
549 eps= (1.00000001-to_phase_eps)*pi;
550 up= where((p0>=eps)&(p1<=-eps));
551 dn= where((p0<=-eps)&(p1>=eps));
552 adj= 0.*p0;
553 if (numberof(up)) adj(up)= 2.*pi;
554 if (numberof(dn)) adj(dn)= -2.*pi;
555 p+= adj(cum) - adj(np/2+1);
556 }
557 if (degrees) p*= 180/pi;
558 return p;
559 }
560
561 to_phase_eps= 0.3;
562
563 /* ------------------------------------------------------------------------ */
564 /* the frequency responses of the Butterworth, Chebyshev, and Cauer
565 * filters have simple analytic forms
566 * - note the frequency scaling rules are also simple */
567
butter(np,w,wc,db)568 func butter(np,w,wc,db)
569 /* DOCUMENT butter(np, w)
570 or butter(np, w, wc, db)
571
572 return frequency response (amplitude) for Butterworth filter;
573 the parameters are the same as for fil_butter.
574
575 SEE ALSO: fil_butter
576 */
577 {
578 if (db) w*= (10.^(0.1*db)-1.)^(0.5/np);
579 if (wc) w/= wc;
580 return 1./sqrt(1.+w^(2*np));
581 }
582
cheby1(np,ripple,w,wc,db)583 func cheby1(np,ripple,w,wc,db)
584 /* DOCUMENT cheby1(np, ripple, w)
585 or cheby1(np, ripple, w, wc, db)
586
587 return frequency response (amplitude) for Chebyshev filter;
588 the parameters are the same as for fil_cheby1.
589
590 SEE ALSO: fil_cheby1
591 */
592 {
593 eps2= 10.^(0.1*ripple)-1.;
594 if (db) w*= cosh(acosh(sqrt((10.^(0.1*db)-1.)/eps2))/np);
595 if (wc) w/= wc;
596 /* 1/sqrt(1+eps2*Tn(w)^2) where Tn is nth Chebyshev polynomial */
597 mask= (w<=1.0);
598 list= where(mask);
599 if (numberof(list)) cn= 1./(1.+eps2*cos(np*acos(w(list)))^2);
600 list= where(!mask);
601 if (numberof(list)) {
602 cnw= sech(np*acosh(w(list)))^2;
603 cnw= cnw/(cnw+eps2);
604 }
605 return sqrt(merge(cn,cnw,mask));
606 }
607
cheby2(np,atten,w,wc,db)608 func cheby2(np,atten,w,wc,db)
609 /* DOCUMENT cheby2(np, atten, w)
610 or cheby2(np, atten, w, wc, db)
611
612 return frequency response (amplitude) for inverse Chebyshev filter;
613 the parameters are the same as for fil_cheby2.
614
615 SEE ALSO: fil_cheby2
616 */
617 {
618 eps2= 1./(10.^(0.1*atten)-1.);
619 if (db) w*= sech(acosh(1./sqrt((10.^(0.1*db)-1.)*eps2))/np);
620 if (wc) w/= wc;
621 /* sqrt(1-1/(1+eps2*Tn(1/w)^2)) where Tn is nth Chebyshev polynomial */
622 w= (1.+(!w)*1.e99)/(w+(!w));
623 mask= (w<=1.0);
624 list= where(mask);
625 if (numberof(list)) cn= 1./(1.+eps2*cos(np*acos(w(list)))^2);
626 list= where(!mask);
627 if (numberof(list)) {
628 cnw= sech(np*acosh(w(list)))^2;
629 cnw= cnw/(cnw+eps2);
630 }
631 return sqrt(1.-merge(cn,cnw,mask));
632 }
633
cauer(n,ripple,atten,w,wc,db)634 func cauer(n, ripple, atten, w, wc, db)
635 /* DOCUMENT cauer(np, ripple, atten, w)
636 or cauer(np, ripple, atten, w, wc, db)
637
638 return frequency response (amplitude) for Cauer filter;
639 the parameters are the same as for fil_cauer.
640
641 SEE ALSO: fil_cauer
642 */
643 {
644 eps2= 10.^(0.1*ripple)-1.;
645 if (atten > 0.) {
646 mb= eps2/(10.^(0.1*atten)-1.);
647 ma= _cauer_parameter(n, mb);
648 } else {
649 ma= 1./(atten*atten);
650 mb= _cauer_parameter(1./n, ma);
651 }
652 ekb= ellip_k(mb);
653 rat= n*ekb/ellip_k(ma);
654 if (db)
655 w/= dn_(ell_am(ell_f(asin(sqrt((1.-eps2/(10.^(0.1*db)-1.))/(1.-mb))),
656 1.-mb) / rat, 1.-ma));
657 if (wc) w/= wc;
658
659 w= abs(w);
660 mask= (w<=1.);
661 list= where(mask);
662 if (numberof(list)) {
663 x= rat*ell_f(asin(w(list)),ma);
664 if (!(n%2)) x+= ekb;
665 f1= sin(ell_am(x,mb));
666 }
667 list= where(!mask);
668 if (numberof(list)) {
669 w= 1./w(list);
670 ka= sqrt(ma);
671 mask2= (w>ka);
672 list= where(mask2);
673 if (numberof(list)) {
674 x= w(list);
675 x= rat*ell_f(asin(sqrt((1.-x*x)/(1.-ma))),1.-ma);
676 f21= 1./dn_(ell_am(x,1.-mb));
677 }
678 list= where(!mask2);
679 if (numberof(list)) {
680 x= rat*ell_f(asin(w(list)/ka),ma);
681 if (!(n%2)) x+= ekb;
682 f22= 1./(sqrt(mb)*sin(ell_am(x,mb)));
683 }
684 f2= merge(f21,f22,mask2);
685 }
686 f= merge(f1,f2,mask);
687
688 return 1./sqrt(1.+eps2*f*f);
689 }
690
691 /* ------------------------------------------------------------------------ */
692 /* also in convol.i */
693
fft_good(n)694 func fft_good(n)
695 /* DOCUMENT fft_good(n)
696
697 returns the smallest number of the form 2^x*3^y*5^z greater
698 than or equal to n. An fft of this length will be much faster
699 than a number with larger prime factors; the speed difference
700 can be an order of magnitude or more.
701
702 For n>100, the worst cases result in a little over a 11% increase
703 in n; for n>1000, the worst are a bit over 6%; still larger n are
704 better yet. The median increase for n<=10000 is about 1.5%.
705
706 SEE ALSO: fft, fft_setup, convol
707 */
708 {
709 if (n<7) return max(n,1);
710 logn= log(n);
711 n5= 5.^indgen(0:long(logn/log(5.) + 1.e-6)); /* exact integers */
712 n3= 3.^indgen(0:long(logn/log(3.) + 1.e-6)); /* exact integers */
713 n35= n3*n5(-,); /* fewer than 300 numbers for n<5e9 */
714 n35= n35(where(n35<=n));
715 n235= 2^long((logn-log(n35))/log(2.) + 0.999999) * long(n35);
716 return min(n235);
717 }
718
719 /* ------------------------------------------------------------------------ */
720