Realtime FFT Analyzer

Share your advanced PureBasic knowledge/code with the community.
Tranquil
Addict
Addict
Posts: 950
Joined: Mon Apr 28, 2003 2:22 pm
Location: Europe

Realtime FFT Analyzer

Post by Tranquil »

Hi all.

As some ppl on the forum has searched for a solution to analyze raw sound datas frequency spectrum, we converted a FFT Routine found on the web.

The following source uses your default soundcard's input and displays its spectrum in a little graph on a window just for testing purpose.

We adjusted the settings in this code to determine a frequency which is at most peak at the moment up to 8 kHz. If you need a higher frequency you only need to adjust the input-frequency. The determined most powerfull frequency is displayed in the Window-Title as the Note that has been found by the FFT. For testing this, I added an archive with some sinus tones. You can download them here: (~9MB Archive)

http://www.danceyourmusic.com/tones.rar

For more accuracy we use 1024 samples in our buffer. If you need lower letancy you can select a lower value (256 oder 512) but it must be a power of two.

Thanks to MrMat for his help on converting the C code.

If you have any hint how to improve Accuracy or speed, just let us know. :)

Have fun.

Code: Select all

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
Tranquil
User avatar
einander
Enthusiast
Enthusiast
Posts: 744
Joined: Thu Jun 26, 2003 2:09 am
Location: Spain (Galicia)

Post by einander »

Very good work!
Thanks for sharing
Dr_Wildrick
User
User
Posts: 36
Joined: Fri Feb 23, 2007 8:00 pm
Location: New York

Removed

Post by Dr_Wildrick »

Removed
Last edited by Dr_Wildrick on Sun Feb 24, 2008 6:05 am, edited 1 time in total.
ABBKlaus
Addict
Addict
Posts: 1143
Joined: Sat Apr 10, 2004 1:20 pm
Location: Germany

Post by ABBKlaus »

Dr_Wildrick wrote:#WAVE_FORMAT_PCM <--- Does not seem to be defined
You should run it with PB4.10 and above :wink:

(i think you are running PB4.02)
hellhound66
Enthusiast
Enthusiast
Posts: 119
Joined: Tue Feb 21, 2006 12:37 pm

Post by hellhound66 »

Removed.
Last edited by hellhound66 on Wed Mar 19, 2008 11:24 pm, edited 1 time in total.
Dr_Wildrick
User
User
Posts: 36
Joined: Fri Feb 23, 2007 8:00 pm
Location: New York

Sine language

Post by Dr_Wildrick »

I updated, Pain in the butt, bit it works now. I noticed that at the center of your frequency power display iit always spike off the top No matter what I sent through it. Anyone elase have similar results? Or do I just need a new sound card. LOL
_Dr. W
va!n
Addict
Addict
Posts: 1104
Joined: Wed Apr 20, 2005 12:48 pm

Post by va!n »

i am very happy tranquil and me managed with help of MrMat (who helped us to fix a conversion bug) a working realtime FFT. Here is a very small optimized and bugfixed version. (Even the Sin and Cos things can be optimized a lot with asm too, not implented yet)

Code: Select all

; // *************************************************************************************
; // *
; // * Project:   FFT (THE FAST FOURIER TRANSFORM) v0.1a    
; // * 
; // * Original version:  logix4u, www.logix4.net 
; // * Coverted version:  Tranquil and Mr.Vain (Thorsten Will) in 2008
; // *
; // * 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) / Log(2) )
    J.l   = ND2

    If Config\buffer = 0 : ProcedureReturn : EndIf

    ; // -------- Clear array values --------

    For pos = 0 To N.w          ; fixed to N.w instead fixed value
        rex( pos ) = 0
        imx( pos ) = 0 
    Next

    ; // -------- Fill array values for analysing --------
    
    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( i >> 1 ) = value / 32767                                 ; // Optimized by doing i >> 1
        imx( i >> 1 ) = 0  
    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 aka Thorsten

Intel i7-980X Extreme Edition, 12 GB DDR3, Radeon 5870 2GB, Windows7 x64,
va!n
Addict
Addict
Posts: 1104
Joined: Wed Apr 20, 2005 12:48 pm

Post by va!n »

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 all

; // *************************************************************************************
; // *
; // * Project:   FFT (THE FAST FOURIER TRANSFORM) v0.1b    
; // * 
; // * Original version:  logix4u, www.logix4.net 
; // * Coverted version:  Tranquil and Mr.Vain (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 aka Thorsten

Intel i7-980X Extreme Edition, 12 GB DDR3, Radeon 5870 2GB, Windows7 x64,
hellhound66
Enthusiast
Enthusiast
Posts: 119
Joined: Tue Feb 21, 2006 12:37 pm

Post by hellhound66 »

Code: Select all

LE.l = Int( Pow( 2, L ) ) 
This should be the same as

Code: Select all

LE.l = 1<<L
Tranquil
Addict
Addict
Posts: 950
Joined: Mon Apr 28, 2003 2:22 pm
Location: Europe

Re: Sine language

Post by Tranquil »

Dr_Wildrick wrote:I updated, Pain in the butt, bit it works now. I noticed that at the center of your frequency power display iit always spike off the top No matter what I sent through it. Anyone elase have similar results? Or do I just need a new sound card. LOL
_Dr. W
What Input source do you use? The code should display the frequency which has the most power as the max. value of the graph. (This is what we need in our projekt).

@hellhound66:

You are right about the bad coding style. We are just happy that it works after converting VERY old basic codes. The linke you mentioned (m-w) whas a typo on conversion. Thanks for spotting this out. :)
Tranquil
Dr_Wildrick
User
User
Posts: 36
Joined: Fri Feb 23, 2007 8:00 pm
Location: New York

Removed

Post by Dr_Wildrick »

Removed
Last edited by Dr_Wildrick on Sun Feb 24, 2008 6:06 am, edited 1 time in total.
va!n
Addict
Addict
Posts: 1104
Joined: Wed Apr 20, 2005 12:48 pm

Post by va!n »

@Hellhound66:
Yes, you are right. Thanks for the "LE.l = 1<<L" tipp. I didnt noticed this before. 8) *thumbs up*

Dr_Wildrick:
Thanks for your feedback! Version 0.1c is a little bit more speed optimized. About the slow drawings... There are ways to speed even the drawing stuff up a lot. But the main feature/function to demonstrate here is the FFT itself. But glad you like the optimized version and that it works now smooth on your machine =)

@all:
Just some very small speed improvements. Special thanks must go to hellhound for the "LE.l = 1<<L" tipp!

@Tranquil:
Please dont use the "Fast Sine & Cosine lookup tables" for the FFT, because they are to slow in comparsion to the inline ASM i will try to add very soon. Hope you like the new optimized versions too ;)

Code: Select all

; // *************************************************************************************
; // *
; // * 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
va!n aka Thorsten

Intel i7-980X Extreme Edition, 12 GB DDR3, Radeon 5870 2GB, Windows7 x64,
dell_jockey
Enthusiast
Enthusiast
Posts: 767
Joined: Sat Jan 24, 2004 6:56 pm

Post by dell_jockey »

Va!n,

I'm very much looking forward to your optimised ASM sin/cos lookup code. Will this be a modular thing that can be dropped into other apps easily? Too bad I'm not yet into ASM (just yet...).
cheers,
dell_jockey
________
http://blog.forex-trading-ideas.com
chris319
Enthusiast
Enthusiast
Posts: 782
Joined: Mon Oct 24, 2005 1:05 pm

Post by chris319 »

In looking over your FFT code I am wondering if this:

Code: Select all

    ; // -------- 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 ) )
should read like this:

Code: Select all

    ; // -------- 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 )  = Sqr(( IMX( cnt ) * IMX( cnt ) )  +  ( REX( cnt ) * REX( cnt )))
taking the square root of those values.
chris319
Enthusiast
Enthusiast
Posts: 782
Joined: Mon Oct 24, 2005 1:05 pm

Post by chris319 »

More optimizations:

Code: Select all

;DECLARE THESE OUTSIDE PROCEDURE
Global N.l   = 1024 ; // Number of samples
Global M.l   = 10 ; // If N.w = 1024 == 10 == Same as:  Int(Log(N) / 0.69314718055994529)

Code: Select all

Procedure record_doFFT( *scope.SCOPE )

    Define.d  TR, TI, SR, SI, UR, UI
    Define.l  J, K, L, 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
    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 ) ) ;NO NEED TO CALL THIS EVERY TIME
        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() ;NO NEED TO CALL THIS EVERY TIME

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

EndProcedure

Code: Select all

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())

StartDrawing( WindowOutput( FFTWnd ) )

Record_Start()

Repeat: Until WaitWindowEvent() = #PB_Event_CloseWindow

StopDrawing()

End
Post Reply