Realtime FFT Analyzer

[eddy]

converted in PB 4.5+

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

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

Structure NoteRange
   Note.i
   FromPos.i
   ToPos.i
EndStructure

Procedure.i ShowNote_Init()
   
   Global Dim NoteRange.NoteRange(53)
   
   For Note=0 To 53
      Read.w FromPos.w
      Read.w 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.s 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.i
   top.i
   width.i
   height.i
   middleY.i
   quarterY.i
EndStructure

Structure CONFIG
   
   hWindow.i           ; Window handle
   size.i              ; Wave buffer size
   buffer.i            ; Wave buffer pointer
   output.i            ; WindowOutput()
   wave.i              ; Address of waveform-audio input device
   
   format.WAVEFORMATEX ; Capturing WaveFormatEx
   lBuf.i              ; Capturing Buffer size
   nBuf.i              ; Capturing Buffer number
   nDev.i              ; Capturing Device identifier
   nBit.i              ; Capturing Resolution (8/16)
   nHertz.i            ; Capturing Frequency  (Hertz)
   nChannel.i          ; 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.i, lpWaveHdr.i)
   *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)
   
   Define.d TR, TI, SR, SI, UR, UI
   Define.i J, K, L, NM1, ND2, cnt
   Define.i 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.i=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.i=1<<L                                    ; // Optimized, instead  LE.i = Int( Pow( 2, L ) )
      LE2.i=LE>>1                                  ; // Optimized, instead  N / 2
      UR=1
      UI=0
      SR=Cos(#PI/LE2)                              ; // Calculate sine & cosine values
      SI=-Sin(#PI/LE2)
      For J.i=1 To LE2                             ; // Loop for each sub DFT
         JM1.i=J-1
         For i=JM1 To NM1                          ; // Loop for each butterfly
            IP.i=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

Procedure record_CallBack(hWnd.i, Msg.i, wParam.i, lParam.i)
   
   Result.i=#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())

StartDrawing(WindowOutput(FFTWnd))
   Record_Start()   
   Repeat : Until WaitWindowEvent()=#PB_Event_CloseWindow   
StopDrawing()
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

[Booger]

This seems to be very inaccurate in vista 64 for some reason or another. I.e. Analyzer reports a C1 when Fmod reports a D2.

[chris319]

Added some constants and made the code more efficient. Now displays approximate frequency in a text gadget.

Code: Select all

#N = 1024 ;NUMBER OF SAMPLES
#NM1 = #N - 1
#ND2 = #N >> 1
#M = 10
#WINDOW_WIDTH = 530
;#WINDOW_WIDTH = #N / 2 ;512
#WINDOW_HEIGHT = 400
#SAMPLE_RATE = 40000 ;44100 ;8000 ;LIMIT 4000 Hz AUDIO
#TRACE_COLOR = #White
#gadText1 = 1

Structure MYWAVEFORMATEX
  wFormatTag.u ;UNSIGNED FOR COMPATIBILITY WITH MICROSOFT FLAG VALUES
  nChannels.w
  nSamplesPerSec.l
  nAvgBytesPerSec.l
  nBlockAlign.w
  wBitsPerSample.w
  cbSize.w
EndStructure

Global Dim rex.f(#N + 1)
Global Dim imx.f(#N + 1)
Global Dim OutPutArray.f(#N + 1)
Global FFTWnd

Structure NoteRange
  Note.i
  FromPos.i
  ToPos.i
EndStructure

Procedure.i ShowNote_Init()
           
           Global Dim NoteRange.NoteRange(53)
           
           For Note=0 To 53
              Read.w FromPos.w
              Read.w 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.s 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.i
           top.i
           width.i
           height.i
           middleY.i
           quarterY.i
        EndStructure

        Structure CONFIG
           hWindow.i           ; Window handle
           size.i              ; Wave buffer size
           buffer.i            ; Wave buffer pointer
           output.i            ; WindowOutput()
           wave.i              ; Address of waveform-audio input device
           
           format.MYWAVEFORMATEX ; Capturing WaveFormatEx
           lBuf.i              ; Capturing Buffer size
           nBuf.i              ; Capturing Buffer number
           nDev.i              ; Capturing Device identifier
           nBit.i              ; Capturing Resolution (8/16)
           nHertz.i            ; Capturing Frequency  (Hertz)
           nChannel.i          ; Capturing Channels number (Mono/Stereo)
           
           LScope.SCOPE        ; Waveform display
           RScope.SCOPE        ; Waveform 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 = #SAMPLE_RATE
           Config\nDev = 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.i, lpWaveHdr.i)
           *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)
       
       Define.d TR, TI, SR, SI, UR, UI
       Define.i J, K, L, cnt, MaxValue
       Define.w value
       
       ; // -------- Init some values for FFT analysing --------

        ;THESE ARE NOW CONSTANTS AS THEY DON'T CHANGE       
       ;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
          MessageRequester("Error", "No buffer available.")
          Goto shutdown
        EndIf
       
       ; // -------- Clear and Fill array values for analysing in just only one loop --------
      buff = Config\buffer       
       For i = 0 To Config\size Step 2 ; // Optimized by merging clear and fill array in one loop
          i2 = i >> 1
          rex(i2) = 0             ; 0 to  512
          imx(i2) = 0             ; 0 to  512
          rex(i2 + #N >> 1) = 0        ; 513 to 1024
          imx(i2 + #N >> 1) = 0        ; 513 to 1024
         
          value.w = PeekW(buff + i)
          ;value.w    = PeekW( Config\buffer + i + *scope\channel * 2 )  ; // Enable this For Stereo Inpus
          rex(i2) = value / 32767   ; // Optimized by doing i >> 1
       Next
       
       ; // -------- Start FFT --------
       
       For i.i = 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.i = 1 << L                                    ; // Optimized, instead  LE.i = Int( Pow( 2, L ) )
          LE2.i = LE >> 1                                  ; // Optimized, instead  N / 2
          UR = 1
          UI = 0
          SR = Cos(#PI / LE2)                              ; // Calculate sine & cosine values
          SI = -Sin(#PI / LE2)

          For J.i = 1 To LE2                             ; // Loop for each sub DFT
             JM1.i = J - 1

             For i = JM1 To #NM1                          ; // Loop for each butterfly
                IP.i = 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 ;1024
              outputarray(cnt) = (IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt))
              If maxvalue < outputarray(cnt)
                 maxvalue = outputarray(cnt)
                 freq = cnt * (#SAMPLE_RATE / 1000)
              EndIf
           Next
           
           ; // -------- Draw FFT --------
           
       Box(0, 0, #WINDOW_WIDTH, #WINDOW_HEIGHT + 2, #Black)

       For cnt = 0 To #WINDOW_WIDTH - 1
          yCoord = (outputArray(cnt) / maxValue) * -#WINDOW_HEIGHT
          If yCoord < 0
            Line(cnt, 400, 1, yCoord, #TRACE_COLOR)
          EndIf
       Next
           
       SetWindowTitle(FFTWnd, ShowNote_Get(record_FindNote(MaxPos)))
       SetGadgetText(#gadText1, "Frequency: " + Str(freq) + " Hz")       
           
EndProcedure

Procedure record_CallBack(hWnd.i, Msg.i, wParam.i, lParam.i)
  If Msg = #MM_WIM_DATA
    record_Read(wParam, lParam): record_doFFT(Config\LScope)
  EndIf

ProcedureReturn #PB_ProcessPureBasicEvents
EndProcedure

FFTWnd = OpenWindow(#PB_Any, 0, 0, #WINDOW_WIDTH, #WINDOW_HEIGHT + 100, "", #PB_Window_ScreenCentered | #PB_Window_SystemMenu)
Config\hWindow=WindowID(FFTWnd)
Config\output=WindowOutput(FFTWnd)

SetWindowCallback(@record_CallBack())

StartDrawing(WindowOutput(FFTWnd))
TextGadget(#gadText1, 10, 410, 110, 20, "Frequency: ") ;DISPLAYS APPROXIMATE FREQUENCY

Record_Start()   
Repeat : Until WaitWindowEvent() = #PB_Event_CloseWindow   

shutdown:
StopDrawing()

For i = 0 To Config\nBuf - 1
  FreeMemory(inHdr(i)\lpData) ;FREE ALLOCATED MEMORY
Next

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

[Simo_na]

@chris319
How i can separate the left and right channels data to draw in two canvas, one for left channel and one for right channel ?

[doctornash]

Some mods:


1. Included a choice of window functions - Hann, Hamming, Blackman-Harris. These reduce spectral leakage (suppressing sidelobes in other bins compared to the main lobe) making the display more 'stable' and the result better suited for frequency analysis. Have left the Hann as default, as it is a good general-purpose window. The Hamming window is better at cancelling the nearest side lobe but worse at cancelling any others.

2. It would seem it is supposed to be:

Code: Select all

outputarray(cnt) = Sqr((IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt)))

not

Code: Select all

outputarray(cnt) = (IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt))

As evidenced in the picture. The reference scope image on the left is of the ' Soundcard Oscilloscope' https://www.zeitnitz.eu/scope_en which is showing the spectrum of a 423 Hz Audacity-generated sawtooth wave, for which the harmonic heights are supposed to be 1/2, then 1/3, then 1/4 etc of the fundamental. With Sqr, the PB spectrum on the right aligns well with the reference.

3. Re the maxvalue Frequency display, in our case where frequency of the waveform is 423 Hz:

Code: Select all

freq = cnt * (#SAMPLE_RATE / 1000)

gives freq = 440 Hz, which is way off.

Code: Select all

freq = cnt * (#SAMPLE_RATE / 1024)

gives freq = 429 Hz, which is closer.

For even better resolution, if we do

Code: Select all

Debug Str(cnt * (#SAMPLE_RATE / 1024)) + ";" + StrD(outputarray(cnt))

we see these as the bins on either side:
390 ; 8.4267539978
429 ; 9.4770374298
468 ; 7.7056999207
We can use Quadratic interpolation between the 3 points to estimate the max frequency, per:
y1 = |X[k - 1]|
y2 = |X[k]|
y3 = |X[k + 1]|
d = (y3 - y1) / (2 * (2 * y2 - y1 - y3))
k' = k + d

So:
y1 = 8.4267539978
y2 = 9.4770374298
y3 = 7.7056999207
d = -0.721/(2*(18.954-8.426754-7.7057))) = -0.721/(2*2.82162) = -0.12776
The bin size is 429-390 = 39
-0.12776*39 = -5
So, Max freq = 429-5 = 424 Hz, and this result aligns well with the 423.67 Hz the reference scope shows for it

Code: Select all

#N = 1024 ;NUMBER OF SAMPLES
#NM1 = #N - 1
#ND2 = #N >> 1
#M = 10
#WINDOW_WIDTH = 530
;#WINDOW_WIDTH = #N / 2 ;512
#WINDOW_HEIGHT = 400
#SAMPLE_RATE = 40000 ;44100 ;8000 ;LIMIT 4000 Hz AUDIO
#TRACE_COLOR = #White
#gadText1 = 1

Structure MYWAVEFORMATEX
  wFormatTag.u ;UNSIGNED FOR COMPATIBILITY WITH MICROSOFT FLAG VALUES
  nChannels.w
  nSamplesPerSec.l
  nAvgBytesPerSec.l
  nBlockAlign.w
  wBitsPerSample.w
  cbSize.w
EndStructure

Global Dim rex.f(#N + 1)
Global Dim imx.f(#N + 1)
Global Dim OutPutArray.f(#N + 1)
Global FFTWnd

Global Pi2.f = 6.28318530717958

Structure NoteRange
  Note.i
  FromPos.i
  ToPos.i
EndStructure

Procedure.i ShowNote_Init()
           
           Global Dim NoteRange.NoteRange(53)
           
           For Note=0 To 53
              Read.w FromPos.w
              Read.w 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.s 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.i
           top.i
           width.i
           height.i
           middleY.i
           quarterY.i
        EndStructure

        Structure CONFIG
           hWindow.i           ; Window handle
           size.i              ; Wave buffer size
           buffer.i            ; Wave buffer pointer
           output.i            ; WindowOutput()
           wave.i              ; Address of waveform-audio input device
           
           format.MYWAVEFORMATEX ; Capturing WaveFormatEx
           lBuf.i              ; Capturing Buffer size
           nBuf.i              ; Capturing Buffer number
           nDev.i              ; Capturing Device identifier
           nBit.i              ; Capturing Resolution (8/16)
           nHertz.i            ; Capturing Frequency  (Hertz)
           nChannel.i          ; Capturing Channels number (Mono/Stereo)
           
           LScope.SCOPE        ; Waveform display
           RScope.SCOPE        ; Waveform 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 = #SAMPLE_RATE
           Config\nDev = 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.i, lpWaveHdr.i)
           *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)
       
       Define.d TR, TI, SR, SI, UR, UI
       Define.i J, K, L, cnt, MaxValue
       Define.w value
       
       Define MaxIndex.i
       
       ; // -------- Init some values for FFT analysing --------

        ;THESE ARE NOW CONSTANTS AS THEY DON'T CHANGE       
       ;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
          MessageRequester("Error", "No buffer available.")
          ;Goto shutdown
        EndIf
       
       ; // -------- Clear and Fill array values for analysing in just only one loop --------
      buff = Config\buffer       
       For i = 0 To Config\size Step 2 ; // Optimized by merging clear and fill array in one loop
          i2 = i >> 1
          rex(i2) = 0             ; 0 to  512
          imx(i2) = 0             ; 0 to  512
          rex(i2 + #N >> 1) = 0        ; 513 to 1024
          imx(i2 + #N >> 1) = 0        ; 513 to 1024
         
          value.w = PeekW(buff + i)
          ;value.w    = PeekW( Config\buffer + i + *scope\channel * 2 )  ; // Enable this For Stereo Inpus
          ;rex(i2) = value / 32767   ; // Optimized by doing i >> 1 ;Rectangular window
          rex(i2) = value*(0.5*(1-Cos(Pi2*(i/(Config\size)))))/32767 ;Hann window
          ;rex(i2) = value*(0.53836 - 0.46164*Cos(Pi2*i/(Config\size)))/32767 ;Hamming window
          ;rex(i2) = value*(0.3635819-0.4891775*Cos(Pi2*i/(Config\size))+0.1365995*Cos(2*Pi2*i/(Config\size))-0.0106411*Cos(3*Pi2*i/(Config\size)))/32767 ;Blackman-Harris window            
          
       Next
       
       ; // -------- Start FFT --------
       
       For i.i = 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.i = 1 << L                                    ; // Optimized, instead  LE.i = Int( Pow( 2, L ) )
          LE2.i = LE >> 1                                  ; // Optimized, instead  N / 2
          UR = 1
          UI = 0
          SR = Cos(#PI / LE2)                              ; // Calculate sine & cosine values
          SI = -Sin(#PI / LE2)

          For J.i = 1 To LE2                             ; // Loop for each sub DFT
             JM1.i = J - 1

             For i = JM1 To #NM1                          ; // Loop for each butterfly
                IP.i = 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  
       MaxIndex = 0
       ; // Optimized by merging calculate Outputarray
           ; // and search MaxValue into just one loop.
           For cnt = 0 To #N ;1024
             outputarray(cnt) = Sqr((IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt)))
             ;outputarray(cnt) = (IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt)) ;this does not seem right
             ;Debug Str(cnt * (#SAMPLE_RATE / 1024)) + ";" + StrD(outputarray(cnt))
             ;Debug Str(cnt * (#SAMPLE_RATE / 1000)) + ";" + StrD(outputarray(cnt)) ;this does not seem right
              If maxvalue < outputarray(cnt)
                maxvalue = outputarray(cnt)
                MaxIndex = cnt
              EndIf
            Next
            
            If (MaxIndex>0) And (MaxIndex<#N-1)   
              y1.d = OutPutArray(MaxIndex-1)
              y2.d = OutPutArray(MaxIndex)
              y3.d = OutPutArray(MaxIndex+1)
              delta.d = (y3 - y1) / (2 * (2 * y2 - y1 - y3))          
              freq = (MaxIndex * (#SAMPLE_RATE / 1024)) + (delta*(#SAMPLE_RATE / 1024))            
            Else   
              freq = MaxIndex * (#SAMPLE_RATE / 1024)
            EndIf

            
           ; // -------- Draw FFT --------
            
       Box(0, 0, #WINDOW_WIDTH, #WINDOW_HEIGHT + 2, #Black)
       For cnt = 0 To #WINDOW_WIDTH - 1
          yCoord = (outputArray(cnt) / maxValue) * -#WINDOW_HEIGHT
          If yCoord < 0
            Line(cnt, 400, 1, yCoord, #TRACE_COLOR)
          EndIf
       Next
           
       SetWindowTitle(FFTWnd, ShowNote_Get(record_FindNote(MaxPos)))
       SetGadgetText(#gadText1, "Frequency: " + Str(freq) + " Hz") 
    
EndProcedure

Procedure record_CallBack(hWnd.i, Msg.i, wParam.i, lParam.i)
  If Msg = #MM_WIM_DATA
    record_Read(wParam, lParam): record_doFFT(Config\LScope)
  EndIf

ProcedureReturn #PB_ProcessPureBasicEvents
EndProcedure

FFTWnd = OpenWindow(#PB_Any, 0, 0, #WINDOW_WIDTH, #WINDOW_HEIGHT + 100, "", #PB_Window_ScreenCentered | #PB_Window_SystemMenu)
Config\hWindow=WindowID(FFTWnd)
Config\output=WindowOutput(FFTWnd)

SetWindowCallback(@record_CallBack())

StartDrawing(WindowOutput(FFTWnd))
TextGadget(#gadText1, 10, 410, 110, 20, "Frequency: ") ;DISPLAYS APPROXIMATE FREQUENCY

Record_Start()   
Repeat : Until WaitWindowEvent() = #PB_Event_CloseWindow   

shutdown:
StopDrawing()

For i = 0 To Config\nBuf - 1
  FreeMemory(inHdr(i)\lpData) ;FREE ALLOCATED MEMORY
Next

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

[infratec]

Hi,

if you use EnableExplicit you'll find that

Code: Select all

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

is wrong.
Maybe it should be

Code: Select all

SetWindowTitle(FFTWnd, ShowNote_Get(record_FindNote(MaxValue)))

Bernd

[Seymour Clufley]

Is there any way to get this to analyse a WAV file, rather than live audio?

[AlieLeite]

Hi All....I am trying to input audio into my soundcard, analyze it with a FFT (for a LED spectrum analyzer) and then output the audio through the soundcard in real time.Can anyone suggest software tools that can accomplish this? I am aware of gstreamer, but the python documentation is minimal or outdated. Does anyone know of decent alternatives or decent documentation?

assemblage PCB

[CELTIC88]

hi ,

@AlieLeite , @Seymour Clufley

is easy with "bass.dll"

look here:
viewtopic.php?f=12&t=70878

goood luck

[Simo_na]

hi ,

@AlieLeite , @Seymour Clufley

is easy with "bass.dll"

look here:
viewtopic.php?f=12&t=70878

goood luck

fft_left=getleftcodedevice0
fft_right=getrightcodedevice0

it would be easy even under Purebasic, wanting to do it....even a 16-bit 44100 ... would suffice