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
