The beauty / magic of math (Vol. 1 - 19 build 2023-05-01)

[UEZ]

A little presumptuous of me.

Might you be able to help me with some mathematics.

Mostly sorted, just a fiddly bit , relates to DFT and CFT .

I would suggest that you open a topic under General and maybe I or someone else can help you.

Currently, I've a COVID infection...

[fxm]

The Fourier Transform (FT) makes it possible to go from the temporal space to the frequency space.
Discrete Fourier Transform (DFT) is used to process a digital (sampled) signal. It is a discrete equivalent of the Continuous Fourier Transform (CFT) used to process an analog (continuous) signal.
The Fast Fourier Transform (FFT) is an efficient algorithm for calculating the DFT.

[dafhi]

[updated]

i copied a 'FFT from someone else and put it in my freebasic softsynth

Code: Select all

const pi = 4 * atn(1)

'' ---------------
'    -~- FFT -~-
' ----------------


''2 ^ 10 = 1024
''2 ^ 18 = 262144
Dim As uByte    lutpow = 14

Dim Shared As Integer MaxPowerOfTwo = 20
If lutpow > MaxPowerOfTwo Then lutpow = MaxPowerOfTwo

Dim As Integer   size_PADLUT = 1 Shl lutpow

Function NumberOfBitsNeeded(ByVal PowerOfTwo as Integer) as Integer
    for I As Integer= 0 to MaxPowerOfTwo
        if (PowerOfTwo and (1 shl I)) <> 0 then
            return I
        End If
    Next
    Return 0
End Function
Function ReverseBits(ByVal Index As Integer,ByVal NumBits As Integer) As Integer
  Dim As Integer Rev
  for I As Integer = 0 to NumBits - 1
    Rev Shl= 1
    Rev Or= Index And 1
    Index Shr= 1
  Next
  Return Rev
End Function

Dim Shared As Integer RevBits(size_PADLUT)

Dim As Integer BitsNeeded = NumberOfBitsNeeded(size_PADLUT)
For I As Integer = 0 To size_PADLUT
  RevBits(I) = ReverseBits(I, BitsNeeded)
Next

Function IsPowerOfTwo(ByVal X as Integer) As Integer
  Dim As Integer  I, Y
  Y = 2
  for I = 1 To MaxPowerOfTwo
    if X = Y then Return TRUE
    Y Shl= 1
  Next
  Return FALSE
End Function

Sub FourierTransform(ByVal AngleNumerator As Single, _
                     ByVal NumSamples As Integer, _
                     lpReal As Single Ptr, _
                     lpImag As Single Ptr, _
                     lpRealOut As Single Ptr, _
                     lpImagOut As Single Ptr)
                     
' http://www.koders.com/delphi/fidB6DD10205DAFD71EFF5093D4916237BB2801DD8D.aspx
  
  Dim As Integer I, J, K, N, BlockSize, BlockEnd', NumBits
  Dim As single  Delta_angle, Delta_ar
  Dim As single  Alpha, Beta
  Dim As single  Tr, Ti, Ar, Ai
  
  if not IsPowerOfTwo(NumSamples) or (NumSamples < 2) then
      ?"Error in procedure Fourier:  NumSamples="; NumSamples
      ?" is not a positive integer power of 2."
      Exit Sub
  End If
  
  'NumBits = NumberOfBitsNeeded(NumSamples)
  for I = 0 to NumSamples
    J = RevBits(I)
    lpRealOut[J] = lpReal[I]
    lpImagOut[J] = lpImag[I]
  Next
  
  BlockEnd = 1
  BlockSize = 2
  while BlockSize <= NumSamples
    Delta_angle = AngleNumerator / BlockSize
    Alpha = Sin(0.5 * Delta_angle)
    Alpha = 2.0 * Alpha * Alpha
    Beta = Sin(Delta_angle)

    I = 0
    while I < NumSamples
      Ar = 1.0 ''cos(0)
      Ai = 0.0 ''sin(0)

      J = I
      for N = 0 to BlockEnd - 1
          K = J + BlockEnd
          Tr = Ar * lpRealOut[K] - Ai * lpImagOut[K]
          Ti = Ar * lpImagOut[K] + Ai * lpRealOut[K]
          lpRealOut[K] = lpRealOut[J] - Tr
          lpImagOut[K] = lpImagOut[J] - Ti
          lpRealOut[J] += Tr
          lpImagOut[J] += Ti
          Delta_ar = Alpha * Ar + Beta * Ai
          Ai -= (Alpha * Ai - Beta * Ar)
          Ar -= Delta_ar
          J += 1
      Next

      I += BlockSize
    Wend
    
    BlockEnd = BlockSize
    BlockSize Shl= 1
  Wend
End Sub



type signal_and_transform
  declare sub resize( as long )
  declare sub FFT( as long )
  as single   sr(any)
  as single   si(any)
  as single   dr(any)
  as single   di(any)
end type

sub signal_and_transform.resize( c as long )
  redim sr(c-1)
  redim si(c-1)
  redim dr(c-1)
  redim di(c-1)
end sub

sub signal_and_transform.FFT( c_samps as long )
      FourierTransform( _
    2 * PI, c_samps, @sr(0), @si(0), @dr(0), @di(0) )
end sub


dim as signal_and_transform   signal

var num_samples = 2 ^ 9
signal.resize num_samples
for i as long = 0 to num_samples - 1
  signal.sr(i) = (rnd - .5)
  signal.si(i) = (rnd - .5)
next

signal.FFT num_samples

var w = 800
var h = 600
screenres w, h
for i as long = 0 to num_samples - 1
  pset ( i, h/2*( 1 + signal.dr(i)/sqr(num_samples) ) )
next

sleep

[srvaldez]

hello dafhi
I get Aborting due to runtime error 12 ("segmentation violation" signal) in softsynth.bas::FOURIERTRANSFORM()

[dodicat]

Line 17
Dim shared As Integer size_PADLUT:size_PADLUT = 1 Shl lutpow
and
sub signal_and_transform.resize( c as long )

redim sr(size_PADLUT)
redim si(size_PADLUT)
redim dr(size_PADLUT)
redim di(size_PADLUT)
end sub
To get it to work.
i.e. the single pointers need to point to these larger arrays.
Maybe not all of them, dafhi should check each one.

[srvaldez]

thanks dodicat
runs ok without crashing
@dafhi, I suggest compiling with -exx option until your program runs without issues

[Luxan]

Firstly the definition of the continuous inverse Fourier Transform:

f(t) = integrate(F(w)*exp(i*w*t),w), where w=2*pi*f, f =frequency .

this is similar to the next equation.

Secondly, the continuous integral we want convert to
an equivalent Fourier Transform representation, then to a
DFT or FFT.

P(t) = integrate(p(x)*exp(i*Gfe*x*t),x)

This implies that x == w.
and

f(t) = integrate(F(x)*exp(i*x*t),x)

this just leaves the multiplier Gfe, in the exponent, what's the
next step in aligning these equations; is there some exponential
or Fourier transform identity we can use ?

[Luxan]

I've been considering down sampling , whereby you skip every other
sample, to align the outputs from the two integrals described.

I used a function :

f(j) = (n-j)*sin(12*j*pi/n)

as the input to a n point DFT, n=256 .

I then sampled the output from this at every
second value, saving this to an array of size
n/2 .

The results were in agreement, in this instance Gfe = 2 .
For Gfe < 1 , a different approach is required.

[dafhi]

thanks for the report guys. i tried -exx it still runs.

looks like Luxan's getting 'er done.

[Luxan]

I examined this using wxmaxima cas, note the restricted domain
of f(j) .


The procedure :

1. f(j) ,[j,0,n-1]
2. An fft of size n*m, equivalent to zero extension
3. Sample the fft at every m ,[ ,1,n*m]
4. inverse fft of size n
5. result, the same as f(j) ,[j,0,n-1], from 1

Gfe==1/m