I may need to start a separate topic to cover everything that I'm sending
in this post .
Okay , continuing with the theme of 3d data points , I have generated a
3d array from a 3d fft of a 3d rectangular object.
Then I exchanged data values , within this array , to produce a centered
spectrum.
The selection levels were chosen to filter out the smaller magnitude values.
This I then graphed using the 3d functions developed previously , the whole
selection was rotated to highlight the features.
The 3d fft hasn't been thoroughly check , however the results appear to
be correct.
Any comment regards this is appreciated.
All of this code to be placed in the same directory.
There should be 6 portions of code .
All code :
(c) Copyright 2016 Luxan ,
sciwise@ihug.co.nz
1
Code: Select all
'
' --------------------------------------------------------------------
'
' fn3d5ald.bi , code as an include file .
'
' (c) Copyright 2016 sciwise@ihug.co.nz
'
' Luxan
'
' Edward.Q.Montague [ alias ]
'
'
' 3d function
'
'
' Screen pages are used .
'
'
' -------------------------------------------------------------------
'
type point
x as single
y as single
z as single
u as single ' possible extension for special coord system
end type
'
' -------------------------------------------------------------------
'
declare function roty(q as point,angy as single) as point
declare function rotx(q as point,angx as single) as point
declare function rotz(q as point,angz as single) as point
declare function persp(q as point,d as single) as point
declare function pally(pal() as long) as integer
declare function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
declare function f3(x as single,y as single,z as single) as single
declare function Generate(b4() as single,n as integer,m as integer,p as integer) as integer
declare function normb(b4() as single , n as integer , m as integer , p as integer) as integer
declare function logb(b4() as single,n as integer,m as integer,p as integer) as integer
declare function efilter(b4() as single , b4c() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
declare function dproj(b4c() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' =====================================
'
const Pi = 4*atn(1)
'dim as integer fg
'
'n = 8
'm = 8
' p = 8
'
redim as single b4(0 to 4,0 to 4 , 0 to 4)
' ----------------------------------------------------------
'
dim pal(15) as long
'
' -----------------------------------------------------------
'
dim as integer np,ne
dim as point p1a()
dim as integer edge()
dim as single l1,l2 ' levels
'
' ----------------------------------------------------------
'
2
Code: Select all
'
' ---------------------------------------------------
'
' fn3d5l.bas
'
'
' Load file for 3d graph
'
'
' ---------------------------------------------------
'
' number of vertices
'
storeA:
data 8
'
' vertex data , easier to keep track of
' data when we use multiple data statements.
'
store1:
data 1,1,1
data -1,1,1
data -1,-1,1
data 1,-1,1
data 1,1,-1
data -1,1,-1
data -1,-1,-1
data 1,-1,-1
'
' Number of edges.
'
storeB:
data 12
'
' edge data
'
store2:
data 1,2
data 1,4
data 1,5
data 2,3
data 2,6
data 3,4
data 3,7
data 4,8
data 5,6
data 5,8
data 6,7
data 7,8
'
' Number of edges.
'
storeC:
data 9
'
' edge data , back faces .
'
store3:
data 1 , 2
data 1 , 4
data 2 , 3
data 3 , 4
data 2 , 6
data 3 , 7
data 6 , 7
data 7 , 8
data 4 , 8
'
' -------------------------------------------------------------------------------
'
function rotx(q as point,angx as single) as point
'
' Rotate around x axis .
'
static as point p
'
p.x = q.x
p.y= q.y*cos(angx)-sin(angx)*q.z
p.z= q.z*cos(angx)+sin(angx)*q.y
'
return p
'
end function
'
' -----------------------------------------------------------------------------
'
function roty(q as point,angy as single) as point
'
' Rotate , a single point , around y axis .
'
static as point p
'
p.x = sin(angy)*q.z + cos(angy)*q.x
p.y = q.y
p.z = cos(angy)*q.z -sin(angy)*q.x
'
return p
'
end function
'
' -----------------------------------------------------------------------------
'
function rotz(q as point,angz as single) as point
'
' Rotate around z axis .
'
static as point p
'
p.x = sin(angz)*q.y + cos(angz)*q.x
p.y = cos(angz)*q.y-sin(angz)*q.x
p.z = q.z
'
return p
'
end function
'
' -----------------------------------------------------------------------------
'
function persp(q as point,d as single) as point
'
' 3d perspective .
'
' Applied to a single point .
'
' Add 2 to the numerator when using any negative z value.
'
' In this instance -1 <= z <= 1 , unit cube .
'
' Therefore 2 is appropriate .
'
static as point p
'
p.x = d*q.x/(q.z+3.5)
p.y = d*q.y/(q.z+3.5)
p.z = d
'
return p
'
end function
'
' ---------------------------------------------------------------
'
function pally(pal() as long) as integer
'
' Assign palette information .
'
'
' re defined 15 colour mode .
'
pal(0)=5
pal(1)=1
pal(2)=9
pal(3)=11
pal(4)=2
pal(5)=10
pal(6)=14
pal(7)=6
pal(8)=4
pal(9)=12
pal(10)=13
pal(11)=15
pal(12)=0
pal(13)=0
pal(14)=0
pal(15)=0
'
'
return 0
'
end function
'
' --------------------------------------------------------------
'
function cube( p1a() as point, ByRef np as integer, edge() as integer,ByRef ne as integer) as integer
'
' Read Cube coordinate points and edge sequence data .
'
static as integer i
'
restore storeA
read np
'
redim p1a(1 to np)
'
restore store1
for i =1 to np
read p1a(i).x
read p1a(i).y
read p1a(i).z
next i
'
'
restore storeC
read ne
'
redim edge(1 to ne,0 to 1)
'
restore store3
for i = 1 to ne
read edge(i,0)
read edge(i,1)
next i
'
return 0
'
end function
'
' ----------------------------------------------------------------
'
function dproj(b4c() as point , count as integer , p1a() as point ,np as integer,edge() as integer,ne as integer, pal() as long , ng as integer , div as integer) as integer
'
' -----------------------------------------------------------
'
' Draw cube frame and
'
' count , number of points in array b3c()
'
'
' p1a() point array for cube
' np number of points for cube array , this is usually a constant value .
' edge() edge array for cube
' ne number of edge pairs for edge() array , this is usually a constant value.
'
' pal() , palette array , this is a constant length .
'
' ng , the number of cycles for the complete 360 deg rotation.
'
' div , sets the size of each discrete rotation ; 2*Pi/div .
'
'
' Generate values from 3d function for each rotation .
'
' -----------------------------------------------------------
'
static as single c1,d1,c2,d2,theta
dim as single x,y,z,u,maxa,thi
dim as single x1,y1,z1,u1,x2,y2,z2,u2
static as integer i,j,k,ni,jt,i1,j1,k1
static as point p1,p4,p2(1 to np)
'
c1=1.2
d1=0.4
c2=1.35
d2=1.3
'
theta = Pi/div
'
for ni=0 to ng
'
for jt = 0 to div
cls
'
' .... Draw color palette ,,,,,,,
'
for i = 0 to 11
y1=(d2-d1)*i/15+d1
y2=(d2-d1)*(i+1)/15+d1
line(c1,y1)-(c2,y2),pal(i),bf
line(c1,y1)-(c2,y2),11,b
next i
'
' .......... Draw outer cube .............
'
thi = jt*theta
'
for k = 1 to np
p2(k) = roty(p1a(k),thi)
p2(k) = rotx(p2(k),-0.75)
p2(k)=persp(p2(k),2)
next k
'
for i1 = 1 to ne
j1 = edge(i1,0)
k1 = edge(i1,1)
x1 = p2(j1).x
y1 = p2(j1).y
' z1 = p2(j1).z
x2 = p2(k1).x
y2 = p2(k1).y
' z2 = p2(k1).z
line(x1,y1)-(x2,y2),14
next i1
'
' .................... 3d point values from array b4c() .................
'
for k=0 to count
p1 = b4c(k)
u = p1.u
p4 = roty(p1,thi)
p4 = rotx(p4,-0.75)
p4 = persp(p4,2)
pset(p4.x,p4.y),pal(u*11)
next k
'
ScreenCopy 1, 0
'
sleep 100
next jt
'
' repeat cycle (ng - ni) imes
'
next ni
'
return 0
'
end function
'
' ----------------------------------------------------------------
'
function Generate(b4() as single,n as integer,m as integer,p as integer) as integer
'f
' Generate data for array b4() from function f3 .
'
static as integer i,j,k
static as single x,y,z,u
'
redim as single b4(0 to n-1,0 to m-1 , 0 to p-1)
'
'
'
' .................... Generate 3d point values .................
'
for k = 0 to p-1
z = -1+2*k/(p-1)
for j = 0 to m-1
y = -1 + 2*j/(m-1)
for i = 0 to n-1
x = -1 + 2*i/(n-1)
u = f3(x,y,z)
u = abs(u)
b4(i,j,k) = u
next i
next j
next k
'
return 0
'
end function
'
' ----------------------------------------------------------------
'
function efilter(b4() as single , b4c() as point , n as integer , m as integer , p as integer,l1 as single,l2 as single) as integer
'
' Filter data from array b4() and selectively store in array b4c()
'
'
static as integer i,j,k,count,cn
static as single x,y,z,u
static as point p1
'
'
u=l1
if l1 > l2 then
l1=l2
l2=u
end if
'
'
' .................... Select 3d point values .................
'
count = 0
for k = 0 to p-1
z = -1+2*k/(p-1)
for j = 0 to m-1
y = -1 + 2*j/(m-1)
for i = 0 to n-1
x = -1 + 2*i/(n-1)
u = b4(i,j,k)
'
' Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then
count = count + 1
end if
'
next i
next j
next k
'
if (count > 0 ) then
'
redim as point b4c(count)
'
cn = 0
for k = 0 to p-1
z = -1+2*k/(p-1)
for j = 0 to m-1
y = -1 + 2*j/(m-1)
for i = 0 to n-1
x = -1 + 2*i/(n-1)
u = b4(i,j,k)
'
' Set level where point becomes visible .
'
if (u >= l1 and u<=l2) then
p1.x = x
p1.y = y
p1.z = z
p1.u = u
b4c(cn) = p1
cn = cn + 1
end if
'
next i
next j
next k
'
end if
'
return count
'
end function
'
' ----------------------------------------------------------------
'
function f3(x as single,y as single,z as single) as single
'
' 3d sinc function
'
'
static as single u,v,w
'
u=1
if (abs(x)>0) then u= sin(4*x*Pi)/(4*x*Pi)
v=1
if (abs(y)>0) then v= sin(2*y*Pi)/(2*y*Pi)
w=1
if (abs(z)>0) then w= sin(3*z*Pi)/(3*z*Pi)
'
u=u*v*w
u=log(abs(u)*12700000 + 1)/log(12700001)
return u
'
end function
'
' -----------------------------------------------------------
'
function normb(b4() as single , n as integer , m as integer , p as integer) as integer
'
' Normalize all data values .
'
'
static as single max,fp
static as integer i,j,k
'
max=0
for k=0 to p-1
for j=0 to m-1
for i=0 to n-1
fp=b4(i,j,k)
fp=abs(fp)
if ( fp > max ) then max = fp
next i
next j
next k
'
if (max =0 ) then max = 1
'
for k=0 to p-1
for j=0 to m-1
for i=0 to n-1
fp=b4(i,j,k)
fp=abs(fp)/max
b4(i,j,k) = fp
next i
next j
next k
'
return 0
'
end function
'
' ----------------------------------------------------
'
function logb(b4() as single,n as integer,m as integer,p as integer) as integer
'
' Natural Logarithm of all data in array b4()
'
'
static as integer i,j,k
static as single u
'
for k=0 to p-1
for j=0 to m-1
for i=0 to n-1
u = b4(i,j,k)
u=log(abs(u)*12700000 + 1)/log(12700001)
b4(i,j,k)=u
next i
next j
next k
'
return 0
'
end function
3
Code: Select all
'
' -------------------------------------------------------------------
'
' 3d graph of centered spectrum from 3d fft .
'
' -------------------------------------------------------------------
'
#include "fourier.bas"
#include "fourtst123.bi"
#include "fn3d5ald.bi"
declare function uxchg(bg3() as single,n as integer,m as integer,p as integer) as integer
'
' ---------------------------------------------------
'
fg = pally(pal() )
fg = cube(p1a() ,np ,edge() ,ne )
'
' -----------------------------------------------------------
'
n=8
m=8
p=8
'fg= ft1test(n )
'fg = ft2test(n , m )
'fg= ft3test(n ,m, p )
redim as single b4(0 to n-1,0 to m-1,0 to p-1)
redim as point b4c(4)
'fg=ft3gendat(b4() ,n ,m , p )
'fg=normb(b4(),n,m,p)
'fg=ft3Prtdat(b4() , n ,m , p )
'sleep
'end
screen 12,8,2
' image for working page #1 (visible page #0)
ScreenSet 1, 0
Cls
window (-1.5,-1.5)-(1.5,1.5)
line (-1.4,-1.4)-(1.4,1.4),11,b
ScreenCopy 1, 0
'
'redim b4c(0 to 4) as point
'
' --------------------------------------------------------------
'
n=128
m=128
p=128
redim as single b4(0 to n-1,0 to m-1,0 to p-1)
'fg = Generate(b4() ,n ,m ,p )
fg=ft3gendat(b4() ,n ,m , p )
fg=normb(b4(),n,m,p)
fg=logb(b4() ,n ,m ,p )
' l1=0.85
' l2=1
'
l1=0.78
l2=1
'
fg=uxchg(b4() , n ,m ,p )
fg= efilter(b4() , b4c() , n , m , p ,l1 ,l2 )
fg= dproj(b4c() , fg , p1a() ,np ,edge() ,ne , pal() , 4 , 32 )
'
redim as single b4(0 to 4,0 to 4 , 0 to 4)
redim b4c(4) as point
'
sleep
end
'
' ================================
'
#include "fourtst123x.bas"
#include "fn3d5l.bas"
function uxchg(bg3() as single , n as integer,m as integer,p as integer) as integer
'
'
' Exchange bg3() values into bg3()
'
'
static as integer i,j,k,i1,j1
'
' [k,[i,j]]
'
redim b1(0 to m-1) as single
redim xb(0 to m-1) as single
for k=0 to p-1
for i=0 to n-1
for j=0 to m-1
b1(j)=bg3(i,j,k)
next j
for j=0 to (m/2)-1
i1=j+(m/2)
xb(j)=b1(i1)
xb(i1)=b1(j)
next j
for j=0 to m-1
bg3(i,j,k)=xb(j)
next j
next i
next k
'exit function
'
' [k,[j,i]]
redim b1(0 to n-1) as single
redim xb(0 to n-1) as single
for k=0 to p-1
for j=0 to m-1
for i=0 to n-1
b1(i)=bg3(i,j,k)
next i
'
for i=0 to (n/2)-1
j1=i+(n/2)
xb(i)=b1(j1)
xb(j1)=b1(i)
next i
for i=0 to n-1
bg3(i,j,k)=xb(i)
next i
next j
next k
' exit function
'
' [j,[i,k]]
'
redim b1(0 to p-1) as single
redim xb(0 to p-1) as single
for j=0 to m-1
for i = 0 to n-1
'
for k = 0 to p-1
b1(k)=bg3(i,j,k)
next k
'
for k = 0 to (p/2)-1
j1=k+(p/2)
xb(k)=b1(j1)
xb(j1)=b1(k)
next k
for k = 0 to p-1
bg3(i,j,k)=xb(k)
next k
next i
next j
'
return 0
'
end function
'
' ---------------------------------------------------
'
4
Code: Select all
' ******************************************************************
' Fast Fourier Transform
' Modified from a Pascal program by Don Cross:
' http://groovit.disjunkt.com/analog/time-domain/fft.html
' ******************************************************************
' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------
#DEFINE FFTOk 0
' No error
#DEFINE FFTPowErr -1
' Number of samples is not a power of 2
#DEFINE FFTBoundErr -2
' Arrays do not start at index 0
' ------------------------------------------------------------------
' Mathematical constant
' ------------------------------------------------------------------
#DEFINE TwoPi 6.283185307179586#
' 2*Pi
' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------
TYPE Complex
X AS DOUBLE
Y AS DOUBLE
END TYPE
' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------
COMMON SHARED ErrCode AS INTEGER
' Error code from the last function evaluation
' ******************************************************************
#DEFINE MaxPower 32
FUNCTION IsPowerOfTwo(BYVAL X AS INTEGER) AS INTEGER
DIM AS INTEGER I, Y
Y = 2
FOR I = 1 TO MaxPower - 1
IF X = Y THEN RETURN -1
Y = Y SHL 1
NEXT I
RETURN 0
END FUNCTION
FUNCTION NumberOfBitsNeeded(BYVAL PowerOfTwo AS INTEGER) AS INTEGER
DIM AS INTEGER I
FOR I = 0 TO MaxPower
IF (PowerOfTwo AND (1 SHL I)) <> 0 THEN RETURN I
NEXT I
RETURN 0
END FUNCTION
FUNCTION ReverseBits(BYVAL Index AS INTEGER, BYVAL NumBits AS INTEGER) AS INTEGER
DIM AS INTEGER I, K, Rev
Rev = 0
I = Index
FOR K = 0 TO NumBits - 1
Rev = (Rev SHL 1) OR (I AND 1)
I = I SHR 1
NEXT K
RETURN Rev
END FUNCTION
SUB FourierTransform(BYVAL AngleNumerator AS DOUBLE, _
InArray() AS Complex, _
OutArray() AS Complex)
IF LBOUND(InArray) <> 0 OR LBOUND(OutArray) <> 0 THEN
ErrCode = FFTBoundErr
EXIT SUB
END IF
DIM AS INTEGER NumSamples
NumSamples = UBOUND(InArray) + 1
IF (NOT IsPowerOfTwo(NumSamples)) OR (NumSamples < 2) THEN
ErrCode = FFTPowErr
EXIT SUB
END IF
DIM AS INTEGER NumBits, I, J, K, N, BlockSize, BlockEnd
DIM AS DOUBLE Delta_angle, Delta_ar, Alpha, Beta, Tr, Ti, Ar, Ai
ErrCode = FFTOk
NumBits = NumberOfBitsNeeded(NumSamples)
FOR I = 0 TO NumSamples - 1
J = ReverseBits(I, NumBits)
OutArray(J).X = InArray(I).X
OutArray(J).Y = InArray(I).Y
NEXT I
BlockEnd = 1
BlockSize = 2
DO WHILE BlockSize <= NumSamples
Delta_angle = AngleNumerator / BlockSize
Alpha = SIN(0.5 * Delta_angle)
Alpha = 2 * Alpha * Alpha
Beta = SIN(Delta_angle)
I = 0
DO WHILE I < NumSamples
Ar = 1 ' cos(0)
Ai = 0 ' sin(0)
J = I
FOR N = 0 TO BlockEnd - 1
K = J + BlockEnd
Tr = Ar * OutArray(K).X - Ai * OutArray(K).Y
Ti = Ar * OutArray(K).Y + Ai * OutArray(K).X
OutArray(K).X = OutArray(J).X - Tr
OutArray(K).Y = OutArray(J).Y - Ti
OutArray(J).X = OutArray(J).X + Tr
OutArray(J).Y = OutArray(J).Y + Ti
Delta_ar = Alpha * Ar + Beta * Ai
Ai = Ai - (Alpha * Ai - Beta * Ar)
Ar = Ar - Delta_ar
J = J + 1
NEXT N
I = I + BlockSize
LOOP
BlockEnd = BlockSize
BlockSize = BlockSize SHL 1
LOOP
END SUB
SUB FFT(InArray() AS Complex, OutArray() AS Complex)
' ------------------------------------------------------------------
' Calculates the Fast Fourier Transform of the array of
' complex numbers 'InArray' to produce the output complex
' numbers in 'OutArray'
' ------------------------------------------------------------------
FourierTransform TwoPi, InArray(), OutArray()
END SUB
SUB IFFT(InArray() AS Complex, OutArray() AS Complex)
' ------------------------------------------------------------------
' Calculates the Inverse Fast Fourier Transform of the array of
' complex numbers represented by 'InArray' to produce the output
' complex numbers in 'OutArray'
' ------------------------------------------------------------------
FourierTransform -TwoPi, InArray(), OutArray()
IF ErrCode <> FFTOk THEN EXIT SUB
DIM AS INTEGER NumSamples, I
NumSamples = UBOUND(InArray) + 1
' Normalize the resulting time samples
FOR I = 0 TO NumSamples - 1
OutArray(I).X = OutArray(I).X / NumSamples
OutArray(I).Y = OutArray(I).Y / NumSamples
NEXT I
END SUB
FUNCTION CalcFrequency(BYVAL FrequencyIndex AS INTEGER, _
InArray() AS Complex) AS Complex
' ------------------------------------------------------------------
' This function calculates the complex frequency sample at a given
' index directly. Use this instead of 'FFT' when you only need one
' or two frequency samples, not the whole spectrum.
'
' It is also useful for calculating the Discrete Fourier Transform
' (DFT) of a number of data which is not an integer power of 2.
' For example, you could calculate the DFT of 100 points instead
' of rounding up to 128 and padding the extra 28 array slots with
' zeroes.
' ------------------------------------------------------------------
IF LBOUND(InArray) <> 0 THEN
ErrCode = FFTBoundErr
EXIT FUNCTION
END IF
ErrCode = FFTOk
DIM AS INTEGER NumSamples, K
DIM AS DOUBLE Cos1, Cos2, Cos3, Theta, Beta, Sin1, Sin2, Sin3
DIM AS Complex Res
NumSamples = UBOUND(InArray) + 1
Res.X = 0
Res.Y = 0
Theta = TwoPi * FrequencyIndex / NumSamples
Sin1 = SIN(- 2 * Theta)
Sin2 = SIN(- Theta)
Cos1 = COS(- 2 * Theta)
Cos2 = COS(- Theta)
Beta = 2 * Cos2
FOR K = 0 TO NumSamples - 1
' Update trig values
Sin3 = Beta * Sin2 - Sin1
Sin1 = Sin2
Sin2 = Sin3
Cos3 = Beta * Cos2 - Cos1
Cos1 = Cos2
Cos2 = Cos3
Res.X = Res.X + InArray(K).X * Cos3 - InArray(K).Y * Sin3
Res.Y = Res.Y + InArray(K).Y * Cos3 + InArray(K).X * Sin3
NEXT K
RETURN Res
END FUNCTION
5
Code: Select all
declare function fft1(a1() as Complex,b1() as Complex) as integer
declare function fft2(a2() as Complex , b2() as Complex) as integer
declare function fft3(a3() as Complex,b3() as Complex) as integer
'
declare function ft1test(n as integer) as integer
declare function ft2test(n as integer,m as integer) as integer
declare function ft3test(n as integer,m as integer, p as integer) as integer
declare function ft3gendat(b4() as single ,n as integer,m as integer, p as integer) as integer
declare function ft3Prtdat(b4() as single , n as integer,m as integer, p as integer) as integer
'
'
' Subroutine to call .
'
'
' FFT(InArray() AS Complex, OutArray() AS Complex)
'
dim as integer n,m,p,fg
6
Code: Select all
function fft1(a1() as Complex,b1() as Complex) as integer
'
' One dimensional , complex fft
'
' Input array a1()
' Output array b1()
'
'
FFT(a1() , b1())
'
return 0
'
end function
'
' ------------------------------------------------------------
'
function fft2(a2() as Complex , b2() as Complex) as integer
'
' Two dimensional , complex fft
'
' Input array a2()
' Output array b2()
'
'
static as integer i,j,n,m
n = ubound(a2)
m= ubound(a2,2)
'
'print " n , m ";n;" , ";m
n=n+1
m=m+1
'
redim as complex a1(0 to n-1),b1(0 to n-1)
'
for j=0 to m-1
for i=0 to n-1
b2(i,j).X = a2(i,j).X
b2(i,j).Y = a2(i,j).Y
next i
next j
'
'
for j=0 to m-1
'
for i=0 to n-1
a1(i).X = b2(i,j).X
a1(i).Y = b2(i,j).Y
next i
'
FFT(a1() , b1())
'
for i=0 to n-1
b2(i,j).X = b1(i).X
b2(i,j).Y = b1(i).Y
next i
'
next j
'
redim as complex a1(0 to m-1),b1(0 to m-1)
'
for i = 0 to n - 1
'
for j = 0 to m - 1
a1(j).X = b2(i,j).X
a1(j).Y = b2(i,j).Y
next j
'
FFT(a1() , b1())
'
for j=0 to m-1
b2(i,j).X = b1(j).X
b2(i,j).Y = b1(j).Y
next j
'
next i
'
return 0
'
end function
'
' ------------------------------------------------------------
'
function fft3(a3() as Complex,b3() as Complex) as integer
'
' Three dimensional , complex fft
'
' Input array a3()
' Output array b3() ' alternately , just return a3()
'
'
'
static as integer i, j , k
static as integer n , m , p
n = ubound(a3)
m = ubound(a3,2)
p = ubound(a3,3)
'
redim as Complex a1(0 to 4) ,b1(0 to 4)
'
n = n + 1
m = m + 1
p = p + 1
'
' Copy a3() to b3()
'
for i=0 to n-1
for j=0 to m-1
for k=0 to p-1
b3(i,j,k).X = a3(i,j,k).X
b3(i,j,k).Y = a3(i,j,k).Y
next k
next j
next i
'
' ---------- [ i , [ j , k ] ] ---------------
'
redim as Complex a1(0 to n-1) ,b1(0 to n-1)
'
for k=0 to p-1
for j=0 to m-1
'
for i=0 to n-1
a1(i)=b3(i,j,k)
next i
'
FFT(a1(),b1())
'
for i=0 to n-1
b3(i,j,k)=b1(i)
next i
'
next j
next k
'
' ................... .[j,[i,k]] {k,i,j}
'
redim as Complex a1(0 to m-1) ,b1(0 to m-1)
'
for k=0 to p-1
for i=0 to n-1
'
for j = 0 to m-1
a1(j) = b3(i,j,k)
next j
'
FFT(a1() , b1())
'
for j = 0 to m-1
b3(i,j,k) = b1(j)
next j
'
next i
next k
'
' ...........
'
redim as Complex a1(0 to p-1) ,b1(0 to p-1)
'
for j=0 to m-1
for i=0 to n-1
'
for k = 0 to p-1
a1(k) = b3(i,j,k)
next k
'
FFT(a1() , b1())
'
for k = 0 to p-1
b3(i,j,k) = b1(k)
next k
'
next i
next j
'
redim as Complex a1(0 to 4) ,b1(0 to 4)
'
return 0
'
end function
'
' -----------------------------------------------------------
'
function ft1test(n as integer) as integer
'
' Test fft1
'
'
static as integer i,fg
redim as Complex a1(0 to n-1), b1(0 to n-1)
'
for i=0 to (n/2)-1
a1(i).X = 1
next i
'
print
for i=0 to n - 1
print " "; a1(i).X;" , ";a1(i).Y
next i
'
print
'
FFT(a1() , b1())
'fg= fft1(a1() ,b1() )
'
print
for i=0 to n - 1
print " "; b1(i).X;" , ";b1(i).Y
next i
'
print
'
redim as Complex a1(0 to 3), b1(0 to 3)
'
return 0
'
end function
'
' ----------------------------------------------------------
'
function ft2test(n as integer,m as integer) as integer
'
' Test fft2
'
'
static as integer i,j,fg
static as single fp
redim as Complex a2(0 to n-1,0 to m-1), b2(0 to n-1,0 to m-1)
'
for j = 0 to (m/2) - 1
for i = 0 to (n/2) - 1
a2(i,j).X = 1
next i
next j
'
print
print " .........................................................."
print " 2d data , real "
print
for j = 0 to m - 1
for i = 0 to n - 1
print " ";a2(i,j).X ;
next i
print
next j
print
'
fg = fft2(a2() ,b2() )
'
print " .........................................................."
print " fft2d , real "
print
for j = 0 to m - 1
for i = 0 to n - 1
fp = b2(i,j).X
fp=int(fp*1000)/1000
print " ";fp ;
next i
print
next j
print
'
print " .........................................................."
print " fft2d , imag "
print
for j = 0 to m - 1
print
for i = 0 to n - 1
fp = b2(i,j).Y
fp=int(fp*1000)/1000
print " ";fp ;
next i
next j
print
'
redim as Complex a2(0 to 3,0 to 3), b2(0 to 3,0 to 3)
'
return 0
'
end function
function ft3test(n as integer,m as integer, p as integer) as integer
'
' Test fft3
'
'
static as integer i,j,k,fg
static as single fp,ep
redim as Complex a3(0 to n-1,0 to m-1,0 to p-1), b3(0 to n-1,0 to m-1,0 to p-1)
'
for k=0 to (p/2)-1
for j = 0 to (m/2) - 1
for i = 0 to (n/2) - 1
a3(i,j,k).X = 1
next i
next j
next k
'
'
'
print
print " ------------------------------------------"
print
print " 3d data , real "
print
for k=0 to p- 1
for j = 0 to m - 1
for i = 0 to n - 1
print " ";a3(i,j,k).X ;
next i
print
next j
print
next k
print
'
fg = fft3(a3() ,b3() )
'
print
print" ..........................................................."
print " fft3d , real "
print
for k=0 to p-1
for j = 0 to m - 1
for i = 0 to n - 1
fp = b3(i,j,k).X
fp=int(fp*1000)/1000
print " ";fp ;
next i
print
next j
print
next k
print
'
print
print " .........................................................."
print " fft3d , imag "
print
for k=0 to p-1
for j = 0 to m - 1
for i = 0 to n - 1
fp = b3(i,j,k).Y
fp=int(fp*1000)/1000
print " ";fp ;
next i
print
next j
print
next k
print
'
'
print
print " ..........................................................."
print " fft3d , mag "
print
for k=0 to p-1
for j = 0 to m - 1
for i = 0 to n - 1
ep = b3(i,j,k).X
fp = b3(i,j,k).Y
fp=ep*ep+fp*fp
if (fp>0) then fp=sqr(fp)
fp=int(fp*1000)/1000
print " ";fp ;
next i
print
next j
print
next k
print
'
redim as Complex a3(0 to 3,0 to 3,0 to 3), b3(0 to 3,0 to 3,0 to 3)
'
return 0
'
end function
'
' ------------------------------------------------------------
'
function ft3gendat(b4() as single , n as integer,m as integer, p as integer) as integer
'
' Test fft3
'
'
static as integer i,j,k,fg
static as single fp,ep
redim as Complex a3(0 to n-1,0 to m-1,0 to p-1), b3(0 to n-1,0 to m-1,0 to p-1)
redim as single b4(0 to n-1,0 to m-1,0 to p-1)
'
for k=0 to (p/8)-1
for j = 0 to (m/3) - 1
for i = 0 to (n/4) - 1
a3(i,j,k).X = 1
next i
next j
next k
'
fg = fft3(a3() ,b3() )
'
for k=0 to p-1
for j = 0 to m - 1
for i = 0 to n - 1
ep = b3(i,j,k).X
fp = b3(i,j,k).Y
fp=ep*ep+fp*fp
if (fp>0) then fp=sqr(fp)
fp=int(fp*1000)/1000
b4(i,j,k) = fp
next i
next j
next k
'
redim as Complex a3(0 to 3,0 to 3,0 to 3)
' , b3(0 to 3,0 to 3,0 to 3)
'
return 0
'
end function
'
' ------------------------------------------------------------
'
function ft3Prtdat(b4() as single , n as integer,m as integer, p as integer) as integer
'
' Print data from array b3()
'
static as single fp
static as integer i,j,k
'
print
print " ..........................................................."
print " ft3Prtdat , fft3d , mag "
print
for k=0 to p-1
for j = 0 to m - 1
for i = 0 to n - 1
fp = b4(i,j,k)
print " ";fp ;
next i
print
next j
print
next k
print
'
return 0
'
end function
'
' --------------------------------------------------------------
'