/[sapaclisp]/sapaclisp/utilities.lisp
ViewVC logotype

Contents of /sapaclisp/utilities.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (vendor branch)
Mon May 9 20:25:10 2005 UTC (8 years, 11 months ago) by mmommer
Branch: MAIN, T
CVS Tags: initial_checkin, HEAD
Changes since 1.1: +0 -0 lines
Initial Checkin
1 ;;;-*- Mode: LISP; Package: :SAPA; Syntax: COMMON-LISP -*-
2 ;-----------------------------------------------------------------------------
3 ; (c) 1993, Donald B. Percival <dbp@apl.washington.edu>
4 ;
5 ; This code is licensed under the terms of the modified BSD license
6 ; ("sans advertising clause"). See the file COPYING for details.
7 ;
8 ; Comments about this software can be addressed to dbp@apl.washington.edu
9 ;-----------------------------------------------------------------------------
10 ;
11 ; utilities.lisp
12 ;
13 ; a collection of Lisp functions for certain basic operations ...
14 ; Note: before compiling and loading utilities.lisp,
15 ; you should compile and load
16 ; sapa-package.lisp
17 ;
18 ;-----------------------------------------------------------------------------
19 (if (not (find-package :SAPA))
20 (defpackage :SAPA (:USE :COMMON-LISP)))
21
22 (in-package :SAPA)
23
24 (export '(;;; two macros like incf for multiplication and division ...
25 multf
26 divf
27
28 ;;; functions dealing with logs and decibels ...
29 convert-to-dB
30 careful-convert-to-dB
31 convert-from-dB
32 log10
33 power-of-2
34 next-power-of-2
35
36 ;;; functions to manipulate sequence(s) and return sequence(s) ...
37 a*x
38 a*x!
39 a*x+b
40 a*x+b!
41 a*x+y
42 a*x+y!
43 x+y
44 x+y!
45 x+b
46 x+b!
47 x-y
48 x-y!
49 x*y
50 x*y!
51 add-sequences
52 multiply-sequences
53 circular-shift-sequence
54 copy-vector
55 cumulative-sums
56 difference
57 transform-a-sequence
58 transform-a-sequence!
59 linear-interpolation!
60
61 ;;; functions to manipulate sequence(s) and return scalar(s) ...
62 dot-product
63 euclidean-norm
64 sum-of-squares
65 sum
66 min-and-max-of-seq
67 min-of-seq
68 max-of-seq
69 binary-search
70 compare-seqs
71
72 ;;; functions to generate sequence(s) ...
73 iota
74 sample-from-a-function
75
76 ;;; miscellaneous functions ...
77 sign
78 sampling-time->Nyquist-frequency
79 Nyquist-frequency->sampling-time
80 ))
81
82 ;-------------------------------------------------------------------------------
83 ;-------------------------------------------------------------------------------
84 ;;; The macros multf
85 ;;; divf
86 ;;; act like incf and decf in Common Lisp, but use multiplication
87 ;;; and division rather than addition and subtraction:
88 ;;; (multf x a) is somewhat like (setf x (* x a))
89 ;;; (divf x a) ... (setf x (/ x a))
90 ;;; These are the ONLY macros used in SAPA.
91 ;-------------------------------------------------------------------------------
92 ;-------------------------------------------------------------------------------
93 (define-modify-macro multf (&optional (delta 1)) * "Like <incf>")
94
95 (define-modify-macro divf (&optional (delta 1)) / "Like <incf>")
96
97 ;-------------------------------------------------------------------------------
98 ;;; The following constants are used with the next collection of functions:
99 (defconstant +sapa-10-over-log-10+ (/ 10.0d0 (log 10.0)))
100 (defconstant +sapa-2-pi+ (* 2 pi))
101 (defconstant +sapa-minus-2-pi+ (* -2 pi))
102
103 ;-------------------------------------------------------------------------------
104 ;-------------------------------------------------------------------------------
105 ;;; The functions convert-to-dB
106 ;;; convert-from-dB
107 ;;; log10
108 ;;; power-of-2
109 ;;; next-power-of-2
110 ;;; all take a single number and perform some sort of a ``log or exponent''
111 ;;; operation on it:
112 ;;; convert-to-dB converts the number to decibels;
113 ;;; convert-from-dB converts the number back from its decibel formulation;
114 ;;; log10 returns the log base 10 of the number
115 ;;; power-of-2 tests the number to see if it is a power of 2; and
116 ;;; next-power-of-2 returns a power of 2 greater than or equal to the number.
117 ;-------------------------------------------------------------------------------
118 ;-------------------------------------------------------------------------------
119 (defun convert-to-dB (num)
120 "given
121 [1] num (required)
122 ==> a number
123 returns
124 [1] 10 log_10(num)
125 ---
126 Note: num is intended to be a positive number,
127 but this is not checked"
128 (* (log num) +sapa-10-over-log-10+))
129
130 ;;; (convert-to-dB 0.0001) ;==> -39.999999999999986
131
132 ;-------------------------------------------------------------------------------
133 (defun careful-convert-to-dB (num &optional (zero-mapping -100.0))
134 "given
135 [1] num (required)
136 ==> a number
137 [1] zero-mapping (required)
138 ==> a number
139 returns
140 [1] 10 log_10(num) if num > 0 or
141 zero-mapping if num <= 0"
142 (if (plusp num)
143 (* (log num) +sapa-10-over-log-10+)
144 zero-mapping))
145
146 #|
147 (careful-convert-to-dB 0.0001)
148 ;==> -39.999999999999986
149 (careful-convert-to-dB 0.0)
150 ;==> -100.0
151 (careful-convert-to-dB 0.0 -1000.0)
152 ;==> -1000.0
153 (careful-convert-to-dB 0.0 nil)
154 ;==> nil
155 |#
156
157 ;-------------------------------------------------------------------------------
158 (defun convert-from-dB (num)
159 "given
160 [1] num (required)
161 ==> a number
162 returns
163 [1] 10^(x/10)"
164 (exp (/ num +sapa-10-over-log-10+)))
165
166 ;;; (convert-from-dB 20.0) ;==> 100.00000000000004
167
168 ;-------------------------------------------------------------------------------
169 (defun log10 (num)
170 "given
171 [1] num (required)
172 ==> a number
173 returns
174 [1] log_10(num)
175 ---
176 Note: shorthand for (log num 10)"
177 (log num 10))
178
179 ;;; (log10 100.0) ;==> 2.0
180
181 ;-------------------------------------------------------------------------------
182 (defun power-of-2 (n)
183 "given
184 [1] n (required)
185 ==> an integer
186 returns
187 [1] nil if n is not a power of 2
188 -- or --
189 m if n = 2^m for a positive integer m"
190 (let ((answer+1 (power-of-2+1 n)))
191 (if (numberp answer+1)
192 (1- answer+1))))
193
194 ;;; (power-of-2 1023) ;==> nil
195 ;;; (power-of-2 32) ;==> 5
196
197 ;-------------------------------------------------------------------------------
198 (defun next-power-of-2 (n)
199 "given
200 [1] n (required)
201 ==> an integer
202 returns
203 [1] an integer power of 2 greater than or
204 equal to n"
205 (expt 2 (ceiling (log n 2))))
206
207 ;;; (next-power-of-2 31) ;==> 32
208 ;;; (next-power-of-2 32) ;==> 32
209 ;;; (next-power-of-2 33) ;==> 64
210
211 ;-------------------------------------------------------------------------------
212 ;-------------------------------------------------------------------------------
213 ;;; The functions a*x
214 ;;; a*x!
215 ;;; a*x+b
216 ;;; a*x+b!
217 ;;; a*x+y
218 ;;; a*x+y!
219 ;;; x+y
220 ;;; x+y!
221 ;;; x+b
222 ;;; x+b!
223 ;;; x-y
224 ;;; x-y!
225 ;;; x*y
226 ;;; x*y!
227 ;;; add-sequences
228 ;;; multiply-sequences
229 ;;; circular-shift-sequence
230 ;;; copy-vector
231 ;;; cumulative-sums
232 ;;; difference
233 ;;; transform-a-sequence
234 ;;; transform-a-sequence!
235 ;;; linear-interpolation!
236 ;;; all take sequence(s) of numbers as input and return sequence(s)
237 ;;; of numbers as output. Function names that end with a ``!''
238 ;;; indicate that the function will modify an input sequence
239 ;;; UNLESS the function has a keyword parameter with the name ``result''
240 ;;; with a value supplied by the user.
241 ;-------------------------------------------------------------------------------
242 ;-------------------------------------------------------------------------------
243 (defun a*x (a x)
244 "given
245 [1] a (required)
246 ==> a number
247 [2] x (required)
248 ==> a sequence of numbers
249 returns
250 [1] a vector of numbers obtained by multiplying
251 each element of x by a
252 ----
253 Note: corresponds to SSCAL in LINPACK's bla."
254 (map 'vector #'(lambda (y) (* a y)) x))
255
256 #|
257 (setf *a* (vector 10.0 20.0 30.0))
258 (a*x pi *a*) ;==> #(31.41592653589793 62.83185307179586 94.24777960769379)
259 *a* ;==> #(10.0 20.0 30.0)
260 |#
261
262 ;-------------------------------------------------------------------------------
263 (defun a*x! (a x &key (result x))
264 "given
265 [1] a (required)
266 ==> a number
267 [2] x (required)
268 ==> a sequence of numbers
269 [2] result (keyword; x)
270 <== a sequence of numbers
271 returns
272 [1] a sequence of numbers obtained by multiplying
273 each element of x by a
274 ----
275 Note: corresponds to SSCAL in LINPACK's bla."
276 (map-into result #'(lambda (y) (* a y)) x))
277
278 #|
279 (setf *a* (vector 10.0 20.0 30.0))
280 (a*x! pi *a*) ;==> #(31.41592653589793 62.83185307179586 94.24777960769379)
281 *a* ;==> #(31.41592653589793 62.83185307179586 94.24777960769379)
282 (setf *a* (vector 10.0 20.0 30.0))
283 (setf *b* `(,nil ,nil ,nil))
284 (a*x! pi *a* :result *b*)
285 *a* ;==> #(10.0 20.0 30.0)
286 *b* ;==> (31.41592653589793 62.83185307179586 94.24777960769379)
287 |#
288
289 ;-------------------------------------------------------------------------------
290 (defun a*x+b (a x b)
291 "given
292 [1] a (required)
293 ==> a number
294 [2] x (required)
295 ==> a sequence of numbers
296 [3] b (required)
297 ==> a number
298 returns
299 [1] a vector of numbers obtained by multiplying
300 each element of x by a and adding b"
301 (map 'vector #'(lambda (y) (+ b (* a y))) x))
302
303 #|
304 (setf *x* (vector 1.0 2.0 3.0))
305 (a*x+b pi *x* 10.0)
306 ;==> #(13.141592653589793 16.283185307179586 19.42477796076938)
307 *x* ;==> #(1.0 2.0 3.0)
308 |#
309
310 ;-------------------------------------------------------------------------------
311 (defun a*x+b! (a x b &key (result x))
312 "given
313 [1] a (required)
314 ==> a number
315 [2] x (required)
316 ==> a sequence of numbers
317 [3] b (required)
318 ==> a number
319 [4] result (keyword; x)
320 <== a sequence of numbers
321 returns
322 [1] a sequence of numbers obtained by multiplying
323 each element of x by a and adding b"
324 (map-into result #'(lambda (y) (+ b (* a y))) x))
325
326 #|
327 (setf *x* (vector 1.0 2.0 3.0))
328 (a*x+b! pi *x* 10.0)
329 ;==> #(13.141592653589793 16.283185307179586 19.42477796076938)
330 *x* ;==> #(13.141592653589793 16.283185307179586 19.42477796076938)
331 |#
332
333 ;-------------------------------------------------------------------------------
334 (defun a*x+y (a x y)
335 "given
336 [1] a (required)
337 ==> a number
338 [2] x (required)
339 ==> a sequence of numbers
340 [3] y (required)
341 ==> a sequence of numbers
342 returns
343 [1] a vector of numbers obtained by multiplying
344 each element of x by a and adding the
345 corresponding element of y
346 ----
347 Note: corresponds to SAXPY in LINPACK's bla."
348 (map 'vector #'(lambda (u v) (+ (* a u) v)) x y))
349
350 #|
351 (setf *x* (vector 1.0 2.0 3.0))
352 (setf *y* (vector 10.0 20.0 30.0))
353 (a*x+y pi *x* *y*)
354 ;==> #(13.141592653589793 26.283185307179586 39.424777960769376)
355 *x* ;==> #(1.0 2.0 3.0)
356 *y* ;==> #(10.0 20.0 30.0)
357 |#
358
359 ;-------------------------------------------------------------------------------
360 (defun a*x+y! (a x y &key (result y))
361 "given
362 [1] a (required)
363 ==> a number
364 [2] x (required)
365 ==> a sequence of numbers
366 [3] y (required)
367 ==> a sequence of numbers
368 [4] result (keyword; y)
369 <== a sequence of numbers
370 returns
371 [1] a sequence of numbers obtained by multiplying
372 each element of x by a and adding the
373 corresponding element of y
374 ----
375 Note: corresponds to SAXPY in LINPACK's bla."
376 (map-into result #'(lambda (u v) (+ (* a u) v)) x y))
377
378 #|
379 (setf *x* (vector 1.0 2.0 3.0))
380 (setf *y* (vector 10.0 20.0 30.0))
381 (a*x+y! pi *x* *y*)
382 ;==> #(13.141592653589793 26.283185307179586 39.424777960769376)
383 *x* ;==> #(1.0 2.0 3.0)
384 *y* ;==> #(13.141592653589793 26.283185307179586 39.424777960769376)
385 |#
386
387 ;-------------------------------------------------------------------------------
388 (defun x+y (x y)
389 "given
390 [1] x (required)
391 ==> a sequence of numbers
392 [2] y (required)
393 ==> a sequence of numbers
394 returns
395 [1] a vector of numbers obtained by adding
396 x and y together element by element"
397 (map 'vector #'+ x y))
398
399 #|
400 (setf *x* (vector 1.0 2.0 3.0))
401 (setf *y* (vector 10.0 20.0 30.0))
402 (x+y *x* *y*)
403 ;==> #(11.0 22.0 33.0)
404 *x* ;==> #(1.0 2.0 3.0)
405 *y* ;==> #(10.0 20.0 30.0)
406 |#
407
408 ;-------------------------------------------------------------------------------
409 (defun x+y! (x y &key (result x))
410 "given
411 [1] x (required)
412 ==> a sequence of numbers
413 [2] y (required)
414 ==> a sequence of numbers
415 [3] result (keyword; x)
416 <== a sequence of numbers
417 returns
418 [1] a sequence of numbers obtained by adding
419 x and y together element by element"
420 (map-into result #'+ x y))
421
422 #|
423 (setf *x* (vector 1.0 2.0 3.0))
424 (setf *y* (vector 10.0 20.0 30.0))
425 (x+y! *x* *y*)
426 ;==> #(11.0 22.0 33.0)
427 *x* ;==> #(11.0 22.0 33.0)
428 *y* ;==> #(10.0 20.0 30.0)
429 |#
430
431 ;-------------------------------------------------------------------------------
432 (defun x+b (x b)
433 "given
434 [1] x (required)
435 ==> a sequence of numbers
436 [2] b (required)
437 ==> a number
438 returns
439 [1] a vector of numbers obtained by adding
440 b to each element of x"
441 (map 'vector #'(lambda (y) (+ b y)) x))
442
443 #|
444 (setf *x* (vector 1.0 2.0 3.0))
445 (x+b *x* 10.0)
446 ;==> #(11.0 12.0 13.0)
447 *x* ;==> #(1.0 2.0 3.0)
448 |#
449
450 ;-------------------------------------------------------------------------------
451 (defun x+b! (x b &key (result x))
452 "given
453 [1] x (required)
454 ==> a sequence of numbers
455 [2] b (required)
456 ==> a number
457 [3] result (keyword; x)
458 <== a sequence of numbers
459 returns
460 [1] a sequence of numbers obtained by adding
461 b to each element of x"
462 (map-into result #'(lambda (y) (+ b y)) x))
463
464 #|
465 (setf *x* (vector 1.0 2.0 3.0))
466 (x+b! *x* 10.0)
467 ;==> #(11.0 12.0 13.0)
468 *x* ;==> #(11.0 12.0 13.0)
469 |#
470
471 ;-------------------------------------------------------------------------------
472 (defun x-y (x y)
473 "given
474 [1] x (required)
475 ==> a sequence of numbers
476 [2] y (required)
477 ==> a sequence of numbers
478 returns
479 [1] a vector of numbers obtained by subtracting
480 elements of y from corresponding elements of x"
481 (map 'vector #'- x y))
482
483 #|
484 (setf *x* (vector 1.0 2.0 3.0))
485 (setf *y* (vector 10.0 20.0 30.0))
486 (x-y *x* *y*)
487 ;==> #(-9.0 -18.0 -27.0)
488 *x* ;==> #(1.0 2.0 3.0)
489 *y* ;==> #(10.0 20.0 30.0)
490 |#
491
492 ;-------------------------------------------------------------------------------
493 (defun x-y! (x y &key (result x))
494 "given
495 [1] x (required)
496 ==> a sequence of numbers
497 [2] y (required)
498 ==> a sequence of numbers
499 [3] result (keyword; x)
500 <== a sequence of numbers
501 returns
502 [1] a sequence of numbers obtained by subtracting
503 elements of y from corresponding elements of x"
504 (map-into result #'- x y))
505
506 #|
507 (setf *x* (vector 1.0 2.0 3.0))
508 (setf *y* (vector 10.0 20.0 30.0))
509 (x-y! *x* *y*)
510 ;==> #(-9.0 -18.0 -27.0)
511 *x* ;==> #(-9.0 -18.0 -27.0)
512 *y* ;==> #(10.0 20.0 30.0)
513 |#
514
515 ;-------------------------------------------------------------------------------
516 (defun x*y (x y)
517 "given
518 [1] x (required)
519 ==> a sequence of numbers
520 [2] y (required)
521 ==> a sequence of numbers
522 returns
523 [1] a vector of numbers obtained by multiplying
524 corresponding elements of x and y"
525 (map 'vector #'* x y))
526
527 #|
528 (setf *x* (vector 1.0 2.0 3.0))
529 (setf *y* (vector 10.0 20.0 30.0))
530 (x*y *x* *y*)
531 ;==> #(10.0 40.0 90.0)
532 *x* ;==> #(1.0 2.0 3.0)
533 *y* ;==> #(10.0 20.0 30.0)
534 |#
535
536 ;-------------------------------------------------------------------------------
537 (defun x*y! (x y &key (result x))
538 "given
539 [1] x (required)
540 ==> a sequence of numbers
541 [2] y (required)
542 ==> a sequence of numbers
543 [3] result (keyword; x)
544 <== a sequence of numbers
545 returns
546 [1] a sequence of numbers obtained by multiplying
547 corresponding elements of x and y"
548 (map-into result #'* x y))
549
550 #|
551 (setf *x* (vector 1.0 2.0 3.0))
552 (setf *y* (vector 10.0 20.0 30.0))
553 (x*y! *x* *y*)
554 ;==> #(10.0 40.0 90.0)
555 *x* ;==> #(10.0 40.0 90.0)
556 *y* ;==> #(10.0 20.0 30.0)
557 |#
558
559 ;-------------------------------------------------------------------------------
560 (defun add-sequences (&rest seqs)
561 "given
562 [1] a arbitrary number of sequences of numbers,
563 all of the same size
564 returns
565 [1] a new sequence formed by adding the sequences
566 together on an element by element basis"
567 (apply #'map (type-of (car seqs)) #'+ seqs))
568
569 ;;; (add-sequences '(1 2 3) '(4 5 6) '(7 8 9)) ;==> (12 15 18)
570 ;;; (add-sequences #(1 2 3) #(4 5 6) #(7 8 9)) ;==> #(12 15 18)
571
572 ;-------------------------------------------------------------------------------
573 (defun multiply-sequences (&rest seqs)
574 "given
575 [1] a arbitrary number of sequences of numbers,
576 all of the same size
577 returns
578 [1] a new sequence formed by multiplying the sequences
579 together on an element by element basis"
580 (apply #'map (type-of (car seqs)) #'* seqs))
581
582 ;;; (multiply-sequences '(1 2 3) '(4 5 6) '(7 8 9)) ;==> (28 80 162)
583 ;;; (multiply-sequences #(1 2 3) #(4 5 6) #(7 8 9)) ;==> #(28 80 162)
584
585 ;-------------------------------------------------------------------------------
586 (defun circular-shift-sequence
587 (a-seq
588 &key
589 (result
590 (make-array (length a-seq))))
591 "given
592 [1] a-seq (required)
593 ==> a sequence
594 [2] result (keyword; vector of same length as a-seq)
595 <== a sequence
596 returns
597 [1] a `circularly' left-shifted sequence; i.e.,
598 maps x(0), x(1), x(2), ..., x(n-2), x(n-1)
599 to x(n-1), x(0), x(1), ..., x(n-3), x(n-2)
600 ----
601 Note: result can be the same sequence as result"
602 (let* ((n (length a-seq))
603 (n-1 (1- n))
604 (temp (elt a-seq n-1)))
605 (dotimes (i (1- n))
606 (setf (elt result n-1)
607 (elt a-seq (decf n-1))))
608 (setf (elt result 0) temp)
609 result))
610
611 #|
612 (setf *a* (vector 10.0 20.0 30.0 40.0 50.0 60.0 70.0))
613 (circular-shift-sequence *a* :result *a*)
614 ;==> #(70.0 10.0 20.0 30.0 40.0 50.0 60.0)
615 (circular-shift-sequence *a*)
616 ;==> #(60.0 70.0 10.0 20.0 30.0 40.0 50.0)
617 *a*
618 ;==> #(70.0 10.0 20.0 30.0 40.0 50.0 60.0)
619 |#
620
621 ;-------------------------------------------------------------------------------
622 ;;; The name of this function is a misleading in that it can accept
623 ;;; any type of sequences (not just vectors).
624 (defun copy-vector
625 (from-vector
626 to-vector
627 &key
628 (start 0)
629 (end (length from-vector)))
630 "given
631 [1] from-vector (required)
632 ==> a vector (actually, works with any sequence)
633 [2] to-vector (required)
634 ==> a vector (actually, works with any sequence)
635 [3] start (keyword; 0)
636 ==> start index of from-vector to be transferred
637 [4] end (keyword; length of from-vector)
638 ==> 1 + end index of from-vector to be transferred
639 returns
640 [1] to-vector, with elements start to end - 1
641 from-vector copied into corresponding
642 elements 0 to end - start - 1 of to-vector
643 ---
644 Note: if to-vector is being created on the spot,
645 might consider using CL's copy-seq instead."
646 (replace to-vector from-vector :start2 start :end2 end))
647
648 #|
649 (copy-vector #(1 2 3) (make-array 3))
650 ;==> #(1 2 3)
651 (copy-vector #(1 2 3) (make-array 6 :initial-element 0) :start 1)
652 ;==> #(2 3 0 0 0 0)
653 (copy-vector #(1 2 3) (make-array 6 :initial-element 0) :end 2)
654 ;==> #(1 2 0 0 0 0)
655 |#
656
657 ;-------------------------------------------------------------------------------
658 (defun cumulative-sums
659 (the-seq
660 &key
661 (start 0)
662 (end (length the-seq))
663 (result (make-array (- end start))))
664 "given
665 [1] the-seq (required)
666 ==> a sequence of numbers
667 [2] start (keyword; 0)
668 ==> start index of sequence to be used
669 [3] end (keyword; length of the-seq)
670 ==> 1 + end index of sequence to be used
671 returns
672 [1] sequence with cumulative sum
673 of elements in specified subsequence"
674 (let ((N (- end start)))
675 (if (plusp N) (setf (elt result 0) (elt the-seq start)))
676 (dotimes (i (1- N) result)
677 (setf (elt result (1+ i)) (+ (elt result i)
678 (elt the-seq (incf start)))))))
679
680 ;;; (cumulative-sums #(10 20 30 40 55)) ;==> #(10 30 60 100 155)
681 ;;; (cumulative-sums #(10 20 30 40 55) :start 2) ;==> #(30 70 125)
682
683 ;-------------------------------------------------------------------------------
684 (defun difference
685 (sequence-of-numbers
686 &key
687 (start 0)
688 (end (length sequence-of-numbers))
689 (lag 1)
690 (result (make-array (- end start lag))))
691 "given
692 [1] sequence-of-numbers (required)
693 ==> a sequence of numbers
694 [2] start (keyword; 0)
695 ==> start index
696 [3] end (keyword; length of sequence-of-numbers)
697 ==> end index plus one
698 [4] lag (keyword; 1)
699 ==> lag for differencing
700 [5] result (keyword; vector of length (- end start lag))
701 <== results
702 return
703 [1] a sequence of length (- end start lag)
704 with lag-th order difference of sequence
705 from index start to index end-1
706 ---
707 Note: See Newton, 1988, page 31."
708 (if (or (not (= start 0))
709 (not (= end (length sequence-of-numbers))))
710 (setf sequence-of-numbers (subseq sequence-of-numbers start end)))
711 (let ((k lag))
712 (dotimes (j (- end start lag) result)
713 (setf (elt result j) (- (elt sequence-of-numbers k)
714 (elt sequence-of-numbers j)))
715 (incf k))))
716
717 ;;; (difference #(1 2 4 5 7 9 11)) ;==> #(1 2 1 2 2 2)
718
719 ;-------------------------------------------------------------------------------
720 (defun transform-a-sequence (a-function the-seq)
721 "given
722 [1] a-function (required)
723 ==> a function with one argument
724 [2] the-seq (required)
725 ==> a sequence of numbers
726 returns
727 [1] a sequence of numbers obtained by applying
728 a-function to each element of the-seq"
729 (map (type-of the-seq) a-function the-seq))
730
731 ;;; (transform-a-sequence #'convert-to-dB '(1 10 100)) ;==> (0.0 10.0 20.0)
732 ;;; (transform-a-sequence #'convert-to-dB #(1 10 100)) ;==> #(0.0 10.0 20.0)
733
734 ;-------------------------------------------------------------------------------
735 (defun transform-a-sequence!
736 (a-function
737 the-seq
738 &key
739 (result the-seq))
740 "given
741 [1] a-function (required)
742 ==> a function with one argument
743 [2] the-seq (required)
744 ==> a sequence of numbers
745 [3] result (keyword; the-seq)
746 <== a sequence of numbers
747 returns
748 [1] a sequence of numbers obtained by applying
749 a-function to each element of the-seq"
750 (map-into result a-function the-seq))
751
752 #|
753 (setf *a* #(1 10 100))
754 (transform-a-sequence! #'convert-to-dB *a*) ;==> #(0.0 10.0 20.0)
755 *a* ;==> #(0.0 10.0 20.0)
756 |#
757
758 ;-------------------------------------------------------------------------------
759 (defun linear-interpolation!
760 (a-sequence
761 &key
762 (interpolation-predicate #'zerop))
763 "given
764 [1] a-sequence (required)
765 <=> a sequence of numbers
766 [2] interpolation-predicate (keyword; #'zerop)
767 ==> a unitary predicate which, when true,
768 indicates that a value in a-sequence
769 is to be replaced by a linear interpolation
770 returns
771 [1] a-sequence, modified with linear interpolates
772 over specified values"
773 (let ((n (length a-sequence)))
774 ;;; check to make sure that the beginning and the end of the sequence
775 ;;; are well-defined and that the sequence is long enough to do something
776 ;;; with
777 (assert (and
778 (> n 2)
779 (not (funcall interpolation-predicate (elt a-sequence 0)))
780 (not (funcall interpolation-predicate (elt a-sequence (1- n))))))
781 (do ((i-start 0)
782 (i-stop (1- n))
783 i-first-bad i-last-bad
784 i-good-before-first-bad
785 i-good-after-last-bad
786 rise run)
787 ((>= i-start i-stop) a-sequence) ;;; exit test
788 (block this-block
789 (dotimes (i (1+ (- i-stop i-start)) (setf i-start (1+ i-stop)))
790 (when (funcall interpolation-predicate
791 (elt a-sequence (+ i i-start)))
792 (setf i-first-bad (+ i i-start)
793 i-good-before-first-bad (1- i-first-bad))
794 (dotimes (j (- i-stop i-first-bad)
795 (error "algorithm failure!!! --- evacuate building!!!"))
796 (when (not (funcall interpolation-predicate
797 (elt a-sequence (+ j i-first-bad 1))))
798 (setf i-last-bad (+ j i-first-bad)
799 i-good-after-last-bad (1+ i-last-bad)
800 rise (- (elt a-sequence i-good-after-last-bad)
801 (elt a-sequence i-good-before-first-bad))
802 run (- i-good-after-last-bad i-good-before-first-bad))
803 (dotimes (k (- i-good-after-last-bad i-first-bad))
804 (setf (elt a-sequence (+ i-first-bad k))
805 (+ (elt a-sequence i-good-before-first-bad)
806 (* (/ (1+ k) run) rise))))
807 (setf i-start (1+ i-good-after-last-bad))
808 (return-from this-block)))))))))
809
810 #|
811 (setf *a* #(1 10 0 100))
812 (linear-interpolation! *a*) ;==> #(1 10 55 100)
813 *a* ;==> #(1 10 55 100)
814 (setf *b* #(1 10 nil 100))
815 (linear-interpolation!
816 *b*
817 :interpolation-predicate #'null) ;==> #(1 10 55 100)
818 *b* ;==> #(1 10 55 100)
819 |#
820
821 ;-------------------------------------------------------------------------------
822 ;-------------------------------------------------------------------------------
823 ;;; The functions dot-product
824 ;;; euclidean-norm
825 ;;; sum-of-squares
826 ;;; sum
827 ;;; min-and-max-of-seq
828 ;;; min-of-seq
829 ;;; max-of-seq
830 ;;; binary-search
831 ;;; compare-seqs
832 ;;; all take sequence(s) of numbers as input and return various scalar(s)
833 ;;; of interest. For example, the dot-product between two sequences
834 ;;; of the same length.
835 ;-------------------------------------------------------------------------------
836 ;-------------------------------------------------------------------------------
837 (defun dot-product (x y)
838 "given
839 [1] x (required)
840 ==> a sequence of numbers
841 [2] y (required)
842 ==> a sequence of numbers
843 returns
844 [1] the dot product of x and y
845 ----
846 Note: corresponds to SDOT in LINPACK's bla"
847 (let ((sum 0.0))
848 (dotimes (i (length x) sum)
849 (incf sum (* (elt x i) (elt y i))))))
850
851 ;;; (dot-product #(1 2 3) #(4 5 6)) ;==> 32.0
852
853 ;-------------------------------------------------------------------------------
854 (defun euclidean-norm (x)
855 "given
856 [1] x (required)
857 ==> a sequence of numbers
858 returns
859 [1] Euclidean norm of x
860 ----
861 Note: corresponds to LINPACK's SNRM2,
862 but this is POOR implementation!"
863 (sqrt (dot-product x x)))
864
865 ;;; (euclidean-norm #(4 3)) ;==> 5.0
866
867 ;-------------------------------------------------------------------------------
868 (defun sum-of-squares
869 (the-seq
870 &key
871 (start 0)
872 (end (length the-seq)))
873 "given
874 [1] the-seq (required)
875 ==> a sequence of numbers
876 [2] start (keyword; 0)
877 ==> start index of sequence to be used
878 [3] end (keyword; length of the-seq)
879 ==> 1 + end index of sequence to be used
880 returns
881 [1] sum of squares of elements
882 in specified subsequence"
883 #-allegro
884 (reduce #'+ the-seq :start start :end end :key #'(lambda (x) (* x x)))
885 #+allegro
886 (let ((SS 0.0)
887 (j start))
888 (dotimes (i (- end start) SS)
889 (incf SS (expt (elt the-seq j) 2))
890 (incf j)
891 ))
892 )
893
894 #|
895 (sum-of-squares #(4 3)) ;==> 25
896 (sum-of-squares #(1.0 2.0 3.0)) ;==> 14.0
897 (sum-of-squares #(1.0 2.0 3.0 4.0)) ;==> 30.0
898 (sum-of-squares #(1.0 2.0 3.0 4.0) :end 3) ;==> 14.0
899 (sum-of-squares #(1.0 2.0 3.0 4.0) :start 1 :end 3) ;==> 13.0
900 |#
901
902 ;-------------------------------------------------------------------------------
903 (defun sum
904 (the-seq
905 &key
906 (start 0)
907 (end (length the-seq)))
908 "given
909 [1] the-seq (required)
910 ==> a sequence of numbers
911 [2] start (keyword; 0)
912 ==> start index of sequence to be used
913 [3] end (keyword; length of the-seq)
914 ==> 1 + end index of sequence to be used
915 returns
916 [1] sum of elements in specified subsequence"
917 (reduce #'+ the-seq :start start :end end))
918
919 ;;; (sum #(4 3)) ;==> 7
920 ;;; (sum #(7 8 3 4 3 10 11)) ;==> 46
921 ;;; (sum #(7 8 3 4 3 10 11) :start 3 :end 5) ;==> 7
922
923 ;-------------------------------------------------------------------------------
924 (defun min-and-max-of-seq
925 (the-seq
926 &key
927 (start 0)
928 (end (length the-seq)))
929 "given
930 [1] the-seq (required)
931 ==> a sequence of real-valued numbers
932 [2] start (keyword; 0)
933 ==> start index of sequence to be used
934 [3] end (keyword; length of the-seq)
935 ==> 1 + end index of sequence to be used
936 returns
937 [1] minimum value in sequence
938 [2] maximum value in sequence"
939 (values (reduce #'min the-seq :start start :end end)
940 (reduce #'max the-seq :start start :end end)))
941
942 ;;; (min-and-max-of-seq #(7 8 3 4 3 10 11)) ;==> 3 and 11
943
944 ;-------------------------------------------------------------------------------
945 (defun min-of-seq
946 (the-seq
947 &key
948 (start 0)
949 (end (length the-seq)))
950 "given
951 [1] the-seq (required)
952 ==> a sequence of real-valued numbers
953 [2] start (keyword; 0)
954 ==> start index of sequence to be used
955 [3] end (keyword; length of the-seq)
956 ==> 1 + end index of sequence to be used
957 returns
958 [1] minimum value in sequence
959 [2] index of minimum value"
960 (let ((the-min (reduce #'min the-seq :start start :end end)))
961 (values the-min (position the-min the-seq :start start :end end))))
962
963 ;;; (min-of-seq #(7 8 3 4 3 10 11)) ;==> 3 and 2
964
965 ;-------------------------------------------------------------------------------
966 (defun max-of-seq
967 (the-seq
968 &key
969 (start 0)
970 (end (length the-seq)))
971 "given
972 [1] the-seq (required)
973 ==> a sequence of real-valued numbers
974 [2] start (keyword; 0)
975 ==> start index of sequence to be used
976 [3] end (keyword; length of the-seq)
977 ==> 1 + end index of sequence to be used
978 returns
979 [1] maximum value in sequence
980 [2] index of maximum value"
981 (let ((the-max (reduce #'max the-seq :start start :end end)))
982 (values the-max (position the-max the-seq :start start :end end))))
983
984 ;;; (max-of-seq #(7 8 3 4 3 10 11)) ;==> 11 and 6
985
986 ;-------------------------------------------------------------------------------
987 (defun binary-search
988 (value
989 ordered-seq
990 &key
991 (lower-test #'<=)
992 (upper-test #'<=))
993 "given
994 [1] value (required)
995 ==> a real-valued number
996 [2] ordered-seq (required)
997 ==> an ordered sequence of real-valued numbers
998 [3] lower-test (keyword; #'<=)
999 ==> a predicate of two arguments
1000 that defines the lower test
1001 [4] upper-test (keyword; #'<=)
1002 ==> a predicate of two arguments
1003 that defines the upper test
1004 returns
1005 [1] if lower-test and upper-test are set
1006 to their default values, this function
1007 returns an index i such that v[i] <= value <= v[i+1],
1008 where 0 <= i <= n-2 (if there is no such i, nil is returned);
1009 if nondefault values are supplied for lower-test and upper-test,
1010 then the condition `v[i] <= value <= v[i+1]'
1011 gets modified in an obvious way
1012 ---
1013 Note: value can be an arbitrary object, and ordered-seq
1014 can be a list of arbitrary objects if the binary predicates
1015 lower-test and upper-test are appropriate set"
1016 (binary-search-internal
1017 value ordered-seq 0 (1- (length ordered-seq)) lower-test upper-test))
1018
1019 ;;; (binary-search 8 #(1 2 4 5 7 9 11)) ;==> 4
1020
1021 ;-------------------------------------------------------------------------------
1022 (defun compare-seqs (one-seq another-seq)
1023 "given
1024 [1] one-seq (required)
1025 ==> any sequence
1026 [2] another-seq (required)
1027 ==> any sequence (must be same length as one-seq)
1028 returns
1029 [1] maximum absolute difference between corresponding
1030 elements of the two sequences
1031 [2] average absolute difference
1032 ---
1033 Note: useful for comparing results that in theory
1034 should be identical"
1035 (assert (= (length one-seq) (length another-seq)))
1036 (let ((max-abs-diff 0.0)
1037 (ave-abs-diff 0.0)
1038 (n (length one-seq))
1039 abs-diff)
1040 (dotimes (i n (values (/ ave-abs-diff n) max-abs-diff))
1041 (setf abs-diff (abs (- (elt one-seq i) (elt another-seq i))))
1042 (when (> abs-diff max-abs-diff) (setf max-abs-diff abs-diff))
1043 (incf ave-abs-diff abs-diff))))
1044
1045 #|
1046 (compare-seqs
1047 #(1 2 3 4 5 6 7 8 9 10)
1048 #(1 2 3 4 5 6.1 7 8 9 10))
1049 ;==>
1050 0.009999999999999964
1051 0.09999999999999964
1052 |#
1053
1054 ;-------------------------------------------------------------------------------
1055 ;-------------------------------------------------------------------------------
1056 ;;; The functions iota
1057 ;;; sample-from-a-function
1058 ;;; both generate a sequence of numbers. In the case of iota, the sequence
1059 ;;; of numbers is equally spaced; in the case of sample-from-a-function,
1060 ;;; the sequence is obtained by sampling a function at a specified set
1061 ;;; of points
1062 ;-------------------------------------------------------------------------------
1063 ;-------------------------------------------------------------------------------
1064 (defun iota
1065 (start
1066 end
1067 &key
1068 (type 'vector)
1069 (float-p nil)
1070 (skip 1)
1071 (result (make-sequence type (1+ (/ (- end start) skip)))))
1072 "given
1073 [1] start (required)
1074 ==> an integer
1075 [2] end (required)
1076 ==> an integer >= start
1077 [3] type (keyword; 'vector)
1078 ==> type of sequence desired
1079 [4] float-p (keyword; nil)
1080 ==> if true, sequence consists of floats
1081 [5] skip (keyword; 1)
1082 ==> desired increment between integers in sequence
1083 [6] result (keyword; sequence of length (1+ (- start end)))
1084 ==> storage space for result
1085 returns
1086 [1] a sequence of length (1+ (/ (- end start) skip))
1087 with numbers start, start + skip, ..., end - skip, end
1088 (these are stored in whatever result points to)"
1089 (if float-p
1090 (dotimes (i (1+ (/ (- end start) skip)) result)
1091 (setf (elt result i) (float (+ start (* i skip)))))
1092 (dotimes (i (1+ (/ (- end start) skip)) result)
1093 (setf (elt result i) (+ start (* i skip))))))
1094
1095 ;;; (iota 7 10 :type 'vector) ;==> #(7 8 9 10)
1096
1097 ;-------------------------------------------------------------------------------
1098 (defun sample-from-a-function
1099 (a-function
1100 &key
1101 (seq-of-x-values '())
1102 (x0 0.0)
1103 (delta-x 1.0)
1104 (n (length seq-of-x-values)))
1105 "given
1106 [1] a-function (required)
1107 ==> a function to be sampled from
1108 [2] seq-of-x-values (keyword; '())
1109 ==> sequence of values at which a-function
1110 is to be sampled (default '())
1111 [3] x0 (keyword; 0.0)
1112 ==> point of first value to be sampled
1113 [4] delta-x (keyword; 1.0)
1114 ==> increment between points
1115 [5] n (keyword; length of seq-of-x-values)
1116 ==> number of samples
1117 returns
1118 [1] an array with the values of a-function
1119 at a specified set of points; and
1120 [2] an array with the specified set of points.
1121 ---
1122 Note that the specified set of points either is in seq-of-x-values
1123 -- if it is non-nil -- or is given by
1124 x0, x0 + delta-x, x0 + 2*delta-x, ..., x0 + (n-1)*delta-x"
1125 (let ((the-samples (make-array n))
1126 (x (if seq-of-x-values
1127 (if (arrayp seq-of-x-values)
1128 seq-of-x-values
1129 (make-array n :initial-contents seq-of-x-values))
1130 (make-array n))))
1131 (cond
1132 (seq-of-x-values
1133 (dotimes (i n (values the-samples x))
1134 (setf (aref the-samples i)
1135 (funcall a-function (elt seq-of-x-values i)))))
1136 (t
1137 (dotimes (i n (values the-samples x))
1138 (setf (aref x i) x0)
1139 (setf (aref the-samples i) (funcall a-function x0))
1140 (incf x0 delta-x))))))
1141
1142 #|
1143 (sample-from-a-function #'cos
1144 :delta-x (/ (* 2 pi) 4)
1145 :n 5)
1146 ;==>
1147 #(1.0 6.123031769111886E-17 -1.0 -1.836909530733566E-16 1.0)
1148 #(0.0 1.5707963267948966 3.141592653589793 4.71238898038469 6.283185307179586)
1149 |#
1150
1151 ;-------------------------------------------------------------------------------
1152 ;-------------------------------------------------------------------------------
1153 ;;; The functions sign
1154 ;;; sampling-time->Nyquist-frequency
1155 ;;; Nyquist-frequency->sampling-time
1156 ;;; are, respectively, a Lisp implementation of the Fortran SIGN function
1157 ;;; (useful because it differs from the Common Lisp function signum);
1158 ;;; a simple function that takes the sampling time (called delta t
1159 ;;; in the SAPA book) and returns the Nyquist frequency; and another
1160 ;;; simple function for going from the Nyquist frequency to the
1161 ;;; sampling time
1162 ;-------------------------------------------------------------------------------
1163 ;-------------------------------------------------------------------------------
1164 (defun sign (a1 a2)
1165 "implementation of Fortran SIGN function"
1166 (if (minusp a2)
1167 (- (abs a1))
1168 (abs a1)))
1169
1170 ;-------------------------------------------------------------------------------
1171 (defun sampling-time->Nyquist-frequency (sampling-time)
1172 "given a sampling time, returns the associated Nyquist frequency
1173 ---
1174 Note: see page 98 of the SAPA book"
1175 (/ (* 2 sampling-time)))
1176
1177 ;-------------------------------------------------------------------------------
1178 (defun Nyquist-frequency->sampling-time (Nyquist-frequency)
1179 "given a Nyquist frequency, returns the associated sampling-time
1180 ---
1181 Note: see page 98 of the SAPA book"
1182 (/ (* 2 Nyquist-frequency)))
1183
1184 ;-------------------------------------------------------------------------------
1185 ;-------------------------------------------------------------------------------
1186 ;;; Everything below here consists of internal symbols in the SAPA package
1187 ;;; and should be regarded as "dirty laundry" ...
1188 ;-------------------------------------------------------------------------------
1189 ;-------------------------------------------------------------------------------
1190 ;;; Used by power-of-2
1191 (defun power-of-2+1 (n)
1192 (cond
1193 ((or (not (numberp n)) (not (plusp n)) (< n 1))
1194 nil)
1195 ((= n 1) 1)
1196 (t
1197 (let ((local-answer (power-of-2+1 (/ n 2))))
1198 (if (numberp local-answer)
1199 (1+ local-answer))))))
1200
1201 ;-------------------------------------------------------------------------------
1202 ;;; NOTE: used by binary-search. This arrangement was done because the way
1203 ;;; that this function is recursively called makes the lambda list
1204 ;;; bulky as a user interface. Originally, we had &optional instead of
1205 ;;; &key in binary-search, and there was no binary-search-internal.
1206 ;;; There is probably a slick way of getting around this is Lisp,
1207 ;;; but I (dbp) don't know what it is as of 7/29/91!
1208 ;;; NOTE: `end' here is one less than the usual Lisp parameter of that name
1209 (defun binary-search-internal
1210 (value
1211 ordered-seq
1212 start
1213 end
1214 lower-test
1215 upper-test)
1216 (cond
1217 ((= (- end start) 1)
1218 (if (and (funcall lower-test (elt ordered-seq start) value)
1219 (funcall upper-test value (elt ordered-seq end)))
1220 start)) ;returns nil if termination condition fails
1221 (t
1222 (let ((middle (truncate (+ start end) 2)))
1223 (if (funcall upper-test value (elt ordered-seq middle))
1224 (binary-search-internal
1225 value ordered-seq start middle lower-test upper-test)
1226 (binary-search-internal
1227 value ordered-seq middle end lower-test upper-test)
1228 )))))

  ViewVC Help
Powered by ViewVC 1.1.5