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

[<ma...@ll...>]

Okay...so I embarrass myself again...I didn't thoroughly check out the vector
case, since I was concertrating on the 2D case.  The requires another modification 
to fft.lisp.  The below DIFF against fft.lisp should be the right one.

sorry

mike



*** src/matlisp/src/fft.lisp	Thu Oct 25 17:48:45 2001
--- src/matlisp-working/src/fft.lisp	Sun Feb  2 12:22:45 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,242 ----
  		       ((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)
! 	  result)
  
!       ;; 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)
--- 290,310 ----
  	(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 ****
--- 345,406 ----
  
        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)
+ 	(progn
+ 	  (zfftf n (store x) wsave)
+ 	  x)
+ 
+       (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)
+ 	  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)
+