diff --git a/documentation/index.html b/documentation/index.html index b57ce50e1c1ba160e3c8929a44b2caedcf7ab5ff..07b695ba3b32b088015e2b5a464f17e183c75d40 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 0000000000000000000000000000000000000000..d3ec0b4ca3c104a84aae5f9e984f287102d7eacf --- /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 b29da782ed0aa6d64640baf7ddab5d4edcba4884..b27f230fb225885ef3982eafa6db3ea16195494a 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)