How to code sound FFT ?

[Jim]

All those squiggles are trying to tell you that the frequency at a given point in time is the sum of all the sub-frequencies in your signal.  Don't let it bother you too much...

In answer to your questions.
a) MaxPartials/partials is how much detail you want to analyse - the number of sub-frequencies you want.  Set it to something like 128.  For every incoming chunk of signal (say, 1/50th second's worth of audio), you get an ARRAY of MaxPartials as the output.
b) yes, you need to zero those arrays.  if pb does that for you then you don't need to do it too.
c) the output is an array.  it's stored in 'amp'.  If you chose 128 partials, then effectively you are dividing the incoming signal into 128 frequencies.  The values in 'amp' are the magnitudes of each of those sub-frequencies.

Code: [Select]

   For i = 0 to framzesize
Typo?

If you plot the array vertically against time, (so, 128 dark to light pixels based on the values in 'amp' vertically per 1/50th of a second) you will get a pretty picture - something you might well recognise.

If it's too slow to start with, drop the partials to 64 or 32.

You don't need an FFT library, you have everything you need!

Jim

[va!n]

Since yesterday i am sitting by my best friend in duisburg and we are trying to code and finish some missing routines for our actually project like a good and fast FFT routine. From yesterday until this morning we have tried to get convert and working a FFT version, which seems to work very nice. *yeah*

So i will post the source(s) to help hopefully other peoples and possible some other guys here may be help to optimize the code in speed =)



Here is the original VB sourcecode we have tried to convert:

Code: [Select]

Dim REX(512) 'REX[ ] holds the real part of the frequency domain
 Dim IMX(512) 'IMX[ ] holds the imaginary part of the frequency domain
 Const N = 512
Public Sub fft()

pi = 3.14159265 'Set constants

1000 'THE FAST FOURIER TRANSFORM
'copyright © 1997-1999 by California Technical Publishing
'published with  permission from Steven W Smith, www.dspguide.com
'GUI by logix4u , www.logix4u.net
'modified by logix4u, www.logix4.net
1010 'Upon entry, N% contains the number of points in the DFT, REX[ ] and
1020 'IMX[ ] contain the real and imaginary parts of the input. Upon return,
1030 'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1.
1060 NM1% = N% - 1
1070 ND2% = N% / 2
1080 M% = CInt(Log(N%) / Log(2))
1090 J% = ND2%
1100 '
1110 For i% = 1 To N% - 2 'Bit reversal sorting
1120 If i% >= J% Then GoTo 1190
1130 TR = REX(J%)
1140 TI = IMX(J%)
1150 REX(J%) = REX(i%)
1160 IMX(J%) = IMX(i%)
1170 REX(i%) = TR
1180 IMX(i%) = TI
1190 K% = ND2%
1200 If K% > J% Then GoTo 1240
1210 J% = J% - K%
1220 K% = K% / 2
1230 GoTo 1200
1240 J% = J% + K%
1250 Next i%
1260 '
1270 For L% = 1 To M% 'Loop for each stage
1280 LE% = CInt(2 ^ L%)
1290 LE2% = LE% / 2
1300 UR = 1
1310 UI = 0
1320 SR = Cos(pi / LE2%) 'Calculate sine & cosine values
1330 SI = -Sin(pi / LE2%)
1340 For J% = 1 To LE2% 'Loop for each sub DFT
1350 JM1% = J% - 1
1360 For i% = JM1% To NM1% Step LE% 'Loop for each butterfly
1370 IP% = i% + LE2%
1380 TR = REX(IP%) * UR - IMX(IP%) * UI 'Butterfly calculation
1390 TI = REX(IP%) * UI + IMX(IP%) * UR
1400 REX(IP%) = REX(i%) - TR
1410 IMX(IP%) = IMX(i%) - TI
1420 REX(i%) = REX(i%) + TR
1430 IMX(i%) = IMX(i%) + TI
1440 Next i%
1450 TR = UR
1460 UR = TR * SR - UI * SI
1470 UI = TR * SI + UI * SR
1480 Next J%
1490 Next L%
1500 '
End Sub



Here is the converted ( /2 instead >> 1 unoptimized) PB version, which seems to work fine:

Code: [Select]

pi.f = 3.14159265
N = 8 ; Number of samples

Dim rex.f( N - 1 )
Dim imx.f( N - 1 )

; Upon entry, N contains the number of points in the DFT, REX() and
; IMX() contain the real and imaginary parts of the input. Upon return,
; REX() and IMX() contain the DFT output. All signals run from 0 to N-1.

NM1.l = N - 1
ND2.l = N / 2
M.l = Int( Log( N ) / Log( 2 ) )
J.l = ND2

For i.l = 1 To N - 2                                ; Bit reversal sorting
    If i < J
        TR.f = REX( J )
        TI.f = IMX( J )
        REX( J ) = REX( i )
        IMX( J ) = IMX( i )
        REX( i ) = TR
        IMX( i ) = TI
    EndIf
    K = ND2
    While K <= J
        J = J - K
        K = K / 2
    Wend
    J = J + K
Next

For L = 1 To M                                      ; Loop for each stage
    LE.l = Int( Pow( 2, L ) )
    LE2.l = LE / 2
    UR.f = 1
    UI.f = 0
    SR.f = Cos( #PI / LE2 )                         ; Calculate sine & cosine values
    SI.f = - Sin( #PI / LE2 )
    For J.l = 1 To LE2                              ; Loop for each sub DFT
        JM1.l = J - 1
        For i = JM1 To NM1                          ; Loop for each butterfly
            IP.l = i + LE2
            TR = REX( IP ) * UR - IMX( IP ) * UI     ; Butterfly calculation
            TI = REX( IP ) * UI + IMX( IP ) * UR
            REX( IP ) = REX( i ) - TR
            IMX( IP ) = IMX( i ) - TI
            REX( i  ) = REX( i ) + TR
            IMX( i  ) = IMX( i ) + TI
            i + LE - 1
        Next
        TR = UR
        UR = TR * SR - UI * SI
        UI = TR * SI + UI * SR
    Next
Next



And at least here is my latest try to convert a small modified version from C to PB:

Code: [Select]

; / ***************************************************************************/
; / * fft.pb                                                                  */
; / ***************************************************************************/
; / * Forward And Reverse Fast Fourier Transform routines.                    */
; / * ====================================================                    */
; / * Original Source by Tobias Dammers (c) 2005.                             */
; / * PB version by Thorsten Will aka va!n in 2008                            */
; / *                                                                         */
; / * Reference: http://wiki.allegro.cc/FFT                                   */
; / *                                                                         */
; / * Feel free To use these routines As you see fit.                         */
; / ***************************************************************************/

; // Fourier Transform is the process of transforming signals from the time domain
; // (levels over time) into the frequency domain (levels over frequency). The dis-
; // crete-time version of this is called the Discrete Fourier Transform (DFT).
; // The standard approach is a rather costly operation, but For blocks of Data
; // whose length is a power of 2, a dramatic optimization is possible, resulting
; // in an algorithm called Fast Fourier Transform, FFT. This is a C implementation
; // of the algorithm, And of its reverse process.


; // The forward FFT. rex holds the real part of the signal, imx the imaginary part.
; // Both arrays must be at least 2^m samples long. The result is returned in the
; // same two arrays, ie the transform is in-place. After transforming, rex holds
; // the real part of the frequency domain, And imx the imaginary part.
; // To calculate the frequency spectrum of an arbitrary block of wave Data, put the
; // input in rex[], And set imx[] To all-zero.

Global Dim rex.f( 512 )
Global Dim imx.f( 512 )

Procedure FFT( m.l )

    Protected nm1.l, nd2.l
    Protected i.l, j.l, l.l, n.l
    Protected tr.f, ti.f
   
    n.l   = 1 << m
    nm1.l = n - 1
    nd2.l = n >> 1
    j.l   = nd2
    tr.f
    ti.f
    i.l
    l.l

    i = 1
    While  i < n   
        k.l
       
        If ( i < j )
            tr = rex( j )
            ti = imx( j )
            rex( j ) = rex( i )
            imx( j ) = imx( i )
            rex( i ) = tr
            imx( i ) = ti
        EndIf

        k = nd2

        While ( ( k <= j ) And ( k > 0 ) )     ; &&
            j = j - k
            k = k >> 1
        Wend
        j = j + k
        i +1
    Wend

 
    l = 1
    While l <= m   
        le.l  = 1 << l
        le2.l = le >> 1
        ur.f  = 1.0
        ui.f  = 0.0
        sr.f  =  Cos( #PI / le2 )
        si.f  = -Sin( #PI / le2 )
                   
        j = 1
        While j <= le2
            i = j-1
            While  i <= nm1   
                ip.l = i + le2
                tr = rex( ip ) * ur - imx( ip ) * ui
                ti = rex( ip ) * ui + imx( ip ) * ur
                rex( ip ) = rex( i ) - tr
                imx( ip ) = imx( i ) - ti
                rex( i  ) = rex( i ) + tr
                imx( i  ) = imx( i ) + ti
                i = i + le
            Wend
            tr = ur
            ur = tr * sr - ui * si
            ui = tr * si + ui * sr
            j + 1
        Wend
        l + 1
    Wend

EndProcedure



; // The inverse FFT. This function takes advantage of the special symmetry of
; // Fourier pairs; all we need to do is scale the input, apply a forward FFT
; // then scale the output.
; // As a test, you can feed a signal through the forward fft, then feed the result
; // back through the inverse fft, giving you the original signal, give Or take
; // some roundoff errors.

Procedure Inv_FFT( m.l )

    n.l = 1 << m
    i.l
   
    i = 0
    While  i < n
        imx( i ) = -imx( i )
        i + 1
    Wend

    FFT( m )

    i = 0
    While i < n
        rex( i ) = rex( i ) / n       ; rex[i] /= (float)n;
        imx( i ) = imx( i ) / -n      ; imx[i] /= -(float)n
        i + 1
    Wend

EndProcedure


If some of you have a benchmark for this both routines (which is the fastest) or how and where to optimize the code (basic/ x86 or MMX or SSEx) it would be really cool. Thanks for any feedback. Hope this source(s) may be helpfull for some guys here.

[taj]

Thats a karma++ moment.

[Jim]

Va!n - there are lots of places in your code where you divide by something.  Often that division could be replaced by a multiply by the reciprocal (ie.  x = a/n, y=b/n  could be re-written  t = 1/n, x=a*t, y=b*t).  And often, that reciprocal calculation can be hauled outside the loops.

Also, you should replace the separate calls to Sin and Cos library functions to use the inline assembly instruction fsincos which calculates both values at the same time.
For instance, I have this code

Code: [Select]

static void sincos(double a, double *s, double *c)
{
double zin,coz;

//a *= 3.141592653589/180.0;
_asm
{
fld [a]
fsincos
fstp [coz]
fstp [zin]
}

*s = zin;
*c = coz;
}
in most of my demos.  You need the multiply if you're working in degrees (which you're not).

Jim

[va!n]

@taj:
Thanks for k++ and feedback =)

@JIM:
Thanks for feedback and the great SIN/COS inline asm hint, i will try to add this later. *thumbs up* K++ for the tipp! =)

@all:
I managed it to optimize the FFT code a little bit more. I will try to optimize the SIN/COS stuff with inline asm later.

Code: [Select]

; // *************************************************************************************
; // *
; // * Project:   FFT (THE FAST FOURIER TRANSFORM) v0.1b   
; // *
; // * Original version:  logix4u, www.logix4.net
; // * Coverted version:  Tranquil and Mr.Vain of Secretly! (Thorsten Will) in 2008
; // * Optimisations by:  va!n aka Thorsten Will
; // *
; // * Reference:
; // * http://logix4u.net/DSP/Fast_Fourier_Transform/Visual_Basic_program_for_Fast_Fourier_Transform.html
; // *
; // * Upon entry, N.w contains the number of points in the DFT, REX() and
; // * IMX() contain the real and imaginary parts of the input. Upon return,
; // * REX() and IMX() contain the DFT output. All signals run from 0 to N.w-1
; // *
; // *************************************************************************************

Procedure record_doFFT( *scope.SCOPE )

    ; // -------- Init some values for FFT analysing --------

    N.w   = 1024                ; // Number of samples
    NM1.l = N - 1
    ND2.l = N >> 1              ; // Optmimized, instead N / 2
    M.l   = 10                  ; // If N.w = 1024 == 10 == Same as:  Int(Log(N) / 0.69314718055994529)
    J.l   = ND2

    If Config\buffer = 0 : ProcedureReturn : EndIf

    ; // -------- Clear and Fill array values for analysing in just only one loop --------
         
    For i = 0 To Config\size Step 2       ; // Optimized by merging clear and fill array in one loop
        rex( i >> 1          ) = 0         ; //   0 to  512
        imx( i >> 1          ) = 0         ; //   0 to  512
        rex( i >> 1 + N >> 1 ) = 0         ; // 513 to 1024
        imx( i >> 1 + N >> 1 ) = 0         ; // 513 to 1024
       
        value.w    = PeekW( Config\buffer + i )
;       value.w    = PeekW( Config\buffer + i + *scope\channel * 2 )  ; // Enable this For Stereo Inpus
        rex( i >> 1 ) = value / 32767                                 ; // Optimized by doing i >> 1
;       imx( i >> 1 ) = 0                                             ; // Can be removed now
    Next

    ; // -------- Start FFT --------

    For i.l = 1 To N - 2                                ; // Bit reversal sorting
        If i < J
            TR.f = REX( J )
            TI.f = IMX( J )
            REX( J ) = REX( i )
            IMX( J ) = IMX( i )
            REX( i ) = TR
            IMX( i ) = TI
        EndIf
       
        K = ND2
       
        While K <= J
            J = J - K
            K = K >> 1                                  ; // Optmimized, instead N / 2
        Wend
        J = J + K
    Next

    For L = 1 To M                                      ; // Loop for each stage
        LE.l = Int( Pow( 2, L ) )
        LE2.l = LE >> 1                                 ; // Optmimized, instead N / 2
        UR.f = 1
        UI.f = 0
        SR.f =   Cos( #PI / LE2 )                       ; // Calculate sine & cosine values
        SI.f = - Sin( #PI / LE2 )
        For J.l = 1 To LE2                              ; // Loop for each sub DFT
            JM1.l = J - 1
            For i = JM1 To NM1                          ; // Loop for each butterfly
                IP.l = i + LE2
                TR = REX( IP ) * UR - IMX( IP ) * UI    ; // Butterfly calculation
                TI = REX( IP ) * UI + IMX( IP ) * UR
                REX( IP ) = REX( i ) - TR
                IMX( IP ) = IMX( i ) - TI
                REX( i  ) = REX( i ) + TR
                IMX( i  ) = IMX( i ) + TI
                i + LE - 1
            Next i
            TR = UR
            UR = TR * SR - UI * SI
            UI = TR * SI + UI * SR
        Next
    Next

    ; // -------- Calculate Outputarray --------
   
    For cnt = 0 To N.w                                      ; fixed to N.w instead wrong fixed value
        outputarray( cnt )  =  ( IMX( cnt ) * IMX( cnt ) )  +  ( REX( cnt ) * REX( cnt ) ) 
    Next

    ; // -------- Search for MaxValue of the Paket --------

    maxvalue = 0
   
    For cnt = 0 To N.w                                      ; // fixed to N.w instead wrong fixed value
        If  maxvalue < outputarray( cnt )
            maxvalue = outputarray( cnt )
        EndIf
    Next
 
    ; // -------- Draw FFT -------- 
 
    StartDrawing( WindowOutput( FFTWnd ) )
        Box( 0, 0, N.w, 500, $0 )
        MaxPeak = 0
        For cnt = 0 To 500 ;                                   ; // Change fixed value to N.w !?
            DiffY = Outputarray( cnt ) / MaxValue * 400
            Box( cnt, 400, 1, -DiffY, $FFFFFF)
            If DiffY > MaxPeak
                MaxPeak = DiffY
                MaxPos  = cnt
            EndIf
        Next
    StopDrawing()

    SetWindowTitle( FFTWnd, ShowNote_Get( record_FindNote( MaxPos ) ) )

EndProcedure


[va!n]

Okay guys... i am sure some of you are interested in FFT coding, so here we did a small, fast test application (sound recorder) to test our FFT routine and its results. Here the full PB sourcecode. An executeable version will be added for testing too.

Sourcecode using our first "FFT v0.1" version:

Code: [Select]

Global Dim rex.f(512*2+1)
Global Dim imx.f(512*2+1)
Global Dim OutPutArray.f(512*2+1)


Global FFTWnd

Structure NoteRange
  Note.l
  FromPos.l
  ToPos.l
EndStructure

Procedure.l ShowNote_Init()

    Global Dim NoteRange.NoteRange(53)
   
    For Note=0 To 53
      Read FromPos.w
      Read ToPos.w
      NoteRange(Note)\FromPos = FromPos
      NoteRange(Note)\ToPos   = ToPos
      NoteRange(Note)\Note    = Note
    Next

    Global Dim g_RealNote.s( 53 )

    For i = 0 To 42+12-1
        Read sRealNote.s
        g_RealNote( i ) = sRealNote.s
    Next
EndProcedure

Procedure.s ShowNote_Get( lValue)
    ProcedureReturn g_RealNote.s( lValue )
EndProcedure

ShowNote_Init()

Structure SCOPE
  channel.b
  left.l
  top.l
  width.l
  height.l
  middleY.l
  quarterY.l
EndStructure

Structure CONFIG
 
  hWindow.l           ; Window handle
  size.l              ; Wave buffer size
  buffer.l            ; Wave buffer pointer
  output.l            ; WindowOutput()
  wave.l              ; Address of waveform-audio input device

  format.WAVEFORMATEX ; Capturing WaveFormatEx
  lBuf.l              ; Capturing Buffer size
  nBuf.l              ; Capturing Buffer number
  nDev.l              ; Capturing Device identifier
  nBit.l              ; Capturing Resolution (8/16)
  nHertz.l            ; Capturing Frequency  (Hertz)
  nChannel.l          ; Capturing Channels number (Mono/Stereo)
 
  LScope.SCOPE        ; Wave form display
  RScope.SCOPE        ; Wave form display

EndStructure

Global Config.CONFIG
Global Dim inHdr.WAVEHDR(16)

Config\format\wFormatTag = #WAVE_FORMAT_PCM

Procedure Record_Start()

  Config\format\nChannels = 1

  Config\format\wBitsPerSample  = 16
  Config\format\nSamplesPerSec  = 8000
  Config\nDev                   = 0 ; (0 default MS Sound Mapper device)
  Config\lBuf                   = 1024
  Config\nBuf                   = 8
  Config\nBit                   = 1
 
  Config\format\nBlockAlign     = (Config\format\nChannels*Config\format\wBitsPerSample)/8
  Config\format\nAvgBytesPerSec = Config\format\nSamplesPerSec*Config\format\nBlockAlign
 
  If #MMSYSERR_NOERROR = waveInOpen_(@Config\wave,#WAVE_MAPPER+Config\nDev,@Config\format,Config\hWindow,#Null,#CALLBACK_WINDOW|#WAVE_FORMAT_DIRECT)
   
    For i=0 To Config\nBuf-1
      inHdr(i)\lpData=AllocateMemory(Config\lBuf)
      inHdr(i)\dwBufferLength=Config\lBuf
      waveInPrepareHeader_(Config\wave,inHdr(i),SizeOf(WAVEHDR))
      waveInAddBuffer_(Config\wave,inHdr(i),SizeOf(WAVEHDR))
    Next
   
    If #MMSYSERR_NOERROR = waveInStart_(Config\wave)
      SetTimer_(Config\hWindow,0,1,0)
    EndIf
   
  EndIf
 
EndProcedure

Procedure Record_Read(hWaveIn.l,lpWaveHdr.l)
  *hWave.WAVEHDR=lpWaveHdr
  Config\buffer=*hWave\lpData
  Config\size=*hWave\dwBytesRecorded
  waveInAddBuffer_(hWaveIn,lpWaveHdr,SizeOf(WAVEHDR))
EndProcedure

Procedure record_FindNote(Value)
  For Note = 0 To 53
    If Value=>NoteRange(Note)\FromPos And Value<=NoteRange(Note)\ToPos
      ProcedureReturn note
    EndIf
  Next
EndProcedure

Procedure record_doFFT(*scope.SCOPE)

  If Config\buffer = 0 : ProcedureReturn : EndIf

  For pos=0 To 1024:rex(pos)=0:imx(pos)=0:Next
  pos = 0
  For i=0 To Config\size Step 2
    value.w=PeekW(Config\buffer+i)
    ;value.w=PeekW(Config\buffer+i+*scope\channel*2) Enable this for Stereo Inpus
    rex(pos) = value/32767
    imx(pos) = 0  :
    pos + 1
  Next

  N.w = 1024    ; Num Samples
  ;N.w = 512
  m-w = 102
  NM1.l = N - 1
  ND2.l = N / 2
  M.l = Int(Log(N) / 0.69314718055994529)
  J.l = ND2

  For i.l = 1 To N - 2 ; Bit reversal sorting
      If i < J
          TR.f = REX(J)
          TI.f = IMX(J)
          REX(J) = REX(i)
          IMX(J) = IMX(i)
          REX(i) = TR
          IMX(i) = TI
      EndIf
      K = ND2
      While K <= J
          J = J - K
          K = K / 2
      Wend
      J = J + K
  Next i

  For L = 1 To M ; Loop for each stage
      LE.l = Int(Pow(2, L))
      ;LE2.l = LE / 2
      LE2.l = LE >> 1
      UR.f = 1
      UI.f = 0
      SR.f = Cos(#PI / LE2) ; Calculate sine & cosine values
      SI.f = - Sin(#PI / LE2)
      For J.l = 1 To LE2 ; Loop for each sub DFT
          JM1.l = J - 1
          For i = JM1 To NM1 ; Loop for each butterfly
              IP.l = i + LE2
              TR = REX(IP) * UR - IMX(IP) * UI ; Butterfly calculation
              TI = REX(IP) * UI + IMX(IP) * UR
              REX(IP) = REX(i) - TR
              IMX(IP) = IMX(i) - TI
              REX(i) = REX(i) + TR
              IMX(i) = IMX(i) + TI
              i + LE - 1
          Next i
          TR = UR
          UR = TR * SR - UI * SI
          UI = TR * SI + UI * SR
      Next J
  Next L

  ; Outputarray berechnen
 
   For cnt=0 To 512
     outputarray(cnt) = (IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt))
   Next cnt

  ;Search for MaxValue of the Paket

    maxvalue = 1
;   
   For cnt = 0 To 1024
       If (maxvalue < outputarray(cnt))
           maxvalue = outputarray(cnt)
       EndIf
   Next cnt
 
 StartDrawing(WindowOutput(FFTWnd))
 Box(0,0,500,500,$0)
  MaxPeak=0
  For cnt = 5 To 512
    DiffY=Outputarray(cnt)/MaxValue*400
    LineXY(cnt,400,cnt,400-DiffY,$FFFFFF)
    If DiffY>MaxPeak
      MaxPeak=DiffY
      AkNote=cnt
    EndIf
  Next cnt
  SetWindowTitle(FFTWnd,Str(record_FindNote(aknote))+" on:"+Str(aknote)+" note:"+ShowNote_Get(record_FindNote(aknote)))
 
  StopDrawing()

EndProcedure

Procedure record_CallBack(hWnd.l,Msg.l,wParam.l,lParam.l)
 
  Result.l=#PB_ProcessPureBasicEvents
 
  Select Msg
    Case #WM_TIMER    : record_doFFT(Config\LScope)
    Case #MM_WIM_DATA : record_Read(wParam,lParam)
  EndSelect
 
  ProcedureReturn Result
 
EndProcedure



FFTWnd = OpenWindow(#PB_Any,0,0,500,500,"",#PB_Window_ScreenCentered | #PB_Window_SystemMenu)

Config\hWindow=WindowID(FFTWnd)
Config\output=WindowOutput(FFTWnd)

SetWindowCallback(@record_CallBack())

Record_Start()

Repeat:Until WaitWindowEvent() = #PB_Event_CloseWindow
End




DataSection
Notes:
        Data.w     14,  17
        Data.w     18,  18
        Data.w     19,  19
        Data.w     20,  20
        Data.w     21,  21
        Data.w     22,  23
        Data.w     24,  24
        Data.w     25,  26
        Data.w     27,  27
        Data.w     28,  29
        Data.w     30,  31
        Data.w     32,  32
       
        Data.w     33,  34
        Data.w     35,  37
        Data.w     38,  39
        Data.w     40,  41
        Data.w     42,  44
        Data.w     45,  46
        Data.w     47,  49
        Data.w     50,  52
        Data.w     53,  55
        Data.w     56,  59
        Data.w     60,  62
        Data.w     63,  66
                       
        Data.w     67,  70
        Data.w     71,  74
        Data.w     75,  79
        Data.w     80,  83
        Data.w     84,  88
        Data.w     89,  94
        Data.w     95,  99
        Data.w    100, 105
        Data.w    106, 112
        Data.w    113, 118
        Data.w    119, 125
        Data.w    126, 133
                       
        Data.w    134, 141
        Data.w    142, 149
        Data.w    150, 158
        Data.w    159, 168
        Data.w    169, 178
        Data.w    179, 188
        Data.w    189, 200
        Data.w    201, 212
        Data.w    213, 224
        Data.w    225, 238
        Data.w    239, 252
        Data.w    253, 267
                       
        Data.w    268, 283
        Data.w    284, 300
        Data.w    301, 318
        Data.w    319, 337
        Data.w    338, 357
        Data.w    358, 375
RealNotes:
        Data.s                     "C0", "C#0", "D0", "D#0", "E0", "F0", "F#0", "G0", "G#0"
        Data.s  "A0", "A#0", "B0", "C1", "C#1", "D1", "D#1", "E1", "F1", "F#1", "G1", "G#1"
        Data.s  "A2", "A#2", "B2", "C2", "C#2", "D2", "D#2", "E2", "F2", "F#2", "G2", "G#2"
        Data.s  "A3", "A#3", "B3", "C3", "C#3", "D3", "D#3", "E3", "F3", "F#3", "G3", "G#3"
        Data.s  "A4", "A#4", "B4", "C4", "C#4", "D4", "D#4", "E4", "F4"
EndDataSection

"FFT v0.1c" is the latest speed optimized version procedure i got managed to work. ASM optimisations are not impented yet.

Code: [Select]


; // *************************************************************************************
; // *
; // * Project:   FFT (THE FAST FOURIER TRANSFORM) v0.1c   
; // *
; // * Original version:  logix4u, www.logix4.net
; // * Coverted version:  Tranquil and Mr.Vain (Thorsten Will) in 2008
; // * Optimisations by:  va!n aka Thorsten Will
; // *
; // * Special thanks 2:  MrMat, Jim, Hellhound
; // *
; // * Reference:
; // * http://logix4u.net/DSP/Fast_Fourier_Transform/Visual_Basic_program_for_Fast_Fourier_Transform.html
; // *
; // * Upon entry, N.w contains the number of points in the DFT, REX() and
; // * IMX() contain the real and imaginary parts of the input. Upon return,
; // * REX() and IMX() contain the DFT output. All signals run from 0 to N.w-1
; // *
; // *************************************************************************************

Procedure record_doFFT( *scope.SCOPE )

    Define.d  TR, TI, SR, SI, UR, UI
    Define.l  J, K, L, M, N, NM1, ND2, cnt
    Define.L  MaxPeak, MaxValue, DiffY 
    Define.w  value

    ; // -------- Init some values for FFT analysing --------

    N   = 1024                ; // Number of samples
    NM1 = N - 1
    ND2 = N >> 1              ; // Optmimized, instead N / 2
    M   = 10                  ; // If N.w = 1024 == 10 == Same as:  Int(Log(N) / 0.69314718055994529)
    J   = ND2

    If Config\buffer = 0 : ProcedureReturn : EndIf

    ; // -------- Clear and Fill array values for analysing in just only one loop --------
         
    For i = 0 To Config\size Step 2       ; // Optimized by merging clear and fill array in one loop
        rex( i >> 1          ) = 0         ; //   0 to  512
        imx( i >> 1          ) = 0         ; //   0 to  512
        rex( i >> 1 + N >> 1 ) = 0         ; // 513 to 1024
        imx( i >> 1 + N >> 1 ) = 0         ; // 513 to 1024
 
        value    = PeekW( Config\buffer + i )
;       value.w    = PeekW( Config\buffer + i + *scope\channel * 2 )  ; // Enable this For Stereo Inpus
        rex( i >> 1 ) = value / 32767                                 ; // Optimized by doing i >> 1
;       imx( i >> 1 ) = 0                                             ; // Can be removed now
    Next

    ; // -------- Start FFT --------

    For i.l = 1 To N - 2                                ; // Bit reversal sorting
        If i < J
            TR = REX( J )
            TI = IMX( J )
            REX( J ) = REX( i )
            IMX( J ) = IMX( i )
            REX( i ) = TR
            IMX( i ) = TI
        EndIf
       
        K = ND2
       
        While K <= J
            J = J - K
            K = K >> 1                                  ; // Optmimized, instead N / 2
        Wend
        J = J + K
    Next

    For L = 1 To M                                      ; // Loop for each stage
        LE.l = 1 << L                                   ; // Optimized, instead  LE.l = Int( Pow( 2, L ) )
        LE2.l = LE >> 1                                 ; // Optimized, instead  N / 2
        UR = 1
        UI = 0
        SR =   Cos( #PI / LE2 )                         ; // Calculate sine & cosine values
        SI = - Sin( #PI / LE2 )
        For J.l = 1 To LE2                              ; // Loop for each sub DFT
            JM1.l = J - 1
            For i = JM1 To NM1                          ; // Loop for each butterfly
                IP.l = i + LE2
                TR = REX( IP ) * UR - IMX( IP ) * UI    ; // Butterfly calculation
                TI = REX( IP ) * UI + IMX( IP ) * UR
                REX( IP ) = REX( i ) - TR
                IMX( IP ) = IMX( i ) - TI
                REX( i  ) = REX( i ) + TR
                IMX( i  ) = IMX( i ) + TI
                i + LE - 1
            Next i
            TR = UR
            UR = TR * SR - UI * SI
            UI = TR * SI + UI * SR
        Next
    Next

    ; // -------- Calculate Outputarray and search for MaxValue of the Paket --------
 
    maxvalue = 0                                           ; // Optimized by merging calculate Outputarray
                                                           ; // and search MaxValue into just one loop.
    For cnt = 0 To N                                      ; // fixed to N.w instead wrong fixed value
        outputarray( cnt )  =  ( IMX( cnt ) * IMX( cnt ) )  +  ( REX( cnt ) * REX( cnt ) ) 
        If  maxvalue < outputarray( cnt )
            maxvalue = outputarray( cnt )
        EndIf
    Next
 
    ; // -------- Draw FFT -------- 
 
    StartDrawing( WindowOutput( FFTWnd ) )
        Box( 0, 0, N, 500, $0 )
        MaxPeak = 0
        For cnt = 0 To 500 ;                                   ; // Change fixed value to N.w !?
            DiffY = Outputarray( cnt ) / MaxValue * 400
            Box( cnt, 400, 1, -DiffY, $FFFFFF)
            If DiffY > MaxPeak
                MaxPeak = DiffY
                MaxPos  = cnt
            EndIf
        Next
    StopDrawing()

    SetWindowTitle( FFTWnd, ShowNote_Get( record_FindNote( MaxPos ) ) )

EndProcedure

ZIP attachement contains "FFT v0.1" (topic source == first version) and my latest "FFT v0.1c" versions with some speed optimisations (still not ASM optimized, yet)

[hellsangel]

great! 
it's what I want to do with a module tracker, but in C.  is it simple to adapt it in C for other player using wavout ?

[va!n]

@[Y]
Yes it should be easy to port to any other language... You just only have to put your output-sound-datas to the FFT input array to have the graph. good luck

[Jim]

Code: [Select]

        rex( i >> 1 + N >> 1 ) = 0         ; // 513 to 1024
        imx( i >> 1 + N >> 1 ) = 0         ; // 513 to 1024
[code]
In C, this would be equivalent to
(i>>(1+N))>>1
where I think you want
(i>>1)+(N>>1)
Might be worth putting the brackets in :)

Your routine will now be spending nearly all its time in Sin and Cos library routines...time to get the asm in ;)

Jim
[/code]