inverse FFT

[Gabi]

@top204

Here is some generic code referencing inverse FFT.
- have added some comments where I though appropriate
- as a side note, all Real and Imaginary values originally using 64bit floats
- I have tested the code in a PC application and it seems ok
- there is no filter function implemented, just inverse FFT calc
- there might be some functions in there that probably need be implemented, ex: Floor(float)

Code Select Expand


' Defines/Variables

dim FFT_binLength As Int = 33 ' the size of FFT bins.  Must be a power of 2 + 1
dim N, N2, N_2, N_21, NU As Int ' integers
dim invertedFFT(), FFTReal(FFT_binLength), FFTImag(FFT_binLength), fReal() As Double ' arrays of double (64bit floating number)
dim tReal, tImag  As Double ' 64bit floating numbers
' note:
' invertedFFT() size will be = (FFT_binLength -1) * 2
' fReal() is an empty array, or say 1 element if you want. Use what you need as this returns
' when N does not match a power of 2





Main:
'
' lame Test to pre-fill Real and Imaginary arrays with some FFT values
' in real app use actual values got via forward FFT, or such
For i = 0 To (FFT_binLength - 1)
FFTReal(i) = 0.04562232
FFTImag(i) = -0.0954 ' the imaginary could probably be passed as zeroes, need to check
Next
'
' Calculate inverse FFT
' and get the Result filling the array that will hold the iFFT (real)values
invertedFFT = Inverse(FFTReal, FFTImag)





'--------------------------------------------------------------
' Calculates the inverse FFT (Fast Fourier Transformation)
' params:
'    Real() = the real part of the frequency (array)
'    Imag() = the imaginary part of the frequency (array)
'    The length of Real and Imag must be a power of 2 + 1 // like 33, 65 etc
' Returns:
'    the array of Real values of the time signal
'--------------------------------------------------------------
Sub Inverse(Real() As Double, Imag() As Double) As Double()
dim i, j As Int
dim ld As Double
dim ang As Double
dim k, l, p, nu1 As Int
dim c, s As Double
dim r As Int

N = (Real.Length - 1) * 2

ld = Logarithm(N, 2) ' get the logarithm in Base:2
If ld - Floor(ld) > 0 Then ' get the Floor function
error_message = ("N must be a power of 2") ' show the error msg
Return fReal ' this array is an empty one, can use whatever error sig you need
End If

NU = ld
N2 = N / 2
N_2 = N2
N_21 = N2 + 1

dim xReal(N), xImag(N), invReal(N), mSin(N_21), mCos(N_21) As Double ' arrays of double (64bit floating number)

xReal(0) = Real(0)
xImag(0) = Imag(0)
For i = 1 To N2 - 1
xReal(i) = Real(i)
xImag(i) = Imag(i)
j = n - i
xReal(j) = Real(i)
xImag(j) = -Imag(i)
Next
xReal(N2) = Real(N2)
xImag(N2) = Imag(N2)

For i = 0 To N2
ang = 2 * cPI * i / N ' cPi = PI constant number 3.14xxxxx
mCos(i) = Cos(ang) ' param of cos  in radians, returns radians
mSin(i) = Sin(ang) ' param of sin  in radians, returns radians
Next

'First phase - calculation
nu1 = NU - 1
k = 0
For l = 1 To NU
Do While k < N
For i = 1 To N2
p = BitReverse(Bit.ShiftRight(k, nu1))
c = mCos(p)
s = mSin(p)
tReal = xReal(k + N2) * c + xImag(k + N2) * s
tImag = xImag(k + N2) * c - xReal(k + N2) * s
xReal(k + N2) = xReal(k) - tReal
xImag(k + N2) = xImag(k) - tImag
xReal(k) = xReal(k) + tReal
xImag(k) = xImag(k) + tImag
k = k + 1
Next
k = k + N2
Loop
k = 0
nu1 = nu1 - 1
N2 = N2 / 2
Next

'Second phase - recombination
k = 0
Do While k < N
r = BitReverse(k)
If r > k Then
tReal = xReal(k)
tImag = xImag(k)
xReal(k) = xReal(r)
xImag(k) = xImag(r)
xReal(r) = tReal
xImag(r) = tImag
End If
k = k + 1
Loop

For i = 0 To N - 1
invReal(i) = xReal(i) / 2
Next

Return invReal
End Sub

'--------------------------------------------------------------
'
'--------------------------------------------------------------
Sub BitReverse(j As Int) As Int
Private j1, j2, k As Int

j1 = j
k = 0
For i = 1 To NU
j2 = j1 / 2
k = 2 * k + j1 - 2 * j2
j1 = j2
Next
Return k
End Sub