Thread: [Matlisp-users] FFT/IFFT Seriously Broke for Matrices

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 422-6466

All,

1.  I noted the following today.  FFT & IFFT only works for vectors.
    2D matricies wont compute correctly.  Try (m- (ifft (fft x)) x)
    on x = (ones 16 2).

2.  The problem seems to be a misplaced closing paren on copying the
    input matrix, X, to a working matrix RESULT.

3.  I have made the following change to the #+CMU portions of the code
    (the diff's are attached)

    a.  Add two new functions, FFT! and IFFT!, and exported them.  These
	are useful themselves for performing fast convolution without
	the need for unneeded cons'ing.  I.e. (setf x (ifft (m.* (fft x) y)))
        can be replaced with (ifft! (m.*! y (fft! x)))

    b.  Removed the loops for scaling in IFFT (too much cons'ing via MATRIX-REF)
        Performed scaling with SCAL!...much faster

    c.  FFT and IFFT call FFT! and IFFT! after create copy of matrix.

4.  Can someone look into making similar mods to the allegro portions of the code?

5.  I recommend updating CVS with the changes soon.

6.  Files modified:  FFT.LISP, DFFTPACK.LISP, PACKAGES.LISP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		   Changes to PACKAGES.LISP follow
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
*** src/matlisp/packages.lisp	Tue Oct  8 19:30:39 2002
--- src/matlisp-working/packages.lisp	Sat Feb  1 08:02:37 2003
***************
*** 207,213 ****
       "EIG"
       "EYE"
       "FFT"
!      "FFT"
       "FILL-MATRIX"
       "FLOAT-MATRIX"
       "FLOAT-MATRIX-ARRAY-TYPE"
--- 207,213 ----
       "EIG"
       "EYE"
       "FFT"
!      "FFT!"
       "FILL-MATRIX"
       "FLOAT-MATRIX"
       "FLOAT-MATRIX-ARRAY-TYPE"
***************
*** 224,229 ****
--- 224,230 ----
       "GETRS!"
       "HELP"
       "IFFT"
+      "IFFT!"
       "IMAG"
       "JOIN"
       "LOAD-BLAS-&-LAPACK-BINARIES"
***************
*** 285,290 ****
--- 286,292 ----
       "NUMBER-OF-ROWS"
       "ONES"
       "PRINT-ELEMENT"
+      "PINV"
       "QR"
       "QR!"
       "GEQR!"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		   Changes to DFFTPACK.LISP follow
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
*** src/matlisp/src/dfftpack.lisp	Tue Jul 11 14:02:03 2000
--- src/matlisp-working/src/dfftpack.lisp	Sat Feb  1 06:07:42 2003
***************
*** 126,132 ****
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :output)
    (wsave (* :double-float) :input))

--- 126,132 ----
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :input-output)
    (wsave (* :double-float) :input))

***************
*** 174,178 ****
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :output)
    (wsave (* :double-float) :input))
--- 174,178 ----
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :input-output)
    (wsave (* :double-float) :input))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		   Changes to FFT.LISP follow
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
*** src/matlisp/src/fft.lisp	Thu Oct 25 17:48:45 2001
--- src/matlisp-working/src/fft.lisp	Sat Feb  1 07:39:51 2003
***************
*** 217,238 ****
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
  	  (zfftf n (store result) wsave))

!       (dotimes (j (ncols x))
! 	(declare (type fixnum j))
! 	 (dotimes (i (nrows x))
! 	   (declare (type fixnum i))
! 	   (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	 (with-vector-data-addresses ((addr-result (store result))
! 				      (addr-wsave wsave))
! 	    (incf-sap :complex-double-float addr-result (* j n))
! 	    (dfftpack::fortran-zfftf n addr-result addr-wsave))))

!       result))

  #+:allegro
--- 217,241 ----
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
+ 
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
  	  (zfftf n (store result) wsave))

!       ;; Do FFT by column...first we copy the array
!       (progn
! 	;; copy X to a working array, RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 	  
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x))) 
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))

! 	;; Then do FFT by column
! 	(fft! result)))))

  #+:allegro
***************
*** 286,313 ****
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (let ((scale-factor (/ (float n 1d0))))
! 	    (dotimes (k n)
! 	      (setf (matrix-ref result k)
! 		    (* scale-factor (matrix-ref result k))))))
! 
! 	(let ((scale-factor (/ (float n 1d0))))
! 	  (dotimes (j (ncols x))
! 	    (declare (type fixnum j))
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	    (with-vector-data-addresses ((addr-result (store result))
! 					 (addr-wsave wsave))
! 	      (incf-sap :complex-double-float addr-result (* j n))
! 	      (dfftpack::fortran-zfftb n addr-result addr-wsave))
! 	    ;; Scale the result
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (* scale-factor (matrix-ref x i j)))))))

!     result))

  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
--- 289,309 ----
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (scal! (/ (float n 1d0)) result))

!       ;; Perform the IFFT by columns...first we copy the array
!       (progn
! 	;; Copy X to the working Array RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x)))
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))

+ 	;; The we do the IFFT
+ 	(ifft! result)))))

  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
***************
*** 348,350 ****
--- 344,403 ----

        result))

+ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+ ;; Modification for destructive FFT's.  The purpose is to eliminate
+ ;; cons'ing for certain application requiring numerous calls to these
+ ;; routines.
+ 
+ (defgeneric fft! (x)
+   (:documentation "See FFT but note that the optional N is NOT permitted.
+   Performs in place FFT modifying X.  This will only work for a complex matrix."))
+ 
+ (defgeneric ifft! (x)
+   (:documentation "See IFFT but note that the optional N is NOT permitted.
+   Performs in place IFFT modifying X.  This will only work for a complex matrix."))
+ 
+ #+:cmu  
+ (defmethod fft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		(max (nrows x) (ncols x))
+ 	      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(zfftf n (store x) wsave)
+ 
+       (dotimes (j (ncols x))
+ 	(declare (type fixnum j))
+ 	(with-vector-data-addresses ((addr-x (store x))
+ 				     (addr-wsave wsave))
+ 	    (incf-sap :complex-double-float addr-x (* j n))
+ 	    (dfftpack::fortran-zfftf n addr-x addr-wsave)))))
+ 
+   ;; return x, the result
+   x)
+ 
+ #+:cmu
+ (defmethod ifft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		      (max (nrows x) (ncols x))
+ 		      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(progn
+ 	  (zfftb n (store x) wsave)
+ 	  (scal! (/ (float n 1d0)) x))
+ 
+       ;; IFFT by column
+       (dotimes (j (ncols x))
+ 	(declare (type fixnum j))
+ 	(with-vector-data-addresses ((addr-x (store x))
+ 				     (addr-wsave wsave))
+ 	      (incf-sap :complex-double-float addr-x (* j n))
+ 	      (dfftpack::fortran-zfftb n addr-x addr-wsave))))
+ 	    
+     ;; Scale the result
+     (scal! (/ (float n 1d0)) x))
+   x)
+ 

Thread: [Matlisp-users] FFT/IFFT Seriously Broke for Matrices

matlisp-users