From 976eac7baeca4828bbcf5200f01d4b0637754d36 Mon Sep 17 00:00:00 2001
From: Sumant Oemrawsingh <soemraws@xs4all.nl>
Date: Fri, 27 Nov 2009 00:25:20 +0100
Subject: [PATCH] Added some extra functions for use with FFTs

fft-frequency-vector can be used to make a vector that contains the actual
sample frequencies (so you don't have to guess or write your own).

fft-shift and fft-inverse-shift shift the zero frequency to the center (for
easy plotting), and the inverse (for proper inverse FFT).

These functions are comparable to fftfreq, fftshift and ifftshift in numpy (if
anyone cares...).
---
 documentation/index.html            | 13 +++++
 fast-fourier-transforms/extras.lisp | 90 +++++++++++++++++++++++++++++
 gsll.asd                            |  1 +
 3 files changed, 104 insertions(+)
 create mode 100644 fast-fourier-transforms/extras.lisp

diff --git a/documentation/index.html b/documentation/index.html
index b57ce50e..07b695ba 100644
--- a/documentation/index.html
+++ b/documentation/index.html
@@ -407,6 +407,19 @@ GSL doesn't have them.
   <li><code>invert-matrix</code> finds the inverse of a matrix and uses
   GSL's LU decomposition functions.</li>
   <li>IEEE floating point number analysis.</li>
+  <li><code>fft-frequency-vector</code> returns a vector where the sample
+  frequencies are contained. If you perform an FFT on a vector of a given size
+  and :sample-size, this vector will contain the sample frequencies in order.
+  If the :shifted keyword is T, then the frequencies are ordered in ascending
+  order.</li>
+  <li><code>fft-shift</code> returns a copy of a vector where the zero
+  frequency has been shifted to the center; the frequency components will be
+  sorted according to their frequency, in ascending order. Optionally, a
+  :stride can be provided.</li>
+  <li><code>fft-inverse-shift</code> performs the inverse action of fft-shift;
+  the zero and positive frequency components are shifted to the beginning,
+  so that the resulting vector is suitable for an inverse FFT. Optionally, a
+  :stride can be provided.</li>
 </ul>
 
 <h2>Status</h2>
diff --git a/fast-fourier-transforms/extras.lisp b/fast-fourier-transforms/extras.lisp
new file mode 100644
index 00000000..d3ec0b4c
--- /dev/null
+++ b/fast-fourier-transforms/extras.lisp
@@ -0,0 +1,90 @@
+;; Extras for FFT.
+;; Sumant Oemrawsingh, Wed Nov 25 2009 - 22:03
+
+(in-package :gsll)
+
+(defun fft-frequency-step (size &optional (sample-spacing 1))
+  "Return the frequency step for the FFT of a vector with the given size and
+  assuming the given sample spacing."
+  (/ (* size sample-spacing)))
+
+(defun fft-highest-frequency (size &optional (sample-spacing 1))
+  "Return the highest frequency for the FFT of a vector with the given size and
+  assuming the given sample spacing."
+  (if (evenp size)
+    (values (/ sample-spacing 2) (/ size 2))
+    (- (/ sample-spacing 2) (/ (fft-frequency-step size sample-spacing) 2))))
+
+(defun fft-frequency-split (size)
+  "Return the index where the highest frequency component is located in a
+  vector of the given size, after applying an FFT.  The second value that is
+  returned, is the most negative frequency divided by the sample frequency
+  step."
+  (if (evenp size)
+    (let ((n (/ size 2)))
+      (values n (- 1 n)))
+    (let ((n (/ (- size 1) 2)))
+      (values n (- n)))))
+
+(export 'fft-frequency-vector)
+(defun fft-frequency-vector (element-type size &key (sample-spacing 1) shifted)
+  "Make and return a vector that contains the sample frequencies of an FFT
+  that has been applied to a vector with the given size and :sample-spacing.
+  
+  If :shifted is T, then the vector will contain the sample frequencies after
+  FFT-SHIFT has been applied to the result of an FFT."
+  (let* ((df (fft-frequency-step size sample-spacing))
+        pos-f
+        neg-f)
+    (multiple-value-bind (max-f-index min-f-index)
+      (fft-frequency-split size)
+      (setf pos-f
+            (loop for i from 0 to max-f-index
+                  collect (coerce (* i df) element-type))
+            neg-f
+            (loop for i from min-f-index below 0
+                  collect (coerce (* i df) element-type)))
+      (make-marray element-type
+                   :initial-contents (if shifted
+                                       (nconc neg-f pos-f)
+                                       (nconc pos-f neg-f))))))
+
+(export 'fft-shift)
+(defun fft-shift (vector &key (stride 1))
+  "Return a copy of a vector that is the result of an FFT, with the zero and
+  positive frequencies shifted to the center and end, so that the data is
+  suitable for e.g. plotting."
+  (let* ((size (size vector))
+         (n (floor size stride))
+         (shifted (copy vector))
+         (split (fft-frequency-split n))
+         (zero-pos (if (evenp n) (- split 1) split)))
+    ;; Positive frequencies
+    (loop for i from zero-pos to (- n 1)
+          do (setf (maref shifted (* i stride))
+                   (maref vector (* (- i zero-pos) stride))))
+    ;; Negative frequencies
+    (loop for i from (+ split 1) to (- n 1)
+          do (setf (maref shifted (* (- i split 1) stride))
+                   (maref vector (* i stride))))
+    shifted))
+
+(export 'fft-inverse-shift)
+(defun fft-inverse-shift (vector &key (stride 1))
+  "Return a copy of a vector where the zero and positive frequency components
+  are shifted to the beginning, i.e. ordered so that it is suitable for an
+  inverse FFT."
+  (let* ((size (size vector))
+         (n (floor size stride))
+         (inv-shifted (copy vector))
+         (split (fft-frequency-split n))
+         (zero-pos (if (evenp n) (- split 1) split)))
+    ;; Positive frequencies
+    (loop for i from zero-pos to (- n 1)
+          do (setf (maref inv-shifted (* (- i zero-pos) stride))
+                   (maref vector (* i stride))))
+    ;; Negative frequencies
+    (loop for i from (+ split 1) to (- n 1)
+          do (setf (maref inv-shifted (* i stride))
+                   (maref vector (* (- i split 1) stride))))
+    inv-shifted))
diff --git a/gsll.asd b/gsll.asd
index b29da782..b27f230f 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -143,6 +143,7 @@
              (:file "select-direction")
              (:file "unpack")
              (:file "discrete")
+             (:file "extras")
 	     (:file "example")))
    (:module random
 	    :depends-on (init data)
-- 
GitLab