example.lisp 6.27 KB
Newer Older
1
2
;; Example FFT: transform a pulse (using the "clean" fft interface)
;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24
3
;; Time-stamp: <2009-11-16 23:10:53EST example.lisp>
4
5
6
7
8
9

(in-package :gsl)

;;;;****************************************************************************
;;;; Pulse test
;;;;****************************************************************************
10
11
12
13
14
15
16
17
18

;;; Here is an example program modelled after the example given in Section
;;; 15.3 of the GSL Manual, which computes the FFT of a short pulse. To make
;;; the resulting fourier transform real the pulse is defined for equal
;;; positive and negative times (-10 ... 10), where the negative times wrap
;;; around the end of the array.
;;; 
;;; The output array from the example in Section 15.3 of the GSL Manual can be
;;; reproduced with:
19
;;; (fft-pulse-test '(complex double-float) 128)
20
21
22
;;;
;;; This example program also yields the same output array as the example
;;; program in Section 15.4 of the GSL Manual:
23
;;; (fft-pulse-test '(complex double-float) 630)
24

25
(defun fft-pulse-test (element-type dimension)
26
27
28
29
30
31
32
33
34
35
  (assert (and (integerp dimension) (> dimension 20)))
  (let ((pulse (make-marray element-type :dimensions dimension))
        (init-value (coerce 1 element-type)))
    (setf (maref pulse 0) init-value)
    (loop for i from 1 to 10
          do (setf (maref pulse i) init-value
                   (maref pulse (- dimension i)) init-value))
    (forward-fourier-transform pulse)))

(save-test
36
  fft-pulse
37
38
39
40
41
42
43
44
  (fft-pulse-test '(complex single-float) 128)
  (fft-pulse-test '(complex single-float) 630)
  (fft-pulse-test '(complex double-float) 128)
  (fft-pulse-test '(complex double-float) 630)
  (fft-pulse-test 'single-float 128)
  (fft-pulse-test 'single-float 630)
  (fft-pulse-test 'double-float 128)
  (fft-pulse-test 'double-float 630))
Liam Healy's avatar
Liam Healy committed
45

46
47
48
49
;;;;****************************************************************************
;;;; Random vector transformations, from the GSL tests
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
50
51
52
53
54
55
56
57
58
59
60
;; From gsl-1.11/fft/urand.c
(let ((urand-seed 1))
  (defun urand ()
    "Generate a random number.  See fft/urand.c."
    (setf urand-seed (logand #x7fffffff (+ 12345 (* urand-seed 1103515245))))
    (/ urand-seed 2147483648.d0))
  (defun reset-urand ()
    (setf urand-seed 1)
    (values)))

;; (make-urand-vector '(complex double-float) 5)
61
(defun make-urand-vector (element-type dimension &key (stride 1))
Liam Healy's avatar
Liam Healy committed
62
  "Make a vector with random elements."
63
64
  (let ((vec (make-marray `(complex ,(component-float-type element-type))
			  :dimensions (list (* stride dimension)))))
65
    (loop for i from 0 below (* stride dimension) by stride
Liam Healy's avatar
Liam Healy committed
66
       do
67
68
69
70
       (setf (maref vec i)
	     (if (subtypep element-type 'complex)
		 (coerce (complex (urand) (urand)) element-type)
		 (complex (coerce (urand) element-type)))))
Liam Healy's avatar
Liam Healy committed
71
72
    vec))

73
74
75
76
77
78
79
80
81
;; (reset-urand)
;; (fft-signal-real-noise '(complex double-float) 10)
(defun fft-signal-real-noise (element-type dimension &key (stride 1))
  (let ((random-vector
          (make-urand-vector element-type dimension :stride stride)))
    (values
      random-vector
      (forward-discrete-fourier-transform random-vector :stride stride))))

82
83
84
85
86
87
88
89
90
91
(defun realpart-vector (complex-vector)
  "The real vector consisting of the real part of the complex vector."
  (let ((real-vector
	 (make-marray
	  (component-float-type (element-type complex-vector))
	  :dimensions (dimensions complex-vector))))
    (loop for i below (total-size complex-vector) do
	 (setf (maref real-vector i)
	       (realpart (maref complex-vector i))))
    real-vector))
Liam Healy's avatar
Liam Healy committed
92

93
94
95
96
97
(defun size-vector-real (vector &key (stride 1))
  "Return the size of a vector while taking the stride into account."
  (coerce (floor (size vector) stride)
          (component-float-type (element-type vector))))

98
99
100
101
102
103
104
105
106
107
108
109
(defun forward-fft-rc (vector complexp &key (stride 1))
  "Return the forward FFT of the vector."
  (let ((forward
	 (forward-fourier-transform
	  (if complexp
	      (copy vector)
	      (realpart-vector vector))
	  :stride stride)))
    (if complexp
	forward
	(unpack forward :unpack-type 'complex :stride stride))))

110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
;(defun test-fft-noise (element-type size &key (stride 1))
;  "Test for FFT with random data (noise); returns the DFT answer and the computed FFT answer.
;   See test_real_radix2 etc. in fft/test.mc."
;  (let ((random-vector (make-urand-vector element-type size :stride stride)))
;    (values
;     (forward-fft-rc
;      random-vector (subtypep element-type 'complex) :stride stride)
;     (forward-discrete-fourier-transform random-vector :stride stride))))

(defun test-fft-noise-r (vector &key (stride 1))
  "Test forward and inverse FFT for a real vector, and return both results in unpacked form."
  (let* ((forward
           (forward-fourier-transform (realpart-vector vector) :stride stride))
         (inverse
           (forward-fourier-transform (copy forward) :half-complex t :stride stride)))
    (values (unpack forward :unpack-type 'complex :stride stride)
            (unpack (elt/ inverse (size-vector-real inverse :stride stride)) :stride stride))))
127

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
(defun test-fft-noise-c (vector &key (stride 1))
  "Test forward, inverse and backward FFT for a complex vector and return all three results."
  (let* ((forward
           (forward-fourier-transform (copy vector) :stride stride))
         (inverse
           (inverse-fourier-transform (copy forward) :stride stride))
         (backward
           (backward-fourier-transform (copy forward) :stride stride)))
    (values forward
            inverse
            backward)))

(defun test-fft-noise (element-type size &key (stride 1))
  (multiple-value-bind
      (c-data fft-c-data)
      (fft-signal-real-noise element-type size :stride stride)
    (if (subtypep element-type 'complex)
      (multiple-value-bind
        (forward inverse backward)
        (test-fft-noise-c c-data :stride stride)
        (values
          fft-c-data
          forward
          c-data
          inverse
          backward ;; just returning this one to see that it is still correct
          (elt/ (copy backward) (size-vector-real backward :stride stride))))
      (multiple-value-bind
        (forward inverse)
        (test-fft-noise-r c-data :stride stride)
        (values fft-c-data forward c-data inverse)))))
159

160
161
162
;; (require :gsll)
;; (test-fft-noise 'double-float 10 :stride 1)
;; (test-fft-noise '(complex double-float) 10 :stride 1)