fast fourier transforms 1D 2D 3D
Posted: Fri Dec 21, 2018 10:28 pm
I needed a set of FFT's for a project, 1D 2D 3D forward and inverse
note for 1D you can use either FFT or FFTCT if you want an in place transform
note for 1D you can use either FFT or FFTCT if you want an in place transform
Code: Select all
; /********************************************************************
; F A S T F O U R I E R T R A N S F O R M P R O G R A M S
;
; by Wang Jian-Sheng 4 Nov 1998, added fft2D(), 11 Apr 2003
; ---------------------------------------------------------------------
; port to PB 5.62 idle 22/12/18
;
; Reference: "Computational Frameworks for the Fast Fourier
; Transform", Charles Van Loan, SIAM, 1992.
;
; There are many FFT algorithms, the most important ones are
; COOLEY-TUKEY: in place, bit reversal
; STOCKHAM AUTOSORT: additional memory size of input Data
; MIXED RADIX: 20% less operations comparing To Cooley-Tukey
; PRIME FACTOR: arbitrary length n
;
; We use a combination of the Stockham autosort algorithm 1.7.2,
; page 57, And multirow Cooley-Tukey (3.1.7), page 124, of the
; reference above.
;
; The discrete Fourier transform is defined by
; y[k] = sum_(j=0,n-1) x[j] Exp(-2 Pi sqrt[-1] j k/n),
; k=0,1,...,n-1. The factor (1/n) is Not included.
; If y[]<-x[]; fft(x,n,1); fft(x,n,-1); then y[]==x[]/n is true.
; Three dimensional transform is generalized straightforwardly.
;
; Interface And usage:
; 1D Fourier transform
; Use: fft(x, n, flag)
; x : an Array of Structure type complex;
; n : the size of Data, must be a power of 2;
; flag : 1 For forward transform, -1 For inverse transform.
;
; 3D Fourier transform
; Use : fft3D(x, n1, n2, n3, flag)
; x : 1D Array of type complex representing 3D Array;
; mapping through C convention, i.e.,
; (i,j,k) -> k + n3*j + n2*n3*i;
; n1, n2, n3 : dimensions in three directions;
; flag : same As in 1D.
;
; 2D FFT is similar but With n1 And n2 only.
Structure complex
Re.d
Im.d
EndStructure
Structure arcomplex
ar.complex[0]
EndStructure
Procedure _stockham(*x.arcomplex,n.i,flag.i,n2.i,*y.arcomplex)
Protected *y_orig.arcomplex
Protected *tmp.complex
Protected i.i, j.i, k.i, k2.i, Ls.i, r.i, jrs.i;
Protected half, m, m2;
Protected wr.d, wi.d, tr.d, ti.d;
*y_orig = *y
half = n >> 1
r = half
Ls = 1
While(r >= n2)
*tmp = *x
*x = *y
*y = *tmp
m = 0
m2 = half
j=0
While j < ls
wr = Cos(#PI*j/Ls)
wi = -flag * Sin(#PI*j/Ls)
jrs = j*(r+r)
k = jrs
While k < jrs+r
k2 = k + r
tr = wr * *y\ar[k2]\Re - wi * *y\ar[k2]\Im
ti = wr * *y\ar[k2]\Im + wi * *y\ar[k2]\Re
*x\ar[m]\Re = *y\ar[k]\Re + tr
*x\ar[m]\Im = *y\ar[k]\Im + ti
*x\ar[m2]\Re = *y\ar[k]\Re - tr
*x\ar[m2]\Im = *y\ar[k]\Im - ti
m+1
m2+1
k+1
Wend
j+1
Wend
r >> 1
Ls << 1
Wend
If (*y <> *y_orig)
For i = 0 To n -1
*y\ar[i] = *x\ar[i]
Next
EndIf
EndProcedure
Procedure _cooley_tukey(*x.arcomplex,n.i,flag.i,n2.i)
Protected c.complex
Protected i.i, j.i, k.i, m.i, p.i, n1.i
Protected Ls.i, ks.i, ms.i, jm.i, dk.i
Protected wr.d, wi.d, tr.d, ti.d
Protected tm.l
n1 = n/n2
k=0
While k < n1
j = 0
m = k
p = 1
While p < n1
j << 1
j + (m & 1)
m >> 1
p << 1
Wend
If j > k
i=0
While i < n2
c\Re = *x\ar[k*n2+i]\Re
c\im = *x\ar[k*n2+i]\Im
*x\ar[k*n2+i]\Re = *x\ar[j*n2+i]\Re
*x\ar[k*n2+i]\im = *x\ar[j*n2+i]\im
*x\ar[j*n2+i]\Re = c\Re
*x\ar[j*n2+i]\im = c\Im
i+1
Wend
EndIf
k+1
Wend
p = 1
While p < n1
Ls = p
p << 1
jm = 0
dk = p*n2
j=0
While j < Ls
wr = Cos((#PI*j/Ls))
wi = -flag * Sin((#PI*j/Ls))
k=jm
While k < n
ks = k + Ls*n2
i=0
While i < n2
m = k + i
ms = ks + i
tr = (wr * *x\ar[ms]\Re) - (wi * *x\ar[ms]\Im)
ti = (wr * *x\ar[ms]\Im) + (wi * *x\ar[ms]\Re)
*x\ar[ms]\Re = *x\ar[m]\Re - tr
*x\ar[ms]\Im = *x\ar[m]\Im - ti
*x\ar[m]\Re + tr
*x\ar[m]\Im + ti
i+1
Wend
k+dk
Wend
jm + n2
j+1
Wend
Wend
EndProcedure
Procedure fft(*x.arcomplex,n.i,flag.i=1)
Protected *y.arcomplex
*y = AllocateMemory((n)*SizeOf(complex))
_stockham(*x, n, flag, 1, *y)
FreeMemory(*y)
EndProcedure
Procedure fft2D(*x.arcomplex,n1.i,n2.i,flag.i=1)
Protected *y.arcomplex;
Protected i, n
n = n1*n2
*y = AllocateMemory(n2*SizeOf(complex))
i=0
While i < n
_stockham(@*x\ar[i],n2,flag, 1,*y)
i+n2
Wend
FreeMemory(*y)
_cooley_tukey(@*x\ar[0], n, flag, n2)
EndProcedure
Procedure fftCT(*x.arcomplex,n.l,flag.l=1)
_cooley_tukey(@*x\ar[0],n,flag,1)
EndProcedure
Procedure fft3D(*x.arcomplex,n1.i,n2.i,n3.i,flag.i=1)
Protected *y.arcomplex;
Protected i, n, n23;
n23 = n2*n3;
n = n1*n23;
*y = AllocateMemory(n23*SizeOf(complex))
For i=0 To n-1
_stockham(@*x\ar[i], n3, flag, 1, *y)
i+n3
Next
For i=0 To n-1
_stockham(@*x\ar[i], n23, flag, n3, *y)
i+n23
Next
FreeMemory(*y)
_cooley_tukey(@*x\ar[0], n, flag, n23)
EndProcedure
CompilerIf #PB_Compiler_IsMainFile
Global Dim inp.complex(7)
Global a,b
For a = 0 To 3
inp(a)\Re = 1
Next
For a = 0 To 7
Debug StrF(inp(a)\Re,3) + " " + StrF(inp(a)\Im,3)
Next
fft(@inp(0),8)
Debug "========================================================="
For a = 0 To 7
Debug StrF(inp(a)\Re,3) + " " + StrF(inp(a)\Im,3)
Next
Debug "========================================================="
fft(@inp(0),8,-1)
For a = 0 To 7
Debug StrF(inp(a)\Re /8 ,3) + " " + StrF(inp(a)\Im,3)
Next
Debug "2D Input======================================================="
Global Dim inp2.complex(3,3)
For a = 1 To 2
For b = 1 To 2
inp2(a,b)\Re = 3
Debug inp2(a,b)\re
Next
Next
Debug "2D FFT========================================================="
fft2d(@inp2(0,0),4,4)
For a = 1 To 2
For b = 1 To 2
Debug StrF(inp2(a,b)\Re,3) + " " + StrF(inp2(a,b)\im,3)
Next
Next
Debug "2D IFFT========================================================="
fft2d(@inp2(0,0),4,4,-1)
For a = 1 To 2
For b = 1 To 2
Debug inp2(a,b)\re / 16
Next
Next
CompilerEndIf
