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

Contents of /sapaclisp/examples.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (vendor branch)
Mon May 9 20:24:42 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: :CL-USER; 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 ; examples.lisp
12 ;
13 ; examples showing how to repeat certain results in the SAPA book
14 ; using the functions in the SAPA package; to run these examples,
15 ; you should first evaluate the "in-package" and "use-package" forms
16 ; and then all the "defvar" forms, after which you can evaluate the
17 ; Lisp forms related to the various figures in the SAPA book
18 ; in any order you choose.
19 ;
20 ;-----------------------------------------------------------------------------
21 (in-package :CL-USER)
22
23 (use-package :SAPA)
24
25 ;-------------------------------------------------------------------------------
26 ;-------------------------------------------------------------------------------
27 ;;; We need to define the following special variables in order to run
28 ;;; the test cases in this file.
29 ;-------------------------------------------------------------------------------
30 ;-------------------------------------------------------------------------------
31 (defvar *rotation-of-earth-data*
32 (vector 71.0 63.0 70.0 88.0 99.0 90.0 110.0 135.0 128.0 154.0 156.0
33 141.0 131.0 132.0 141.0 104.0 136.0 146.0 124.0 129.0 110.0 119.0
34 97.0 115.0 119.0 113.0 115.0 127.0 122.0 147.0 123.0 130.0 149.0
35 183.0 186.0 185.0 189.0 203.0 191.0 217.0 220.0 245.0 213.0 227.0
36 247.0 266.0 221.0 253.0 231.0 236.0 243.0 246.0 238.0 261.0 248.0
37 264.0 270.0 282.0 255.0 286.0 277.0 260.0 261.0 260.0 289.0 297.0
38 319.0 311.0 314.0 309.0 309.0 310.0 315.0 310.0 296.0 264.0 284.0
39 260.0 286.0 271.0 271.0 259.0 279.0 268.0 296.0 280.0 310.0 265.0
40 277.0 257.0 296.0 307.0 306.0 266.0 285.0 279.0 271.0 257.0 270.0
41 232.0))
42
43 (defvar *centered-rotation-of-earth-data*
44 (x+b *rotation-of-earth-data* (- (sample-mean *rotation-of-earth-data*))))
45
46 (defvar *the-acvs*)
47 (defvar *C_h*)
48
49 (defvar *eigenvalues-NW-4*)
50 (defvar *eigenspectra-NW-4*)
51 (defvar *freqs*)
52
53 (defvar *ar-5-coeffs*)
54 (defvar *forward-pred-errors*)
55
56 ;-------------------------------------------------------------------------------
57 ;-------------------------------------------------------------------------------
58 ;;; Some of the examples below recreate analysis on time series as described
59 ;;; in the SAPA book. These time series can be obtained from StatLib
60 ;;; by sending the message
61 ;;; send sapa from datasets
62 ;;; (for details, see page xxvii of the SAPA book). Once these time
63 ;;; series have been obtained from StatLib, it is necessary to strip
64 ;;; out the pertinent numbers for each times series, place them into an
65 ;;; appropriate ASCII file, and then read each file to load the appropriate
66 ;;; time series values into Lisp. The actual details of reading these files
67 ;;; depend on which implementation of Common Lisp is used. Here we show
68 ;;; how it is done in Macintosh Common Lisp 2.0 (MCL 2.0) and in Symbolics
69 ;;; Genera 8.1.1 running on a Symbolics MacIvory model 3. Evaluation of each
70 ;;; of the defvar forms below loads in a time series. For example, in MCL 2.0,
71 ;;; "ccl:SAPA;ocean-wave" refers to a file named ocean-wave in a folder
72 ;;; called SAPA that lives in the folder containing the MCL 2.0 application.
73 ;;; The file "ccl:SAPA;ocean-wave" consists of the 1024 lines from StatLib
74 ;;; containing the ocean wave time series (one number each line).
75 ;-------------------------------------------------------------------------------
76 ;-------------------------------------------------------------------------------
77 (defvar *ocean-wave-data*
78 (let ((the-array (make-array 1024)))
79 (with-open-file (in-stream
80 #+mcl "ccl:SAPA;ocean-wave"
81 #+genera "L:>dbp>ocean-wave..newest"
82 #+allegro "ocean-wave"
83 :direction :input)
84 (dotimes (index (length the-array) the-array)
85 (setf (svref the-array index)
86 (float (read in-stream)))))))
87
88 ;-------------------------------------------------------------------------------
89 (defvar *rough-ice-profile-data*
90 (let ((the-array (make-array 1121)))
91 (with-open-file (in-stream
92 #+mcl "ccl:SAPA;rough-ice-profile"
93 #+genera "L:>dbp>rough-ice-profile..newest"
94 #+allegro "rough-ice-profile"
95 :direction :input)
96 (dotimes (index (length the-array) the-array)
97 (setf (svref the-array index)
98 (float (read in-stream)))))))
99
100 ;-------------------------------------------------------------------------------
101 (defvar *smooth-ice-profile-data*
102 (let ((the-array (make-array 288)))
103 (with-open-file (in-stream
104 #+mcl "ccl:SAPA;smooth-ice-profile"
105 #+genera "L:>dbp>smooth-ice-profile..newest"
106 #+allegro "smooth-ice-profile"
107 :direction :input)
108 (dotimes (index (length the-array) the-array)
109 (setf (svref the-array index)
110 (float (read in-stream)))))))
111
112 ;-------------------------------------------------------------------------------
113 (defvar *Willamette-River-data*
114 (let ((the-array (make-array 395)))
115 (with-open-file (in-stream
116 #+mcl "ccl:SAPA;Willamette-River"
117 #+genera "L:>dbp>Willamette-River..newest"
118 #+allegro "Willamette-River"
119 :direction :input)
120 (dotimes (index (length the-array) the-array)
121 (setf (svref the-array index)
122 (float (read in-stream)))))))
123
124 ;-------------------------------------------------------------------------------
125 (defvar *AR-2*
126 (let ((the-array (make-array 1024)))
127 (with-open-file (in-stream
128 #+mcl "ccl:SAPA;AR-2"
129 #+genera "L:>dbp>AR-2..newest"
130 #+allegro "AR-2"
131 :direction :input)
132 (dotimes (index (length the-array) the-array)
133 (setf (svref the-array index)
134 (float (read in-stream)))))))
135
136 #|
137 (svref *ocean-wave-data* 0) ;==> 477.0
138 (svref *ocean-wave-data* 1023) ;==> -113.0
139
140 (svref *rough-ice-profile-data* 0) ;==> -24.6
141 (svref *rough-ice-profile-data* 1120) ;==> -22.5
142
143 (svref *smooth-ice-profile-data* 0) ;==> -9.1
144 (svref *smooth-ice-profile-data* 287) ;==> -8.9
145
146 (svref *Willamette-River-data* 0) ;==> 8.95361
147 (svref *Willamette-River-data* 394) ;==> 9.06933
148
149 (svref *AR-2* 0) ;==> 1.619873110842541
150 (svref *AR-2* 1023) ;==> 1.081772237497269
151 |#
152
153 ;-------------------------------------------------------------------------------
154 ;-------------------------------------------------------------------------------
155 ;;; Figure 172: filtering of rotation of earth data using simple
156 ;;; 3 point filter with coefficients 1/4, 1/2 and 1/4.
157 ;-------------------------------------------------------------------------------
158 ;-------------------------------------------------------------------------------
159 (let* ((filtered-series (filter-time-series
160 *centered-rotation-of-earth-data*
161 #(1/4 1/2 1/4)))
162 (residuals (x-y (make-array
163 (length filtered-series)
164 :displaced-to *centered-rotation-of-earth-data*
165 :displaced-index-offset 1)
166 filtered-series)))
167 (format t "~&Figure 172 ...")
168 (dotimes (i 5)
169 (format t "~&~8,2F ~8,2F ~8,2F"
170 (svref *centered-rotation-of-earth-data* (1+ i))
171 (svref filtered-series i)
172 (svref residuals i)))
173 (format t "~&...")
174 (let ((N-fs (length filtered-series)))
175 (dotimes (i 5)
176 (format t "~&~8,2F ~8,2F ~8,2F"
177 (svref *centered-rotation-of-earth-data* (1+ (+ i (- N-fs 5))))
178 (svref filtered-series (+ i (- N-fs 5)))
179 (svref residuals (+ i (- N-fs 5))))))
180 (format t "~&... Figure 172")
181 (values))
182
183 #|
184 ;;; Here is what is printed when the above form is evaluated:
185 Figure 172 ...
186 -152.73 -148.98 -3.75
187 -145.73 -142.98 -2.75
188 -127.73 -129.48 1.75
189 -116.73 -121.73 5.00
190 -125.73 -118.48 -7.25
191 ...
192 69.27 63.02 6.25
193 63.27 62.77 0.50
194 55.27 53.77 1.50
195 41.27 48.02 -6.75
196 54.27 41.52 12.75
197 ... Figure 172
198 |#
199
200 ;-------------------------------------------------------------------------------
201 ;-------------------------------------------------------------------------------
202 ;;; Figure 177: squared modulus of transfer function for
203 ;; the Kth order least squares approximation
204 ;;; to an ideal low-pass filter
205 ;-------------------------------------------------------------------------------
206 ;-------------------------------------------------------------------------------
207 (let ((65-point-ls-filter
208 (create-least-squares-low-pass-filter 65 0.1)))
209 (format t "~&Figure 177 ...")
210 (dolist (i '(0 1))
211 (format t "~&~3D: ~F" (- i 32) (svref 65-point-ls-filter i)))
212 (format t "~&...")
213 (dolist (i '(31 32 33))
214 (format t "~&~3D: ~F" (- i 32) (svref 65-point-ls-filter i)))
215 (format t "~&...")
216 (dolist (i '(63 64))
217 (format t "~&~3D: ~F" (- i 32) (svref 65-point-ls-filter i)))
218 (format t "~& sum = ~F" (sum 65-point-ls-filter))
219 (multiple-value-bind (mod-sq-trans-func freqs)
220 (transfer-function-for-filter
221 65-point-ls-filter
222 :return-frequencies-p t)
223 (dolist (i '(0 1 2 3 4 5))
224 (format t "~&~6,4F: ~F"
225 (svref freqs i)
226 (svref mod-sq-trans-func i)))
227 (format t "~&...")
228 (dolist (i '(51 52 53))
229 (format t "~&~6,4F: ~F"
230 (svref freqs i)
231 (svref mod-sq-trans-func i)))
232 (format t "~&...")
233 (dolist (i '(254 255 256))
234 (format t "~&~6,4F: ~F"
235 (svref freqs i)
236 (svref mod-sq-trans-func i))))
237 (format t "~&... Figure 177")
238 (values))
239
240 #|
241 ;;; Here is what is printed when the above form is evaluated:
242 Figure 177 ...
243 -32: 0.009474353340135194
244 -31: 0.00604435859161769
245 ...
246 -1: 0.18737511634014858
247 0: 0.2002963792180472
248 1: 0.18737511634014858
249 ...
250 31: 0.00604435859161769
251 32: 0.009474353340135194
252 sum = 1.0000000000000002
253 0.0000: 1.9286549331065735E-15
254 0.0020: -9.97589282841272E-4
255 0.0039: -0.003518706679958886
256 0.0059: -0.006268132747772469
257 0.0078: -0.0074626548387638734
258 0.0098: -0.005314313156846105
259 ...
260 0.0996: -5.4415656590533095
261 0.1016: -7.797097385966063
262 0.1035: -10.914299436697938
263 ...
264 0.4961: -42.85386826084243
265 0.4980: -40.44459647657306
266 0.5000: -39.73441540905796
267 ... Figure 177
268 |#
269
270 ;-------------------------------------------------------------------------------
271 ;-------------------------------------------------------------------------------
272 ;;; Figure 178: squared modulus of transfer function for
273 ;;; the Kth order least squares
274 ;;; + triangular convergence factors approximation
275 ;;; to an ideal low-pass filter
276 ;-------------------------------------------------------------------------------
277 ;-------------------------------------------------------------------------------
278 (let ((65-point-ls-filter
279 (create-least-squares-low-pass-filter
280 65 0.1
281 :convergence-factors
282 #'(lambda (k) (triangular-convergence-factors k 65)))))
283 (format t "~&Figure 178 ...")
284 (dolist (i '(0 1))
285 (format t "~&~3D: ~F" (- i 32) (svref 65-point-ls-filter i)))
286 (format t "~&...")
287 (dolist (i '(31 32 33))
288 (format t "~&~3D: ~F" (- i 32) (svref 65-point-ls-filter i)))
289 (format t "~&...")
290 (dolist (i '(63 64))
291 (format t "~&~3D: ~F" (- i 32) (svref 65-point-ls-filter i)))
292 (format t "~& sum = ~F" (sum 65-point-ls-filter))
293 (multiple-value-bind (mod-sq-trans-func freqs)
294 (transfer-function-for-filter
295 65-point-ls-filter
296 :return-frequencies-p t)
297 (dolist (i '(0 1 2 3 4 5))
298 (format t "~&~6,4F: ~F"
299 (svref freqs i)
300 (svref mod-sq-trans-func i)))
301 (format t "~&...")
302 (dolist (i '(51 52 53))
303 (format t "~&~6,4F: ~F"
304 (svref freqs i)
305 (svref mod-sq-trans-func i)))
306 (format t "~&...")
307 (dolist (i '(254 255 256))
308 (format t "~&~6,4F: ~F"
309 (svref freqs i)
310 (svref mod-sq-trans-func i)))
311 )
312 (format t "~&... Figure 178")
313 (values))
314
315 #|
316 ;;; Here is what is printed when the above form is evaluated:
317 Figure 178 ...
318 -32: 2.958988584176635E-4
319 -31: 3.7754952616136017E-4
320 ...
321 -1: 0.18726456497603497
322 0: 0.20643377319025294
323 1: 0.18726456497603497
324 ...
325 31: 3.7754952616136017E-4
326 32: 2.958988584176635E-4
327 sum = 0.9999999999999998
328 0.0000: -1.9286549331065743E-15
329 0.0020: 9.19973377882888E-4
330 0.0039: 0.0035297560721751367
331 0.0059: 0.007399570616709966
332 0.0078: 0.011879123556194419
333 0.0098: 0.016187305184551266
334 ...
335 0.0996: -5.636216008278292
336 0.1016: -6.812966030344137
337 0.1035: -8.143286403474598
338 ...
339 0.4961: -49.51447433638478
340 0.4980: -49.421175373585505
341 0.5000: -49.38851214695436
342 ... Figure 178
343 |#
344
345 ;-------------------------------------------------------------------------------
346 ;-------------------------------------------------------------------------------
347 ;;; Figure 181: squared modulus of transfer function for
348 ;;; dpss approximations to an ideal low-pass filter
349 ;-------------------------------------------------------------------------------
350 ;-------------------------------------------------------------------------------
351 (let ((17-point-dpss
352 (dpss-data-taper!
353 (make-array 17 :initial-element 1.0)
354 :taper-parameter (* 0.1 17)))
355 (65-point-dpss
356 (dpss-data-taper!
357 (make-array 65 :initial-element 1.0)
358 :taper-parameter (* 0.1 65))))
359 (a*x! (/ (sum 17-point-dpss)) 17-point-dpss)
360 (a*x! (/ (sum 65-point-dpss)) 65-point-dpss)
361 (format t "~&Figure 181 ...")
362 (multiple-value-bind (mod-sq-trans-func-17 freqs)
363 (transfer-function-for-filter
364 17-point-dpss
365 :return-frequencies-p t
366 :N-fft 512)
367 (let ((mod-sq-trans-func-65
368 (transfer-function-for-filter
369 65-point-dpss
370 :return-frequencies-p t
371 :N-fft 512)))
372 (dolist (i '(0 1 2 3 4 5 6 7))
373 (format t "~&~6,4F: ~7,3F, ~7,3F"
374 (svref freqs i)
375 (svref mod-sq-trans-func-17 i)
376 (svref mod-sq-trans-func-65 i)))
377 (format t "~&...")
378 (dolist (i '(50 51 52 53 54))
379 (format t "~&~6,4F: ~7,3F, ~7,3F"
380 (svref freqs i)
381 (svref mod-sq-trans-func-17 i)
382 (svref mod-sq-trans-func-65 i)))
383 (format t "~&...")
384 (dolist (i '(254 255 256))
385 (format t "~&~6,4F: ~7,3F, ~7,3F"
386 (svref freqs i)
387 (svref mod-sq-trans-func-17 i)
388 (svref mod-sq-trans-func-65 i)))
389 ))
390 (format t "~&... Figure 181")
391 (values))
392
393 #|
394 ;;; Here is what is printed when the above form is evaluated:
395 Figure 181 ...
396 0.0000: 0.000, -0.0000
397 0.0020: -0.008, -0.033
398 0.0039: -0.030, -0.131
399 0.0059: -0.068, -0.294
400 0.0078: -0.121, -0.523
401 0.0098: -0.189, -0.818
402 0.0117: -0.272, -1.179
403 0.0137: -0.371, -1.607
404 ...
405 0.0977: -26.612, -134.511
406 0.0996: -28.537, -176.340
407 0.1016: -30.737, -157.182
408 0.1035: -33.328, -172.494
409 0.1055: -36.530, -160.819
410 ...
411 0.4961: -49.526, -184.556
412 0.4980: -49.369, -182.150
413 0.5000: -49.318, -181.440
414 ... Figure 181
415 ;;; Note: the above listing is from Macintosh Common Lisp. When this form
416 ;;; is evaluated under Genera 8.1 or Allegro Common Lisp, the last
417 ;;; column differs somewhat at low dB levels:
418 0.0977: -26.612, -135.267
419 0.0996: -28.537, -157.133
420 0.1016: -30.737, -150.720
421 0.1035: -33.328, -157.869
422 0.1055: -36.530, -164.541
423 ...
424 0.4961: -49.526, -163.929
425 0.4980: -49.369, -165.401
426 0.5000: -49.318, -165.964
427 ;;; This difference is evidently due to different default precisions
428 ;;; for numerical computations.
429 |#
430
431 ;-------------------------------------------------------------------------------
432 ;-------------------------------------------------------------------------------
433 ;;; Figure 182: squared modulus of transfer function for
434 ;;; approximations to an ideal low-pass filter
435 ;;; using dpss as convergence factors to least squares
436 ;;; approximations
437 ;-------------------------------------------------------------------------------
438 ;-------------------------------------------------------------------------------
439 (let ((delta-point-04
440 (create-dpss-low-pass-filter
441 65 0.04 0.1))
442 (delta-point-01
443 (create-dpss-low-pass-filter
444 65 0.01 0.1)))
445 (a*x! (/ (sum delta-point-04)) delta-point-04)
446 (a*x! (/ (sum delta-point-01)) delta-point-01)
447 (format t "~&Figure 182 ...")
448 (print (sum delta-point-04))
449 (print (sum delta-point-01))
450 (multiple-value-bind (mod-sq-trans-func-04 freqs)
451 (transfer-function-for-filter
452 delta-point-04
453 :return-frequencies-p t
454 :N-fft 512)
455 (let ((mod-sq-trans-func-01
456 (transfer-function-for-filter
457 delta-point-01
458 :return-frequencies-p t
459 :N-fft 512)))
460 (dolist (i '(0 1 2 3 4 5 6))
461 (format t "~&~6,4F: ~7,3F, ~7,3F"
462 (svref freqs i)
463 (svref mod-sq-trans-func-04 i)
464 (svref mod-sq-trans-func-01 i)))
465 (format t "~&...")
466 (dolist (i '(50 51 52 53 54))
467 (format t "~&~6,4F: ~7,3F, ~7,3F"
468 (svref freqs i)
469 (svref mod-sq-trans-func-04 i)
470 (svref mod-sq-trans-func-01 i)))
471 (format t "~&...")
472 (dolist (i '(254 255 256))
473 (format t "~&~6,4F: ~7,3F, ~7,3F"
474 (svref freqs i)
475 (svref mod-sq-trans-func-04 i)
476 (svref mod-sq-trans-func-01 i)))))
477 (format t "~&... Figure 182")
478 (values))
479
480 #|
481 ;;; Here is what is printed when the above form is evaluated:
482 Figure 182 ...
483 0.9999999999999999
484 1.0
485 0.0000: 0.000, 0.000
486 0.0020: 0.0001, 0.0005
487 0.0039: 0.0002, 0.002
488 0.0059: 0.0004, 0.005
489 0.0078: 0.001, 0.009
490 0.0098: 0.001, 0.014
491 0.0117: 0.001, 0.021
492 ...
493 0.0977: -4.944, -4.095
494 0.0996: -5.830, -5.601
495 0.1016: -6.821, -7.452
496 0.1035: -7.923, -9.733
497 0.1055: -9.144, -12.575
498 ...
499 0.4961: -96.614, -50.080
500 0.4980: -94.150, -47.668
501 0.5000: -93.426, -46.957
502 ... Figure 182
503 |#
504
505 ;-------------------------------------------------------------------------------
506 ;-------------------------------------------------------------------------------
507 ;;; Figure 200: Fejer's kernel for N = 4, 16 and 64
508 ;-------------------------------------------------------------------------------
509 ;-------------------------------------------------------------------------------
510 (multiple-value-bind (the-spec-wind freqs N-f)
511 (spectral-window-for-direct-spectral-estimate
512 4)
513 (format t "~&Figure 200, top plot ...")
514 (dotimes (i N-f)
515 (format t "~&~6,4F: ~8,4F"
516 (svref freqs i)
517 (svref the-spec-wind i)))
518 (format t "~&... Figure 200, top plot")
519 (values))
520
521 #|
522 ;;; Here is what is printed when the above form is evaluated:
523 Figure 200, top plot ...
524 0.0000: 6.0206
525 0.0625: 5.1644
526 0.1250: 2.3226
527 0.1875: -3.9257
528 0.2500: -100.0000
529 0.3125: -7.4278
530 0.3750: -5.3329
531 0.4375: -8.8624
532 0.5000: -100.0000
533 ... Figure 200, top plot
534 |#
535
536 ;;; Here we repeat the above calculation on a much finer grid of frequencies ...
537 (multiple-value-bind (the-spec-wind freqs N-f)
538 (spectral-window-for-direct-spectral-estimate
539 4
540 :n-nonzero-freqs 512)
541 (format t "~&Figure 200, top plot ...")
542 (dotimes (i 7)
543 (format t "~&~6,4F: ~8,4F"
544 (svref freqs i)
545 (svref the-spec-wind i)))
546 (format t "~&...")
547 (dotimes (i 7)
548 (format t "~&~6,4F: ~8,4F"
549 (svref freqs (+ i (- N-f 7)))
550 (svref the-spec-wind (+ i (- N-f 7)))))
551 (format t "~&... Figure 200, top plot")
552 (values))
553
554 #|
555 ;;; Here is what is printed when the above form is evaluated:
556 Figure 200, top plot ...
557 0.0000: 6.0206
558 0.0010: 6.0204
559 0.0020: 6.0198
560 0.0029: 6.0188
561 0.0039: 6.0173
562 0.0049: 6.0155
563 0.0059: 6.0132
564 ...
565 0.4941: -28.6858
566 0.4951: -30.2674
567 0.4961: -32.2040
568 0.4971: -34.7016
569 0.4980: -38.2225
570 0.4990: -44.2426
571 0.5000: -100.0000
572 ... Figure 200, top plot
573 |#
574
575 (multiple-value-bind (the-spec-wind freqs N-f)
576 (spectral-window-for-direct-spectral-estimate
577 16)
578 (format t "~&Figure 200, middle plot ...")
579 (dotimes (i 7)
580 (format t "~&~6,4F: ~8,4F"
581 (svref freqs i)
582 (svref the-spec-wind i)))
583 (format t "~&...")
584 (dotimes (i 7)
585 (format t "~&~6,4F: ~8,4F"
586 (svref freqs (+ i (- N-f 7)))
587 (svref the-spec-wind (+ i (- N-f 7)))))
588 (format t "~&... Figure 200, middle plot")
589 (values))
590
591 #|
592 ;;; Here is what is printed when the above form is evaluated:
593 Figure 200, middle plot ...
594 0.0000: 12.0412
595 0.0156: 11.1326
596 0.0312: 8.1328
597 0.0469: 1.6181
598 0.0625: -100.0000
599 0.0781: -2.7629
600 0.0937: -1.2977
601 ...
602 0.4062: -11.6589
603 0.4219: -14.7872
604 0.4375: -100.0000
605 0.4531: -14.9570
606 0.4687: -11.9993
607 0.4844: -15.0410
608 0.5000: -100.0000
609 ... Figure 200, middle plot
610 |#
611
612 (multiple-value-bind (the-spec-wind freqs N-f)
613 (spectral-window-for-direct-spectral-estimate
614 64)
615 (format t "~&Figure 200, bottom plot ...")
616 (dotimes (i 7)
617 (format t "~&~6,4F: ~8,4F"
618 (svref freqs i)
619 (svref the-spec-wind i)))
620 (format t "~&...")
621 (dotimes (i 7)
622 (format t "~&~6,4F: ~8,4F"
623 (svref freqs (+ i (- N-f 7)))
624 (svref the-spec-wind (+ i (- N-f 7)))))
625 (format t "~&... Figure 200, bottom plot")
626 (values))
627
628 #|
629 ;;; Here is what is printed when the above form is evaluated:
630 Figure 200, bottom plot ...
631 0.0000: 18.0618
632 0.0039: 17.1499
633 0.0078: 14.1403
634 0.0117: 7.6092
635 0.0156: -100.0000
636 0.0195: 3.1758
637 0.0234: 4.6048
638 ...
639 0.4766: -18.0382
640 0.4805: -21.0557
641 0.4844: -100.0000
642 0.4883: -21.0662
643 0.4922: -18.0592
644 0.4961: -21.0714
645 0.5000: -100.0000
646 ... Figure 200, bottom plot
647 |#
648
649 ;-------------------------------------------------------------------------------
650 ;-------------------------------------------------------------------------------
651 ;;; Figure 211: dpss data tapers and spectral windows for N = 64,
652 ;;; NW = 1 and 2
653 ;-------------------------------------------------------------------------------
654 ;-------------------------------------------------------------------------------
655 (let ((the-taper (dpss-data-taper!
656 (make-array 64 :initial-element 1.0)
657 :taper-parameter 1.0)))
658 (format t "~&sum of squares of taper elements = ~8,4F"
659 (sum-of-squares the-taper))
660 (format t "~&Figure 211, top left-hand plot ...")
661 (dotimes (i 7)
662 (format t "~&~2D: ~8,4F"
663 (1+ i)
664 (svref the-taper i)))
665 (format t "~&...")
666 (dotimes (i 7)
667 (format t "~&~2D: ~8,4F"
668 (+ (1+ i) (- 64 7))
669 (svref the-taper (+ i (- 64 7)))))
670 (format t "~&... Figure 211, top left-hand plot")
671 (values))
672
673 #|
674 ;;; Here is what is printed when the above form is evaluated:
675 sum of squares of taper elements = 1.0000
676 Figure 211, top left-hand plot ...
677 1: 0.0353
678 2: 0.0404
679 3: 0.0457
680 4: 0.0512
681 5: 0.0568
682 6: 0.0626
683 7: 0.0684
684 ...
685 58: 0.0684
686 59: 0.0626
687 60: 0.0568
688 61: 0.0512
689 62: 0.0457
690 63: 0.0404
691 64: 0.0353
692 ... Figure 211, top left-hand plot
693 |#
694
695 (let ((the-taper (dpss-data-taper!
696 (make-array 64 :initial-element 1.0)
697 :taper-parameter 2.0)))
698 (format t "~&sum of squares of taper elements = ~8,4F"
699 (sum-of-squares the-taper))
700 (format t "~&Figure 211, top right-hand plot ...")
701 (dotimes (i 7)
702 (format t "~&~2D: ~8,4F"
703 (1+ i)
704 (svref the-taper i)))
705 (format t "~&...")
706 (dotimes (i 7)
707 (format t "~&~2D: ~8,4F"
708 (+ (1+ i) (- 64 7))
709 (svref the-taper (+ i (- 64 7)))))
710 (format t "~&... Figure 211, top right-hand plot")
711 (values))
712
713 #|
714 ;;; Here is what is printed when the above form is evaluated:
715 sum of squares of taper elements = 1.0000
716 Figure 211, top right-hand plot ...
717 1: 0.0034
718 2: 0.0055
719 3: 0.0079
720 4: 0.0110
721 5: 0.0146
722 6: 0.0188
723 7: 0.0236
724 ...
725 58: 0.0236
726 59: 0.0188
727 60: 0.0146
728 61: 0.0110
729 62: 0.0079
730 63: 0.0055
731 64: 0.0034
732 ... Figure 211, top right-hand plot
733 |#
734
735 (multiple-value-bind (the-spec-wind freqs N-f)
736 (spectral-window-for-direct-spectral-estimate
737 64
738 :data-taper #'dpss-data-taper!
739 :data-taper-parameters 1.0)
740 (format t "~&Figure 211, bottom left-hand plot ...")
741 (dotimes (i 7)
742 (format t "~&~6,4F: ~8,4F"
743 (svref freqs i)
744 (svref the-spec-wind i)))
745 (format t "~&...")
746 (dotimes (i 7)
747 (format t "~&~6,4F: ~8,4F"
748 (svref freqs (+ i (- N-f 7)))
749 (svref the-spec-wind (+ i (- N-f 7)))))
750 (format t "~&... Figure 211, bottom left-hand plot")
751 (values))
752
753 #|
754 ;;; Here is what is printed when the above form is evaluated:
755 Figure 211, bottom left-hand plot ...
756 0.0000: 17.4690
757 0.0039: 16.8738
758 0.0078: 15.0172
759 0.0117: 11.6357
760 0.0156: 6.0027
761 0.0195: -4.8372
762 0.0234: -12.4591
763 ...
764 0.4766: -29.6551
765 0.4805: -32.6327
766 0.4844: -78.3605
767 0.4883: -32.7070
768 0.4922: -29.6760
769 0.4961: -32.6802
770 0.5000: -100.0000
771 ... Figure 211, bottom left-hand plot
772 |#
773
774 (multiple-value-bind (the-spec-wind freqs N-f)
775 (spectral-window-for-direct-spectral-estimate
776 64
777 :data-taper #'dpss-data-taper!
778 :data-taper-parameters 2.0)
779 (format t "~&Figure 211, bottom right-hand plot ...")
780 (dotimes (i 7)
781 (format t "~&~6,4F: ~8,4F"
782 (svref freqs i)
783 (svref the-spec-wind i)))
784 (format t "~&...")
785 (dotimes (i 7)
786 (format t "~&~6,4F: ~8,4F"
787 (svref freqs (+ i (- N-f 7)))
788 (svref the-spec-wind (+ i (- N-f 7)))))
789 (format t "~&... Figure 211, bottom right-hand plot")
790 (values))
791
792 #|
793 ;;; Here is what is printed when the above form is evaluated:
794 Figure 211, bottom right-hand plot ...
795 0.0000: 16.3413
796 0.0039: 15.9770
797 0.0078: 14.8693
798 0.0117: 12.9709
799 0.0156: 10.1904
800 0.0195: 6.3636
801 0.0234: 1.1838
802 ...
803 0.4766: -51.8112
804 0.4805: -54.6679
805 0.4844: -88.3687
806 0.4883: -54.9352
807 0.4922: -51.8305
808 0.4961: -54.8103
809 0.5000: -100.0000
810 ... Figure 211, bottom right-hand plot
811 |#
812
813 ;-------------------------------------------------------------------------------
814 ;-------------------------------------------------------------------------------
815 ;;; Figure 234: normalized cumulative periodogram for first 32 points
816 ;;; of AR(2) series
817 ;-------------------------------------------------------------------------------
818 ;-------------------------------------------------------------------------------
819 ;;; For this example we use the first 64 points of the AR(2) time series
820 ;;; shown in Figure 45 of the SAPA book.
821 (let* ((32-pt-ts (make-array 32 :displaced-to *AR-2*))
822 (cum-per-results (multiple-value-list
823 (cumulative-periodogram 32-pt-ts)))
824 (cum-per (elt cum-per-results 0))
825 (cum-per-freqs (elt cum-per-results 1))
826 (KS-results (subseq cum-per-results 2))
827 (the-periodogram (periodogram 32-pt-ts
828 :N-nonzero-freqs :Fourier
829 :sdf-transformation nil)))
830 (format t "~&Figure 234, left-hand plot, right-hand plot ...")
831 (dotimes (i 5)
832 (format t "~&~6,4F: ~8,4F ~8,4F"
833 (svref cum-per-freqs i)
834 (svref the-periodogram i)
835 (svref cum-per i)))
836 (format t "~&...")
837 (let ((N-cum-per (length cum-per)))
838 (dotimes (i 5)
839 (format t "~&~6,4F: ~8,4F ~8,4F"
840 (svref cum-per-freqs (+ i (- N-cum-per 5)))
841 (svref the-periodogram (+ i (- N-cum-per 5)))
842 (svref cum-per (+ i (- N-cum-per 5))))))
843 (apply #'format t "~&... Figure 234
844 Kolmogorov test statistic = ~A
845 index of maximum deviation = ~A
846 frequency of maximum deviation = ~A
847 quantile of Kolmogorov test statistic = ~A
848 reject/fail to reject null hypothesis = ~A
849 slope of upper and lower lines = ~A
850 intercept of upper line = ~A
851 intercept of lower line = ~A"
852 KS-results)
853 (values))
854
855 #|
856 ;;; Here is what is printed when the above form is evaluated:
857 Figure 234, left-hand plot, right-hand plot ...
858 0.0312: 0.2003 0.0055
859 0.0625: 2.9395 0.0866
860 0.0937: 7.8185 0.3022
861 0.1250: 7.9642 0.5219
862 0.1562: 5.6365 0.6773
863 ...
864 0.3437: 2.0928 0.9507
865 0.3750: 0.6617 0.9689
866 0.4062: 0.1675 0.9735
867 0.4375: 0.7901 0.9953
868 0.4687: 0.1696 1.0000
869 ... Figure 234
870 Kolmogorov test statistic = 0.3916144875655616
871 index of maximum deviation = 4
872 frequency of maximum deviation = 0.15625
873 quantile of Kolmogorov test statistic = 0.349
874 reject/fail to reject null hypothesis = reject
875 slope of upper and lower lines = 2.2857142857142856
876 intercept of upper line = 0.2775714285714286
877 intercept of lower line = -0.349
878 |#
879
880 ;-------------------------------------------------------------------------------
881 ;-------------------------------------------------------------------------------
882 ;;; Figure 265: Parzen lag window, smoothing window and spectral windows
883 ;;; using default and dpss data tapers
884 ;-------------------------------------------------------------------------------
885 ;-------------------------------------------------------------------------------
886 (progn
887 (format t "~&Figure 265, top left-hand plot ...")
888 (dotimes (i 7)
889 (format t "~&~2D: ~7,5F"
890 i
891 (parzen-lag-window i 37)))
892 (format t "~&...")
893 (dotimes (i 7)
894 (format t "~&~2D: ~7,5F"
895 (+ i 33)
896 (parzen-lag-window (+ i 33) 37)))
897 (format t "~&... Figure 265, top left-hand plot")
898 (values))
899
900 #|
901 ;;; Here is what is printed when the above form is evaluated:
902 0: 1.00000
903 1: 0.99574
904 2: 0.98342
905 3: 0.96375
906 4: 0.93746
907 5: 0.90524
908 6: 0.86781
909 ...
910 33: 0.00253
911 34: 0.00107
912 35: 0.00032
913 36: 0.00004
914 37: 0.00000
915 38: 0.00000
916 39: 0.00000
917 ... Figure 265, top left-hand plot
918 |#
919
920 (multiple-value-bind (the-smooth-wind freqs N-f)
921 (smoothing-window-for-lag-window-spectral-estimate
922 64
923 #'(lambda (tau)
924 (parzen-lag-window tau 37)))
925 (format t "~&Figure 265, top right-hand plot ...")
926 (dotimes (i 7)
927 (format t "~&~6,4F: ~8,4F"
928 (svref freqs i)
929 (svref the-smooth-wind i)))
930 (format t "~&...")
931 (dotimes (i 7)
932 (format t "~&~6,4F: ~8,4F"
933 (svref freqs (+ i (- N-f 7)))
934 (svref the-smooth-wind (+ i (- N-f 7)))))
935 (format t "~&... Figure 265, top right-hand plot")
936 (values))
937
938 #|
939 ;;; Here is what is printed when the above form is evaluated:
940 Figure 265, top right-hand plot ...
941 0.0000: 14.4326
942 0.0039: 14.2831
943 0.0078: 13.8316
944 0.0117: 13.0682
945 0.0156: 11.9756
946 0.0195: 10.5271
947 0.0234: 8.6825
948 ...
949 0.4766: -46.8807
950 0.4805: -45.6696
951 0.4844: -44.6896
952 0.4883: -44.5244
953 0.4922: -45.2033
954 0.4961: -46.3784
955 0.5000: -47.0461
956 ... Figure 265, top right-hand plot
957 |#
958
959 (multiple-value-bind (the-spec-wind freqs N-f)
960 (spectral-window-for-lag-window-spectral-estimate
961 64
962 #'(lambda (tau)
963 (parzen-lag-window tau 37)))
964 (format t "~&Figure 265, bottom left-hand plot ...")
965 (dotimes (i 7)
966 (format t "~&~6,4F: ~8,4F"
967 (svref freqs i)
968 (svref the-spec-wind i)))
969 (format t "~&...")
970 (dotimes (i 7)
971 (format t "~&~6,4F: ~8,4F"
972 (svref freqs (+ i (- N-f 7)))
973 (svref the-spec-wind (+ i (- N-f 7)))))
974 (format t "~&... Figure 265, bottom left-hand plot")
975 (values))
976
977 #|
978 ;;; Here is what is printed when the above form is evaluated:
979 Figure 265, bottom left-hand plot ...
980 0.0000: 13.8038
981 0.0039: 13.6754
982 0.0078: 13.2890
983 0.0117: 12.6417
984 0.0156: 11.7287
985 0.0195: 10.5453
986 0.0234: 9.0893
987 ...
988 0.4766: -21.0246
989 0.4805: -21.0307
990 0.4844: -21.0353
991 0.4883: -21.0396
992 0.4922: -21.0440
993 0.4961: -21.0475
994 0.5000: -21.0489
995 ... Figure 265, bottom left-hand plot
996 |#
997
998 (multiple-value-bind (the-spec-wind freqs N-f)
999 (spectral-window-for-lag-window-spectral-estimate
1000 64
1001 #'(lambda (tau)
1002 (parzen-lag-window tau 37))
1003 :data-taper #'dpss-data-taper!
1004 :data-taper-parameters 4.0)
1005 (format t "~&Figure 265, bottom right-hand plot ...")
1006 (dotimes (i 7)
1007 (format t "~&~6,4F: ~8,4F"
1008 (svref freqs i)
1009 (svref the-spec-wind i)))
1010 (format t "~&...")
1011 (dotimes (i 7)
1012 (format t "~&~6,4F: ~8,4F"
1013 (svref freqs (+ i (- N-f 7)))
1014 (svref the-spec-wind (+ i (- N-f 7)))))
1015 (format t "~&... Figure 265, bottom right-hand plot")
1016 (values))
1017
1018 #|
1019 ;;; Here is what is printed when the above form is evaluated:
1020 Figure 265, bottom right-hand plot ...
1021 0.0000: 13.2324
1022 0.0039: 13.1422
1023 0.0078: 12.8715
1024 0.0117: 12.4186
1025 0.0156: 11.7814
1026 0.0195: 10.9566
1027 0.0234: 9.9399
1028 ...
1029 0.4766: -45.0465
1030 0.4805: -45.2264
1031 0.4844: -45.3652
1032 0.4883: -45.4628
1033 0.4922: -45.5266
1034 0.4961: -45.5628
1035 0.5000: -45.5747
1036 ... Figure 265, bottom right-hand plot
1037 |#
1038
1039 ;-------------------------------------------------------------------------------
1040 ;-------------------------------------------------------------------------------
1041 ;;; Figure 296: periodogram and other direct spectral estimates
1042 ;;; for ocean wave data
1043 ;;; Note: in what follows we assume that the 1024 values for
1044 ;;; ocean wave data have been placed into the array pointed
1045 ;;; to by *ocean-wave-data* -- at the beginning of this file
1046 ;;; are examples of how to load numbers from an ASCII file
1047 ;;; into a Lisp vector.
1048 ;-------------------------------------------------------------------------------
1049 ;-------------------------------------------------------------------------------
1050 (sample-mean-and-variance *ocean-wave-data*)
1051 ;==> 209.103515625
1052 ; 143954.35647201538
1053
1054 (multiple-value-bind (the-periodogram freqs N-f)
1055 (periodogram
1056 *ocean-wave-data*
1057 :sampling-time 1/4)
1058 (format t "~&Figure 296, top plot ...")
1059 (dotimes (i 7)
1060 (format t "~&~6,4F: ~8,4F"
1061 (svref freqs i)
1062 (svref the-periodogram i)))
1063 (format t "~&...")
1064 (dotimes (i 7)
1065 (format t "~&~6,4F: ~8,4F"
1066 (svref freqs (+ i (- N-f 7)))
1067 (svref the-periodogram (+ i (- N-f 7)))))
1068 (format t "~&... Figure 296, top plot")
1069 (values))
1070
1071 #|
1072 ;;; Here is what is printed when the above form is evaluated:
1073 Figure 296, top plot ...
1074 0.0039: 47.8344
1075 0.0078: 39.8411
1076 0.0117: 43.0482
1077 0.0156: 43.9995
1078 0.0195: 50.3106
1079 0.0234: 54.7701
1080 0.0273: 54.0996
1081 ...
1082 1.9766: 13.2352
1083 1.9805: 11.5889
1084 1.9844: 11.8137
1085 1.9883: 12.3749
1086 1.9922: 14.7094
1087 1.9961: 12.9344
1088 2.0000: 12.2424
1089 ... Figure 296, top plot
1090 |#
1091
1092 (multiple-value-bind (the-direct-spectral-estimate freqs N-f)
1093 (direct-spectral-estimate
1094 *ocean-wave-data*
1095 :sampling-time 1/4
1096 :data-taper #'dpss-data-taper!
1097 :data-taper-parameters 1.0)
1098 (format t "~&Figure 296, second plot ...")
1099 (dotimes (i 7)
1100 (format t "~&~6,4F: ~8,4F"
1101 (svref freqs i)
1102 (svref the-direct-spectral-estimate i)))
1103 (format t "~&...")
1104 (dotimes (i 7)
1105 (format t "~&~6,4F: ~8,4F"
1106 (svref freqs (+ i (- N-f 7)))
1107 (svref the-direct-spectral-estimate (+ i (- N-f 7)))))
1108 (format t "~&... Figure 296, second plot")
1109 (values))
1110
1111 #|
1112 ;;; Here is what is printed when the above form is evaluated:
1113 Figure 296, second plot ...
1114 0.0039: 46.2041
1115 0.0078: 23.2523
1116 0.0117: 42.4188
1117 0.0156: 47.8951
1118 0.0195: 51.3775
1119 0.0234: 54.3370
1120 0.0273: 52.6524
1121 ...
1122 1.9766: 2.3258
1123 1.9805: -3.7210
1124 1.9844: -0.6480
1125 1.9883: -1.8448
1126 1.9922: 6.4357
1127 1.9961: -0.5324
1128 2.0000: -2.6338
1129 ... Figure 296, second plot
1130 |#
1131
1132 (multiple-value-bind (the-direct-spectral-estimate freqs N-f)
1133 (direct-spectral-estimate
1134 *ocean-wave-data*
1135 :sampling-time 1/4
1136 :data-taper #'dpss-data-taper!
1137 :data-taper-parameters 2.0)
1138 (format t "~&Figure 296, third plot ...")
1139 (dotimes (i 7)
1140 (format t "~&~6,4F: ~8,4F"
1141 (svref freqs i)
1142 (svref the-direct-spectral-estimate i)))
1143 (format t "~&...")
1144 (dotimes (i 7)
1145 (format t "~&~6,4F: ~8,4F"
1146 (svref freqs (+ i (- N-f 7)))
1147 (svref the-direct-spectral-estimate (+ i (- N-f 7)))))
1148 (format t "~&... Figure 296, third plot")
1149 (values))
1150
1151 #|
1152 ;;; Here is what is printed when the above form is evaluated:
1153 Figure 296, third plot ...
1154 0.0039: 43.8021
1155 0.0078: 36.5615
1156 0.0117: 42.9830
1157 0.0156: 49.5223
1158 0.0195: 51.8054
1159 0.0234: 53.3146
1160 0.0273: 50.7869
1161 ...
1162 1.9766: -11.0994
1163 1.9805: -8.6132
1164 1.9844: -5.4962
1165 1.9883: -6.0073
1166 1.9922: -0.1094
1167 1.9961: -14.6655
1168 2.0000: -16.8471
1169 ... Figure 296, third plot
1170 |#
1171
1172 (multiple-value-bind (the-direct-spectral-estimate freqs N-f)
1173 (direct-spectral-estimate
1174 *ocean-wave-data*
1175 :sampling-time 1/4
1176 :data-taper #'dpss-data-taper!
1177 :data-taper-parameters 4.0)
1178 (format t "~&Figure 296, bottom plot ...")
1179 (dotimes (i 7)
1180 (format t "~&~6,4F: ~8,4F"
1181 (svref freqs i)
1182 (svref the-direct-spectral-estimate i)))
1183 (format t "~&...")
1184 (dotimes (i 7)
1185 (format t "~&~6,4F: ~8,4F"
1186 (svref freqs (+ i (- N-f 7)))
1187 (svref the-direct-spectral-estimate (+ i (- N-f 7)))))
1188 (format t "~&... Figure 296, bottom plot")
1189 (values))
1190
1191 #|
1192 ;;; Here is what is printed when the above form is evaluated:
1193 Figure 296, bottom plot ...
1194 0.0039: 40.9897
1195 0.0078: 42.3745
1196 0.0117: 46.2331
1197 0.0156: 50.5033
1198 0.0195: 51.8640
1199 0.0234: 52.3855
1200 0.0273: 50.8260
1201 ...
1202 1.9766: -19.6861
1203 1.9805: -8.4228
1204 1.9844: -5.0683
1205 1.9883: -3.3882
1206 1.9922: -2.3835
1207 1.9961: -8.3321
1208 2.0000: -22.7452
1209 ... Figure 296, bottom plot
1210 |#
1211
1212 ;-------------------------------------------------------------------------------
1213 ;-------------------------------------------------------------------------------
1214 ;;; Figure 298: periodogram for ocean wave data
1215 ;;; at a finer grid of frequencies ...
1216 ;-------------------------------------------------------------------------------
1217 ;-------------------------------------------------------------------------------
1218 (multiple-value-bind (the-periodogram freqs N-f)
1219 (periodogram
1220 *ocean-wave-data*
1221 :N-nonzero-freqs :next-power-of-2
1222 :sampling-time 1/4)
1223 (format t "~&Figure 298 ...")
1224 (dotimes (i 7)
1225 (format t "~&~6,4F: ~8,4F"
1226 (svref freqs i)
1227 (svref the-periodogram i)))
1228 (format t "~&...")
1229 (dotimes (i 7)
1230 (format t "~&~6,4F: ~8,4F"
1231 (svref freqs (+ i (- N-f 7)))
1232 (svref the-periodogram (+ i (- N-f 7)))))
1233 (format t "~&... Figure 298")
1234 (values))
1235
1236 #|
1237 ;;; Here is what is printed when the above form is evaluated:
1238 Figure 298 ...
1239 0.0020: 50.5514
1240 0.0039: 47.8344
1241 0.0059: 37.1331
1242 0.0078: 39.8411
1243 0.0098: 39.0109
1244 0.0117: 43.0482
1245 0.0137: 39.6363
1246 ...
1247 1.9883: 12.3749
1248 1.9902: 0.1867
1249 1.9922: 14.7094
1250 1.9941: -5.7852
1251 1.9961: 12.9344
1252 1.9980: -10.8475
1253 2.0000: 12.2424
1254 ... Figure 298
1255 |#
1256
1257 ;-------------------------------------------------------------------------------
1258 ;-------------------------------------------------------------------------------
1259 ;;; Figure 301: lag window spectral estimates for ocean wave data
1260 ;;; based upon direct spectral estimate with NW=2 dpss data taper
1261 ;-------------------------------------------------------------------------------
1262 ;-------------------------------------------------------------------------------
1263 (multiple-value-bind (the-direct-spectral-estimate freqs N-f C_h the-acvs)
1264 (direct-spectral-estimate
1265 *ocean-wave-data*
1266 :sampling-time 1/4
1267 :return-acvs-p t
1268 :data-taper #'dpss-data-taper!
1269 :data-taper-parameters 2.0)
1270 (declare (ignore the-direct-spectral-estimate freqs N-f))
1271 (format t "~&C_h = ~5,3F" C_h)
1272 (format t "~& lag acvs")
1273 (dotimes (i 7)
1274 (format t "~&~4D: ~12,4F"
1275 i
1276 (svref the-acvs i)))
1277 (format t "~&...")
1278 (dotimes (i 7)
1279 (format t "~&~4D: ~12,4F"
1280 (+ i (- 1024 7))
1281 (svref the-acvs (+ i (- 1024 7)))))
1282 (setf *the-acvs* the-acvs
1283 *C_h* C_h)
1284 (values))
1285
1286 #|
1287 ;;; Here is what is printed when the above form is evaluated:
1288 C_h = 2.005
1289 lag acvs
1290 0: 143954.3565
1291 1: 138834.3641
1292 2: 124334.1234
1293 3: 102687.0485
1294 4: 76732.5318
1295 5: 49174.6362
1296 6: 22178.5321
1297 ...
1298 1017: -0.3232
1299 1018: -0.2984
1300 1019: -0.2523
1301 1020: -0.1911
1302 1021: -0.1258
1303 1022: -0.0679
1304 1023: -0.0250
1305 |#
1306
1307 (multiple-value-bind (Parzen-m-150 freqs N-f nu B_W)
1308 (lag-window-spectral-estimate
1309 *the-acvs*
1310 #'(lambda (lag)
1311 (parzen-lag-window lag 150))
1312 :sampling-time 1/4
1313 :C_h *C_h*)
1314 (format t "~&Figure 301, top plot ...")
1315 (dotimes (i 7)
1316 (format t "~&~6,4F: ~8,4F"
1317 (svref freqs i)
1318 (svref Parzen-m-150 i)))
1319 (format t "~&...")
1320 (dotimes (i 7)
1321 (format t "~&~6,4F: ~8,4F"
1322 (svref freqs (+ i (- N-f 7)))
1323 (svref Parzen-m-150 (+ i (- N-f 7)))))
1324 (format t "~&... Figure 301, top plot")
1325 (format t "~&equivalent degrees of freedom = ~5,1F" nu)
1326 (format t "~&smoothing window bandwidth = ~6,4F" B_W)
1327 (values))
1328
1329 #|
1330 ;;; Here is what is printed when the above form is evaluated:
1331 Figure 301, top plot ...
1332 0.0039: 47.5847
1333 0.0078: 48.2507
1334 0.0117: 49.1616
1335 0.0156: 50.1528
1336 0.0195: 51.1108
1337 0.0234: 51.9681
1338 0.0273: 52.6864
1339 ...
1340 1.9766: -5.8887
1341 1.9805: -5.7822
1342 1.9844: -5.5922
1343 1.9883: -5.3665
1344 1.9922: -5.1577
1345 1.9961: -5.0117
1346 2.0000: -4.9595
1347 ... Figure 301, top plot
1348 equivalent degrees of freedom = 12.6
1349 smoothing window bandwidth = 0.0494
1350 |#
1351
1352 (multiple-value-bind (Parzen-m-55 freqs N-f nu B_W)
1353 (lag-window-spectral-estimate
1354 *the-acvs*
1355 #'(lambda (lag)
1356 (parzen-lag-window lag 55))
1357 :sampling-time 1/4
1358 :C_h *C_h*)
1359 (format t "~&Figure 301, middle plot ...")
1360 (dotimes (i 7)
1361 (format t "~&~6,4F: ~8,4F"
1362 (svref freqs i)
1363 (svref Parzen-m-55 i)))
1364 (format t "~&...")
1365 (dotimes (i 7)
1366 (format t "~&~6,4F: ~8,4F"
1367 (svref freqs (+ i (- N-f 7)))
1368 (svref Parzen-m-55 (+ i (- N-f 7)))))
1369 (format t "~&... Figure 301, middle plot")
1370 (format t "~&equivalent degrees of freedom = ~5,1F" nu)
1371 (format t "~&smoothing window bandwidth = ~6,4F" B_W)
1372 (values))
1373
1374 #|
1375 ;;; Here is what is printed when the above form is evaluated:
1376 Figure 301, middle plot ...
1377 0.0039: 51.7292
1378 0.0078: 51.7425
1379 0.0117: 51.7636
1380 0.0156: 51.7913
1381 0.0195: 51.8242
1382 0.0234: 51.8604
1383 0.0273: 51.8985
1384 ...
1385 1.9766: -1.8860
1386 1.9805: -1.8645
1387 1.9844: -1.8520
1388 1.9883: -1.8457
1389 1.9922: -1.8433
1390 1.9961: -1.8428
1391 2.0000: -1.8428
1392 ... Figure 301, middle plot
1393 equivalent degrees of freedom = 34.4
1394 smoothing window bandwidth = 0.1349
1395 |#
1396
1397 (multiple-value-bind (Daniell-m-30 freqs N-f nu B_W)
1398 (lag-window-spectral-estimate
1399 *the-acvs*
1400 #'(lambda (lag)
1401 (daniell-lag-window lag 30))
1402 :sampling-time 1/4
1403 :C_h *C_h*)
1404 (format t "~&Figure 301, bottom plot ...")
1405 (dotimes (i 7)
1406 (format t "~&~6,4F: ~8,4F"
1407 (svref freqs i)
1408 (svref Daniell-m-30 i)))
1409 (format t "~&...")
1410 (dotimes (i 7)
1411 (format t "~&~6,4F: ~8,4F"
1412 (svref freqs (+ i (- N-f 7)))
1413 (svref Daniell-m-30 (+ i (- N-f 7)))))
1414 (format t "~&... Figure 301, bottom plot")
1415 (format t "~&equivalent degrees of freedom = ~5,1F" nu)
1416 (format t "~&smoothing window bandwidth = ~6,4F" B_W)
1417 (values))
1418
1419 #|
1420 ;;; Here is what is printed when the above form is evaluated:
1421 Figure 301, bottom plot ...
1422 0.0039: 52.1040
1423 0.0078: 52.2497
1424 0.0117: 52.4965
1425 0.0156: 52.5851
1426 0.0195: 52.4657
1427 0.0234: 52.0096
1428 0.0273: 51.6623
1429 ...
1430 1.9766: -6.0432
1431 1.9805: -5.9070
1432 1.9844: -5.8789
1433 1.9883: -5.7424
1434 1.9922: -5.4674
1435 1.9961: -5.2132
1436 2.0000: -5.1368
1437 ... Figure 301, bottom plot
1438 equivalent degrees of freedom = 34.1
1439 smoothing window bandwidth = 0.1337
1440 ;;; Note: page 302 of the SAPA book gives the degrees of freedom
1441 ;;; for these three cases as 13, 35 and 35 -- the above
1442 ;;; calculations differs slightly because we have obtained
1443 ;;; C_h using Equation (251b) rather from Table 248.
1444 |#
1445
1446 ;;; Here we show how to compute and relocate the spectral window
1447 ;;; shown in the lower left-hand corner of the top plot ...
1448 (multiple-value-bind (the-smooth-wind freqs N-f)
1449 (smoothing-window-for-lag-window-spectral-estimate
1450 1024
1451 #'(lambda (tau)
1452 (parzen-lag-window tau 150))
1453 :sampling-time 1/4
1454 :N-nonzero-freqs :half-next-power-of-2)
1455 (let ((dB-offset (+ (- (svref the-smooth-wind 0)) -10.0)))
1456 (setf the-smooth-wind (one-sided-sdf->two-sided-sdf the-smooth-wind)
1457 freqs (one-sided-freq->two-sided-freq freqs)
1458 N-f (1- (* 2 N-f)))
1459 (x+b! the-smooth-wind dB-offset)
1460 (x+b! freqs 0.25)
1461 (dotimes (i N-f (values))
1462 (if (<= 0.2 (svref freqs i) 0.3)
1463 (format t "~&~6,4F: ~8,4F"
1464 (svref freqs i)
1465 (svref the-smooth-wind i))))))
1466
1467 #|
1468 ;;; Here is what is printed when the above form is evaluated:
1469 0.2031: -44.8539
1470 0.2070: -35.7970
1471 0.2109: -29.5888
1472 0.2148: -24.9154
1473 0.2187: -21.2410
1474 0.2227: -18.2942
1475 0.2266: -15.9187
1476 0.2305: -14.0171
1477 0.2344: -12.5259
1478 0.2383: -11.4022
1479 0.2422: -10.6175
1480 0.2461: -10.1536
1481 0.2500: -10.0000
1482 0.2539: -10.1536
1483 0.2578: -10.6175
1484 0.2617: -11.4022
1485 0.2656: -12.5259
1486 0.2695: -14.0171
1487 0.2734: -15.9187
1488 0.2773: -18.2942
1489 0.2812: -21.2410
1490 0.2852: -24.9154
1491 0.2891: -29.5888
1492 0.2930: -35.7970
1493 0.2969: -44.8539
1494 |#
1495
1496 ;-------------------------------------------------------------------------------
1497 ;-------------------------------------------------------------------------------
1498 ;;; Figure 305: periodograms for rough and smooth ice profile data
1499 ;;; Note: in what follows we assume that the 1121 values for
1500 ;;; rough ice profile data have been placed into the array pointed
1501 ;;; to by *rough-ice-profile-data* and that the 288 values for
1502 ;;; smooth ice profile data have been placed into the array pointed
1503 ;;; to by *smooth-ice-profile-data* -- at the beginning of this file
1504 ;;; are examples of how to load numbers from an ASCII file
1505 ;;; into a Lisp vector.
1506 ;-------------------------------------------------------------------------------
1507 ;-------------------------------------------------------------------------------
1508 (sample-mean *rough-ice-profile-data*) ;==> -19.916859946476382
1509
1510 (multiple-value-bind (the-periodogram freqs N-f)
1511 (periodogram
1512 *rough-ice-profile-data*
1513 :sampling-time 1.7712)
1514 (format t "~&Figure 305, top plot ...")
1515 (dotimes (i 7)
1516 (format t "~&~6,4F: ~8,4F"
1517 (svref freqs i)
1518 (svref the-periodogram i)))
1519 (format t "~&...")
1520 (dotimes (i 7)
1521 (format t "~&~6,4F: ~8,4F"
1522 (svref freqs (+ i (- N-f 7)))
1523 (svref the-periodogram (+ i (- N-f 7)))))
1524 (format t "~&... Figure 305, top plot")
1525 (values))
1526
1527 #|
1528 ;;; Here is what is printed when the above form is evaluated:
1529 Figure 305, top plot ...
1530 0.0003: 19.3206
1531 0.0006: 29.8243
1532 0.0008: 33.4726
1533 0.0011: 33.9632
1534 0.0014: 34.7072
1535 0.0017: 35.0793
1536 0.0019: 36.2898
1537 ...
1538 0.2806: -0.3305
1539 0.2809: 0.2546
1540 0.2812: -1.4455
1541 0.2815: 2.7054
1542 0.2817: 4.2203
1543 0.2820: -5.3819
1544 0.2823: -0.1654
1545 ... Figure 305, top plot
1546 |#
1547
1548 (sample-mean *smooth-ice-profile-data*) ;==> -8.934027777777777
1549
1550 (multiple-value-bind (the-periodogram freqs N-f)
1551 (periodogram
1552 *smooth-ice-profile-data*
1553 :sampling-time 1.7712)
1554 (format t "~&Figure 305, bottom plot ...")
1555 (dotimes (i 7)
1556 (format t "~&~6,4F: ~8,4F"
1557 (svref freqs i)
1558 (svref the-periodogram i)))
1559 (format t "~&...")
1560 (dotimes (i 7)
1561 (format t "~&~6,4F: ~8,4F"
1562 (svref freqs (+ i (- N-f 7)))
1563 (svref the-periodogram (+ i (- N-f 7)))))
1564 (format t "~&... Figure 305, bottom plot")
1565 (values))
1566
1567 #|
1568 ;;; Here is what is printed when the above form is evaluated:
1569 Figure 305, bottom plot ...
1570 0.0011: 5.7414
1571 0.0022: -0.3703
1572 0.0033: 0.9999
1573 0.0044: 4.3387
1574 0.0055: 7.5413
1575 0.0066: 4.8875
1576 0.0077: 1.8933
1577 ...
1578 0.2757: -11.9119
1579 0.2768: -17.4537
1580 0.2779: -22.4416
1581 0.2790: -15.1718
1582 0.2801: -9.8917
1583 0.2812: -10.6674
1584 0.2823: -16.0906
1585 ... Figure 305, bottom plot
1586 |#
1587
1588 ;-------------------------------------------------------------------------------
1589 ;-------------------------------------------------------------------------------
1590 ;;; Figure 307: lag window spectral estimates for rough and smooth
1591 ;;; ice profile data
1592 ;-------------------------------------------------------------------------------
1593 ;-------------------------------------------------------------------------------
1594 (let ((rough-acvs (acvs *rough-ice-profile-data*)))
1595 (multiple-value-bind (Parzen-m-88 freqs N-f nu B_W)
1596 (lag-window-spectral-estimate
1597 rough-acvs
1598 #'(lambda (lag)
1599 (parzen-lag-window lag 88))
1600 :sampling-time 1.7712)
1601 (format t "~&Figure 307, top plot ...")
1602 (dotimes (i 7)
1603 (format t "~&~6,4F: ~8,4F"
1604 (svref freqs i)
1605 (svref Parzen-m-88 i)))
1606 (format t "~&...")
1607 (dotimes (i 7)
1608 (format t "~&~6,4F: ~8,4F"
1609 (svref freqs (+ i (- N-f 7)))
1610 (svref Parzen-m-88 (+ i (- N-f 7)))))
1611 (format t "~&... Figure 307, top plot")
1612 (format t "~&equivalent degrees of freedom = ~5,1F" nu)
1613 (format t "~&smoothing window bandwidth = ~6,4F" B_W)
1614 (format t "~&time series bandwidth = ~6,4F"
1615 (time-series-bandwidth rough-acvs :sampling-time 1.7712))))
1616 #|
1617 ;;; Here is what is printed when the above form is evaluated:
1618 Figure 307, top plot ...
1619 0.0003: 32.2615
1620 0.0006: 32.2518
1621 0.0008: 32.2357
1622 0.0011: 32.2130
1623 0.0014: 32.1838
1624 0.0017: 32.1480
1625 0.0019: 32.1056
1626 ...
1627 0.2806: 0.6920
1628 0.2809: 0.6401
1629 0.2812: 0.5964
1630 0.2815: 0.5615
1631 0.2817: 0.5360
1632 0.2820: 0.5206
1633 0.2823: 0.5154
1634 ... Figure 307, top plot
1635 equivalent degrees of freedom = 47.2
1636 smoothing window bandwidth = 0.0119
1637 time series bandwidth = 0.0239
1638 |#
1639
1640 (let ((smooth-acvs (acvs *smooth-ice-profile-data*)))
1641 (multiple-value-bind (Parzen-m-20 freqs N-f nu B_W)
1642 (lag-window-spectral-estimate
1643 smooth-acvs
1644 #'(lambda (lag)
1645 (parzen-lag-window lag 20))
1646 :sampling-time 1.7712)
1647 (format t "~&Figure 307, middle plot ...")
1648 (dotimes (i 7)
1649 (format t "~&~6,4F: ~8,4F"
1650 (svref freqs i)
1651 (svref Parzen-m-20 i)))
1652 (format t "~&...")
1653 (dotimes (i 7)
1654 (format t "~&~6,4F: ~8,4F"
1655 (svref freqs (+ i (- N-f 7)))
1656 (svref Parzen-m-20 (+ i (- N-f 7)))))
1657 (format t "~&... Figure 307, middle plot")
1658 (format t "~&equivalent degrees of freedom = ~5,1F" nu)
1659 (format t "~&smoothing window bandwidth = ~6,4F" B_W)
1660 (format t "~&time series bandwidth = ~6,4F"
1661 (time-series-bandwidth smooth-acvs :sampling-time 1.7712)))
1662 (values))
1663
1664 #|
1665 ;;; Here is what is printed when the above form is evaluated:
1666 Figure 307, middle plot ...
1667 0.0011: 3.2732
1668 0.0022: 3.2655
1669 0.0033: 3.2527
1670 0.0044: 3.2349
1671 0.0055: 3.2120
1672 0.0066: 3.1841
1673 0.0077: 3.1513
1674 ...
1675 0.2757: -13.2314
1676 0.2768: -13.2507
1677 0.2779: -13.2666
1678 0.2790: -13.2789
1679 0.2801: -13.2877
1680 0.2812: -13.2930
1681 0.2823: -13.2948
1682 ... Figure 307, middle plot
1683 equivalent degrees of freedom = 53.4
1684 smoothing window bandwidth = 0.0523
1685 time series bandwidth = 0.1069
1686 |#
1687
1688 ;-------------------------------------------------------------------------------
1689 ;-------------------------------------------------------------------------------
1690 ;;; Figure 312: WOSA spectral estimates of AR(2) series
1691 ;;; Note: in what follows we assume that the 1024 values for AR(20 series
1692 ;;; have been placed into a vector pointed to by *AR-2* -- at the
1693 ;;; beginning of this file are examples of how to load numbers
1694 ;;; from an ASCII file into a Lisp vector.
1695 ;-------------------------------------------------------------------------------
1696 ;-------------------------------------------------------------------------------
1697 ;;; block size = 8 points
1698 (multiple-value-bind (wosa freqs dof)
1699 (wosa-spectral-estimate
1700 *AR-2*
1701 8
1702 :oversampling-factor 8)
1703 (format t "~&Figure 307, top plot ...")
1704 (dotimes (i 7)
1705 (format t "~&~6,4F: ~8,4F"
1706 (svref freqs i)
1707 (svref wosa i)))
1708 (format t "~&...")
1709 (let ((N-f (length freqs)))
1710 (dotimes (i 7)
1711 (format t "~&~6,4F: ~8,4F"
1712 (svref freqs (+ i (- N-f 7)))
1713 (svref wosa (+ i (- N-f 7))))))
1714 (format t "~&... Figure 307, top plot")
1715 (format t "~&equivalent degrees of freedom = ~5,1F" dof)
1716 (values))
1717
1718 #|
1719 ;;; Here is what is printed when the above form is evaluated:
1720 Figure 307, top plot ...
1721 0.0000: 3.1104
1722 0.0156: 3.1860
1723 0.0312: 3.3992
1724 0.0469: 3.7143
1725 0.0625: 4.0842
1726 0.0781: 4.4621
1727 0.0937: 4.8075
1728 ...
1729 0.4062: -4.2055
1730 0.4219: -4.5725
1731 0.4375: -4.8838
1732 0.4531: -5.1387
1733 0.4687: -5.3309
1734 0.4844: -5.4517
1735 0.5000: -5.4930
1736 ... Figure 307, top plot
1737 equivalent degrees of freedom = 453.2
1738 |#
1739
1740 ;;; block size = 16 points
1741 (multiple-value-bind (wosa freqs dof)
1742 (wosa-spectral-estimate
1743 *AR-2*
1744 16
1745 :oversampling-factor 8)
1746 (format t "~&Figure 307, middle plot ...")
1747 (dotimes (i 7)
1748 (format t "~&~6,4F: ~8,4F"
1749 (svref freqs i)
1750 (svref wosa i)))
1751 (format t "~&...")
1752 (let ((N-f (length freqs)))
1753 (dotimes (i 7)
1754 (format t "~&~6,4F: ~8,4F"
1755 (svref freqs (+ i (- N-f 7)))
1756 (svref wosa (+ i (- N-f 7))))))
1757 (format t "~&... Figure 307, middle plot")
1758 (format t "~&equivalent degrees of freedom = ~5,1F" dof)
1759 (values))
1760
1761 #|
1762 ;;; Here is what is printed when the above form is evaluated:
1763 Figure 307, middle plot ...
1764 0.0000: 2.7104
1765 0.0078: 2.7175
1766 0.0156: 2.7407
1767 0.0234: 2.7860
1768 0.0312: 2.8624
1769 0.0391: 2.9807
1770 0.0469: 3.1503
1771 ...
1772 0.4531: -5.7505
1773 0.4609: -5.9087
1774 0.4687: -6.0551
1775 0.4766: -6.1807
1776 0.4844: -6.2774
1777 0.4922: -6.3384
1778 0.5000: -6.3593
1779 ... Figure 307, middle plot
1780 equivalent degrees of freedom = 233.8
1781 |#
1782
1783 ;;; block size = 64 points
1784 (multiple-value-bind (wosa freqs dof)
1785 (wosa-spectral-estimate
1786 *AR-2*
1787 64
1788 :oversampling-factor 8)
1789 (format t "~&Figure 307, bottom plot ...")
1790 (dotimes (i 7)
1791 (format t "~&~6,4F: ~8,4F"
1792 (svref freqs i)
1793 (svref wosa i)))
1794 (format t "~&...")
1795 (let ((N-f (length freqs)))
1796 (dotimes (i 7)
1797 (format t "~&~6,4F: ~8,4F"
1798 (svref freqs (+ i (- N-f 7)))
1799 (svref wosa (+ i (- N-f 7))))))
1800 (format t "~&... Figure 307, bottom plot")
1801 (format t "~&equivalent degrees of freedom = ~5,1F" dof)
1802 (values))
1803
1804 #|
1805 ;;; Here is what is printed when the above form is evaluated:
1806 Figure 307, bottom plot ...
1807 0.0000: 2.6662
1808 0.0020: 2.6924
1809 0.0039: 2.7645
1810 0.0059: 2.8654
1811 0.0078: 2.9718
1812 0.0098: 3.0595
1813 0.0117: 3.1077
1814 ...
1815 0.4883: -8.0203
1816 0.4902: -8.1908
1817 0.4922: -8.3616
1818 0.4941: -8.5197
1819 0.4961: -8.6494
1820 0.4980: -8.7351
1821 0.5000: -8.7650
1822 ... Figure 307, bottom plot
1823 equivalent degrees of freedom = 58.5
1824 ;;; Note: the above values for the equivalent degrees of freedom
1825 ;;; differ from what is stated on page 311 of the SAPA book
1826 ;;; (the book has 483.3, 240.7 and 58.8 rather than 453.2, 233.8 and 58.5).
1827 ;;; This discrepancy is due to the fact that wosa-spectral-estimate
1828 ;;; uses Equation (292b), whereas the values quoted in the book
1829 ;;; are based upon Equation (294) (an approximation to Equation (292b)).
1830 ;;; Also, the spectra plotted in Figure 312 were actually computed
1831 ;;; using a slightly different definition for the Hanning data taper
1832 ;;; (essentially the one in Bloomfield's book). The spectra computed
1833 ;;; above thus are slightly different from what are plotted in the SAPA
1834 ;;; book (the biggest changes are a 0.5 dB difference at f = 0 and
1835 ;;; f = 0.5 for the 8 point block size, i.e., the top plot of Figure 312).
1836 |#
1837
1838 ;-------------------------------------------------------------------------------
1839 ;-------------------------------------------------------------------------------
1840 ;;; Figure 373: multitaper spectral estimates for ocean wave data
1841 ;-------------------------------------------------------------------------------
1842 ;-------------------------------------------------------------------------------
1843 (multiple-value-bind (dpss-NW-4 eigenvalues-NW-4)
1844 (dpss-tapers-tri-diag
1845 (length *ocean-wave-data*)
1846 7
1847 :taper-parameter 4.0
1848 :compute-true-eigenvalues-p t)
1849 (setf *eigenvalues-NW-4* eigenvalues-NW-4)
1850 (format t "~&order eigenvalue")
1851 (dotimes (k 7)
1852 (format t "~&k = ~1D: ~F" k (svref eigenvalues-NW-4 k)))
1853 ;;; We compute 7 eigenspectra but use only 6 of them to form
1854 ;;; the simple multitaper spectral estimate in plot (a).
1855 ;;; Note: for Figure 373, the eigenspectra are computed by setting
1856 ;;; the keyword options recenter-after-tapering-p and
1857 ;;; restore-power-option-p to nil. This produces a spectral
1858 ;;; estimate in agreement with the equations in Chapter 7,
1859 ;;; but arguably it might be better to just leave these keywords
1860 ;;; set to their default values of t.
1861 (multiple-value-bind
1862 (simple-multitaper-NW-4-0-to-6 freqs N-f eigenspectra-NW-4)
1863 (multitaper-spectral-estimate
1864 *ocean-wave-data*
1865 dpss-NW-4
1866 :sampling-time 1/4
1867 :recenter-after-tapering-p nil
1868 :restore-power-option-p nil)
1869 (declare (ignore simple-multitaper-NW-4-0-to-6))
1870 (setf *eigenspectra-NW-4* eigenspectra-NW-4
1871 *freqs* freqs)
1872 (let ((simple-multitaper-NW-4-0-to-5
1873 (eigenspectra->multitaper-spectral-estimate
1874 eigenspectra-NW-4 :N-eigenspectra 6)))
1875 (format t "~&Figure 373, plot (a) ...")
1876 (dotimes (i 7)
1877 (format t "~&~6,4F: ~8,4F"
1878 (svref freqs i)
1879 (svref simple-multitaper-NW-4-0-to-5 i)))
1880 (format t "~&...")
1881 (dotimes (i 7)
1882 (format t "~&~6,4F: ~8,4F"
1883 (svref freqs (+ i (- N-f 7)))
1884 (svref simple-multitaper-NW-4-0-to-5 (+ i (- N-f 7)))))
1885 (format t "~&... Figure 373, plot (a)")
1886 (values))))
1887
1888 #|
1889 ;;; Here is what is printed when the above form is evaluated:
1890 order eigenvalue
1891 k = 0: 0.9999999997056523
1892 k = 1: 0.9999999723287881
1893 k = 2: 0.9999987902598706
1894 k = 3: 0.9999675626065103
1895 k = 4: 0.9994101803916527
1896 k = 5: 0.9925053051988856
1897 k = 6: 0.9366554082073846
1898 Figure 373, plot (a) ...
1899 0.0039: 45.0265
1900 0.0078: 46.7926
1901 0.0117: 48.5376
1902 0.0156: 50.5297
1903 0.0195: 51.1797
1904 0.0234: 53.3410
1905 0.0273: 53.9492
1906 ...
1907 1.9766: -5.0610
1908 1.9805: -4.9186
1909 1.9844: -4.7007
1910 1.9883: -3.5230
1911 1.9922: -3.0773
1912 1.9961: -3.3919
1913 2.0000: -3.3314
1914 ... Figure 373, plot (a)
1915 |#
1916
1917 (multiple-value-bind (adaptive-multitaper-NW-4 dof)
1918 (eigenspectra->adaptive-multitaper-spectral-estimate
1919 *eigenspectra-NW-4*
1920 *eigenvalues-NW-4*
1921 (sample-variance *ocean-wave-data*)
1922 :sampling-time 1/4)
1923 (multiple-value-bind (upper lower)
1924 (create-ci-for-amt-sdf-estimate
1925 adaptive-multitaper-NW-4
1926 dof)
1927 (format t "~&Figure 373, plot (c) ... plot (d) ...")
1928 (dotimes (i 7)
1929 (format t "~&~6,4F: ~8,4F ~8,4F ~8,4F ~8,4F"
1930 (svref *freqs* i)
1931 (svref upper i)
1932 (svref adaptive-multitaper-NW-4 i)
1933 (svref lower i)
1934 (svref dof i)))
1935 (format t "~&...")
1936 (let ((N-f (length *freqs*)))
1937 (dotimes (i 7)
1938 (format t "~&~6,4F: ~8,4F ~8,4F ~8,4F ~8,4F"
1939 (svref *freqs* (+ i (- N-f 7)))
1940 (svref upper (+ i (- N-f 7)))
1941 (svref adaptive-multitaper-NW-4 (+ i (- N-f 7)))
1942 (svref lower (+ i (- N-f 7)))
1943 (svref dof (+ i (- N-f 7))))))
1944 (format t "~&... Figure 373, plots (c) and (d)")
1945 (values)))
1946
1947 #|
1948 ;;; Here is what is printed when the above form is evaluated:
1949 Figure 373, plot (c) ... plot (d) ...
1950 0.0039: 48.5041 44.5444 41.8350 13.9860
1951 0.0078: 50.9711 47.0137 44.3053 13.9985
1952 0.0117: 53.7760 49.8188 47.1104 13.9996
1953 0.0156: 54.1360 50.1787 47.4704 13.9994
1954 0.0195: 55.6164 51.6589 48.9505 13.9982
1955 0.0234: 56.6806 52.7230 50.0145 13.9973
1956 0.0273: 57.6636 53.7058 50.9973 13.9966
1957 ...
1958 1.9766: 1.0163 -5.8131 -9.6216 6.0346
1959 1.9805: -0.5008 -7.4550 -11.3023 5.8786
1960 1.9844: 2.9011 -3.7747 -7.5345 6.2375
1961 1.9883: 2.5575 -4.1425 -7.9100 6.2047
1962 1.9922: 3.3604 -3.2706 -7.0160 6.2991
1963 1.9961: 2.2524 -4.4749 -8.2511 6.1681
1964 2.0000: 4.0624 -2.5117 -6.2389 6.3789
1965 ... Figure 373, plots (c) and (d)
1966 |#
1967
1968 (multiple-value-bind (simple-multitaper-NW-6-0-to-9 freqs N-f)
1969 (multitaper-spectral-estimate
1970 *ocean-wave-data*
1971 (dpss-tapers-tri-diag
1972 (length *ocean-wave-data*)
1973 10
1974 :taper-parameter 6.0)
1975 :sampling-time 1/4
1976 :recenter-after-tapering-p nil
1977 :restore-power-option-p nil)
1978 (format t "~&Figure 373, plot (b) ...")
1979 (dotimes (i 7)
1980 (format t "~&~6,4F: ~8,4F"
1981 (svref freqs i)
1982 (svref simple-multitaper-NW-6-0-to-9 i)))
1983 (format t "~&...")
1984 (dotimes (i 7)
1985 (format t "~&~6,4F: ~8,4F"
1986 (svref freqs (+ i (- N-f 7)))
1987 (svref simple-multitaper-NW-6-0-to-9 (+ i (- N-f 7)))))
1988 (format t "~&... Figure 373, plot (b)")
1989 (values))
1990
1991 #|
1992 Figure 373, plot (b) ...
1993 0.0039: 48.0430
1994 0.0078: 49.0666
1995 0.0117: 49.7126
1996 0.0156: 51.5444
1997 0.0195: 51.9033
1998 0.0234: 53.5943
1999 0.0273: 53.7669
2000 ...
2001 1.9766: -4.2660
2002 1.9805: -4.8461
2003 1.9844: -4.4768
2004 1.9883: -4.8493
2005 1.9922: -3.9439
2006 1.9961: -2.6031
2007 2.0000: -1.7041
2008 ... Figure 373, plot (b)
2009 |#
2010
2011 ;-------------------------------------------------------------------------------
2012 ;-------------------------------------------------------------------------------
2013 ;;; Figures 440-1: designing a prewhitening filter for ocean wave data
2014 ;-------------------------------------------------------------------------------
2015 ;-------------------------------------------------------------------------------
2016 ;;; estimation of AR(5) model using the Yule-Walker method,
2017 ;;; from which we compute a parametric sdf estimate ...
2018 (multiple-value-bind (ar-5-coeffs innovations-variance)
2019 (yule-walker-algorithm-given-data
2020 *ocean-wave-data*
2021 5)
2022 (multiple-value-bind (Y-W-5-sdf freqs N-f)
2023 (ar-coeffs->sdf ar-5-coeffs
2024 innovations-variance
2025 :sampling-time 1/4
2026 :return-frequencies-p t)
2027 (format t "~&Figure 440, thin curve, top plot ...")
2028 (dotimes (i 7)
2029 (format t "~&~6,4F: ~8,4F"
2030 (svref freqs i)
2031 (svref Y-W-5-sdf i)))
2032 (format t "~&...")
2033 (dotimes (i 7)
2034 (format t "~&~6,4F: ~8,4F"
2035 (svref freqs (+ i (- N-f 7)))
2036 (svref Y-W-5-sdf (+ i (- N-f 7)))))
2037 (format t "~&... Figure 440, thin curve, top plot")
2038 (values)))
2039
2040 #|
2041 ;;; Here is what is printed when the above form is evaluated:
2042 Figure 440, thin curve, top plot ...
2043 0.0000: 51.2480
2044 0.0078: 51.2659
2045 0.0156: 51.3197
2046 0.0234: 51.4097
2047 0.0312: 51.5365
2048 0.0391: 51.7008
2049 0.0469: 51.9034
2050 ...
2051 1.9531: 10.5844
2052 1.9609: 10.5978
2053 1.9687: 10.6089
2054 1.9766: 10.6175
2055 1.9844: 10.6237
2056 1.9922: 10.6275
2057 2.0000: 10.6287
2058 ... Figure 440, thin curve, top plot
2059 |#
2060
2061 ;;; estimation of AR(27) model using Burg's algorithm,
2062 ;;; from which we compute a parametric sdf estimate ...
2063 (multiple-value-bind (ar-27-coeffs innovations-variance)
2064 (Burg-algorithm
2065 *ocean-wave-data*
2066 27)
2067 (multiple-value-bind (Burg-27-sdf freqs N-f)
2068 (ar-coeffs->sdf ar-27-coeffs
2069 innovations-variance
2070 :sampling-time 1/4
2071 :return-frequencies-p t)
2072 (format t "~&Figure 440, thin curve, bottom plot ...")
2073 (dotimes (i 7)
2074 (format t "~&~6,4F: ~8,4F"
2075 (svref freqs i)
2076 (svref Burg-27-sdf i)))
2077 (format t "~&...")
2078 (dotimes (i 7)
2079 (format t "~&~6,4F: ~8,4F"
2080 (svref freqs (+ i (- N-f 7)))
2081 (svref Burg-27-sdf (+ i (- N-f 7)))))
2082 (format t "~&... Figure 440, thin curve, bottom plot")
2083 (values)))
2084
2085 #|
2086 ;;; Here is what is printed when the above form is evaluated:
2087 Figure 440, thin curve, bottom plot ...
2088 0.0000: 51.8188
2089 0.0078: 51.8771
2090 0.0156: 52.0414
2091 0.0234: 52.2796
2092 0.0312: 52.5356
2093 0.0391: 52.7342
2094 0.0469: 52.7986
2095 ...
2096 1.9531: -6.2274
2097 1.9609: -5.8871
2098 1.9687: -5.5495
2099 1.9766: -5.2402
2100 1.9844: -4.9885
2101 1.9922: -4.8231
2102 2.0000: -4.7653
2103 ... Figure 440, thin curve, bottom plot
2104 |#
2105
2106 ;;; estimation of AR(5) model using Burg's algorithm,
2107 ;;; from which we [a] compute a parametric sdf estimate
2108 ;;; [b] compute the periodogram for the forward prediction errors
2109 ;;; [c] smooth the periodogram using a Parzen lag window and
2110 ;;; then postcolor this lag window estimate using the
2111 ;;; AR prewhitening filter
2112 (multiple-value-bind (ar-5-coeffs innovations-variance
2113 junk-1 junk-2
2114 forward-pred-errors)
2115 (Burg-algorithm *ocean-wave-data* 5)
2116 (declare (ignore junk-1 junk-2))
2117 (setf *ar-5-coeffs* ar-5-coeffs
2118 *forward-pred-errors* forward-pred-errors)
2119 (multiple-value-bind (Burg-5-sdf freqs N-f)
2120 (ar-coeffs->sdf ar-5-coeffs
2121 innovations-variance
2122 :sampling-time 1/4
2123 :return-frequencies-p t)
2124 (format t "~&Figure 440, thick curve, both plots ...")
2125 (dotimes (i 7)
2126 (format t "~&~6,4F: ~8,4F"
2127 (svref freqs i)
2128 (svref Burg-5-sdf i)))
2129 (format t "~&...")
2130 (dotimes (i 7)
2131 (format t "~&~6,4F: ~8,4F"
2132 (svref freqs (+ i (- N-f 7)))
2133 (svref Burg-5-sdf (+ i (- N-f 7)))))
2134 (format t "~&... Figure 440, thick curve, both plots")
2135 (values)))
2136
2137 #|
2138 ;;; Here is what is printed when the above form is evaluated:
2139 Figure 440, thick curve, both plots ...
2140 0.0000: 51.2758
2141 0.0078: 51.2919
2142 0.0156: 51.3403
2143 0.0234: 51.4215
2144 0.0312: 51.5361
2145 0.0391: 51.6849
2146 0.0469: 51.8690
2147 ...
2148 1.9531: -7.9316
2149 1.9609: -7.9411
2150 1.9687: -7.9489
2151 1.9766: -7.9549
2152 1.9844: -7.9592
2153 1.9922: -7.9618
2154 2.0000: -7.9627
2155 ... Figure 440, thick curve, both plots
2156 |#
2157
2158 (progn
2159 (format t "~&Figure 441, plot (a) ...")
2160 (dotimes (i 7)
2161 (format t "~&~4D: ~8,4F"
2162 (+ i 6)
2163 (svref *forward-pred-errors* i)))
2164 (format t "~&...")
2165 (let ((N-fpe (length *forward-pred-errors*)))
2166 (dotimes (i 7)
2167 (format t "~&~4D: ~8,4F"
2168 (+ i 6 (- N-fpe 7))
2169 (svref *forward-pred-errors* (+ i (- N-fpe 7))))))
2170 (format t "~&... Figure 441, plot (a)")
2171 (values))
2172
2173 #|
2174 ;;; Here is what is printed when the above form is evaluated:
2175 Figure 441, plot (a) ...
2176 6: 6.6713
2177 7: 14.7491
2178 8: -5.1810
2179 9: 2.9906
2180 10: -9.3390
2181 11: 21.8462
2182 12: -0.5174
2183 ...
2184 1018: -12.8457
2185 1019: 0.9300
2186 1020: -6.7758
2187 1021: 3.0929
2188 1022: 3.7940
2189 1023: 3.0392
2190 1024: 0.2890
2191 ... Figure 441, plot (a)
2192 |#
2193
2194 (multiple-value-bind (the-periodogram freqs N-f)
2195 (periodogram
2196 *forward-pred-errors*
2197 :sampling-time 1/4)
2198 (format t "~&Figure 441, plot (b) ...")
2199 (dotimes (i 7)
2200 (format t "~&~6,4F: ~8,4F"
2201 (svref freqs i)
2202 (svref the-periodogram i)))
2203 (format t "~&...")
2204 (dotimes (i 7)
2205 (format t "~&~6,4F: ~8,4F"
2206 (svref freqs (+ i (- N-f 7)))
2207 (svref the-periodogram (+ i (- N-f 7)))))
2208 (format t "~&... Figure 441, plot (b)")
2209 (values))
2210
2211 #|
2212 ;;; Here is what is printed when the above form is evaluated:
2213 Figure 441, plot (b) ...
2214 0.0039: 12.3693
2215 0.0078: 4.4997
2216 0.0117: 7.5193
2217 0.0156: 8.1898
2218 0.0195: 14.8114
2219 0.0234: 19.1417
2220 0.0273: 18.2599
2221 ...
2222 1.9766: 17.3194
2223 1.9805: 17.9742
2224 1.9844: 18.6913
2225 1.9883: 10.7566
2226 1.9922: 25.4055
2227 1.9961: 10.0394
2228 2.0000: 6.1584
2229 ... Figure 441, plot (b)
2230 |#
2231
2232 (multiple-value-bind (Parzen-m-55 freqs N-f nu B_W)
2233 (lag-window-spectral-estimate
2234 (acvs *forward-pred-errors*)
2235 #'(lambda (lag)
2236 (parzen-lag-window lag 55))
2237 :sampling-time 1/4
2238 :sdf-transformation nil)
2239 (let ((postcolored-Parzen-m-55 (postcolor-spectral-estimate
2240 Parzen-m-55
2241 (ar-coeffs->prewhitening-filter *ar-5-coeffs*)
2242 (length *forward-pred-errors*))))
2243 (format t "~&Figure 441, thick curve, plot (c) ...")
2244 (dotimes (i 7)
2245 (format t "~&~6,4F: ~8,4F"
2246 (svref freqs i)
2247 (svref postcolored-Parzen-m-55 i)))
2248 (format t "~&...")
2249 (dotimes (i 7)
2250 (format t "~&~6,4F: ~8,4F"
2251 (svref freqs (+ i (- N-f 7)))
2252 (svref postcolored-Parzen-m-55 (+ i (- N-f 7)))))
2253 (format t "~&... Figure 441, thick curve, plot (c)")
2254 (format t "~&equivalent degrees of freedom = ~5,1F" nu)
2255 (format t "~&smoothing window bandwidth = ~6,4F" B_W))
2256 (values))
2257
2258 #|
2259 ;;; Here is what is printed when the above form is evaluated:
2260 Figure 441, thick curve, plot (c) ...
2261 0.0039: 51.6943
2262 0.0078: 51.7057
2263 0.0117: 51.7240
2264 0.0156: 51.7487
2265 0.0195: 51.7789
2266 0.0234: 51.8134
2267 0.0273: 51.8512
2268 ...
2269 1.9766: -5.6470
2270 1.9805: -5.5772
2271 1.9844: -5.5188
2272 1.9883: -5.4726
2273 1.9922: -5.4392
2274 1.9961: -5.4190
2275 2.0000: -5.4122
2276 ... Figure 441, thick curve, plot (c)
2277 equivalent degrees of freedom = 68.7
2278 smoothing window bandwidth = 0.1349
2279 |#

  ViewVC Help
Powered by ViewVC 1.1.5