SyntaxBomb - Indie Coders

This Tutorial will show a FFT algorithm in Basic. Starting from the pure kernel I will develope a FFT-Class and a bundle of usefull functions around FFT in a module.

To keep the Tutorial clean, please write your questions and comments into this correponding thread:
https://www.syntaxbomb.com/blitzmax-blitzmax-ng/discussion-to-fft-module/msg347053648/#msg347053648




The pure FFT:
The following FFT-Code is an update from my FFT-algo from 2012. Now the code runs on BlitzMax NG and uses double precism.
Function FastFourier(Real:Double[] Var)
' F F T - F A S T   F O U R I E R   T R A N S F O R M A T I O N
'
' Author Midimaster www.midimaster.de
'
' Copyright: Public Domain
' Version 1.00
'
' The FFT needs an array REAL[] with the audio-samples
' The array needs a defined LENGTH
' LENGTH need to be 2^n (128, 256, 512, ... ,32768)
' As bigger LENGTH is as more precise is the result,
' but a long audio often contains often more than one tone
' so a good compromize is LENGTH=2^12, this means window is 85msec


Local Length:Int= Real.Length
Local Imag:Double[Length]
Local Middle:Int = Length/2

'  BIT REVERSING:
Local j:Int = Middle

For Local I:Int = 1 To Length - 2
If I < j Then
Local loc:Double = Real[j]
Real[j] = Real[I]
Real[i] = loc
' absolet, wenn Imag[]=0
loc     = Imag[j]
Imag[j] = Imag[I]
Imag[i] = loc
End If
Local K:Int  = Middle
While K <= j
j = j - K
K = K / 2
Wend
j = j + K
Next

'  FFT CALCULATION

Local Potenz:Int=1
For Local L:Int  = 1 To Int(Log(Length) / Log(2))   
Potenz:Int = Potenz Shl 1

Local Ur:Double = 1.0
Local Ui:Double = 0
Local arc:Double=  Double(180.0)/ Double(Potenz)
Local Sr:Double =  Cos(arc)
Local Si:Double = -Sin(arc)

For Local j:Int = 0 To Potenz-1
Local i:Int = j
While I < Length
Local Ip:Int    = I + Potenz
Local TR:Double = Real[ip] * Ur - Imag[ip] * Ui
Local TI:Double = Real[ip] * Ui + Imag[ip] * Ur
Real[ip] = Real[i] - TR
Imag[ip] = Imag[i] - TI
Real[i]  = Real[i] + TR
Imag[i]  = Imag[i] + TI
I        = I + 2*Potenz
Wend
Local loc:Double = Ur
Ur = loc * Sr - Ui * Si
Ui = loc * Si + Ui * Sr
Next
Next
End Function
 


And here is a complete runable BlitzMAX NG example:

Analyses a complex audiosignal build from three tones. Use: Move the mouse up/down to play around to variy audio frequencies.

' F F T - F A S T   F O U R I E R   T R A N S F O R M A T I O N
'
' Author Midimaster www.midimaster.de
'
' Copyright: Public Domain
' Version 0.99
'
'
'
' This demo shows:
' how to discover all frequencies in a complex audio sample   
' The FFT take a look (window) into the TAudioSample
' Size is the lenghth of the window (in samples) and need to be 2^n (128, 256, 512, ... ,32768)
' As bigger size is as more precise is the result, but a big size already contains often more than one tone
' so a good compromize is size=2^12, this means window is 85msec


SuperStrict
Graphics 800,600

Const SAMPLE_RATE:Int = 48000
Global  Size:Int=2^12

Global Music:Double[48000] ' 1 second of audio material in 48Khz
Global Real:Double[size]   ' holds the real part of the frequency domain
Global Imag:Double[size]   ' holds the imaginary part of the frequency domain
Global  Out:Double[size/2] ' holds the power results


'  move the MouseY up and down to change tone frequency
'  build your own audio in function CreateMusic()
' 
'
Global PeakAt:Double, Hertz:Int
Local LastMouse:Int=-1
MoveMouse 100,100
Flip
Repeat
Cls
If MouseY()<>LastMouse
LastMouse=MouseY()
Hertz=LastMouse^1.35
CreateMusic Hertz
FastFourier size, Real, Imag
EndIf
Paint
Flip 1
Until AppTerminate()


Function CreateMusic(Hertz:Int)
' fills an array with sample-values (64 bit floating point, values from -1 to +1)
Local f:Double = Hertz*360.0/SAMPLE_RATE

' first tone
For Local I:Int = 0 Until SAMPLE_RATE
Music[i]=Sin(f*i)
Next

'second tone: double frequency but half volume
For Local I:Int = 0 Until SAMPLE_RATE
Music[i]  =Music[i]+ Sin(2*f*i)/2
Next

'third tone: 4x frequency but 1/4 volume
For Local I:Int = 0 Until SAMPLE_RATE
Music[i] = Music[i]+ Sin(4*f*i)/4
Next
End Function





Function FastFourier(Length:Int , Real:Double[] Var ,  Imag:Double[] Var)
Local Zeit%=MilliSecs()

' STEP I: FETCH DATA
PrepareWindow Size

' STEP II: BIT REVERSING:

Local Middle:Int = Length/2
Local j:Int = Middle

For Local I:Int = 1 To Length - 2
If I < j Then
Local loc:Double = Real[j]
Real[j] = Real[I]
Real[i] = loc
End If
Local K:Int  = Middle
While K <= j
j = j - K
K = K / 2
Wend
j = j + K
Next

' STEP III: FFT CALCULATION

Local Potenz:Int=1
For Local L:Int  = 1 To Int(Log(Length) / Log(2))   
Potenz:Int = Potenz Shl 1

Local Ur:Double = 1.0
Local Ui:Double = 0
Local arc:Double=  Double(180.0)/ Double(Potenz)
Local Sr:Double =  Cos(arc)
Local Si:Double = -Sin(arc)

For Local j:Int = 0 To Potenz-1
Local i:Int = j
While I < Length
Local Ip:Int    = I + Potenz
Local TR:Double = Real[ip] * Ur - Imag[ip] * Ui
Local TI:Double = Real[ip] * Ui + Imag[ip] * Ur
Real[ip] = Real[i] - TR
Imag[ip] = Imag[i] - TI
Real[i]  = Real[i] + TR
Imag[i]  = Imag[i] + TI
I        = I + 2*Potenz
Wend
Local loc:Double = Ur
Ur = loc * Sr - Ui * Si
Ui = loc * Si + Ui * Sr
Next
Next

' STEP IV: GET RESULTS
CalculateResult
Print "Zeit = " + (MilliSecs()-zeit)

End Function


Function PrepareWindow(Quality:Int)
' fills the two FFT-arrays: one with the audio-window, the other is  reseted only
For Local I%=0 Until  size
Real[i]=Music[i]
Imag[i]=0
Next
End Function


Function CalculateResult()
' combines the Real and the Imaginary Array to a result
For Local i:Int=1 Until size/2
out[i] = Sqr((Imag[i] * Imag[i]) + (Real[i] * Real[i]))
Next

'find the highest value in the array
Local maxi:Double
For Local i%=1 Until size/2
If out[i] > maxi
maxi= out[i]
PeakAt=i
EndIf
Next
Print "PeakAt" + peakat + " " + maxi + " " + size
End Function



'  only graphics things for the demo app:

Function Paint()
SetColor 88,88,88
DrawRect 5,5,790,140
SetColor 0,255,0
For Local i%=1 Until 800
DrawLine  i-1, 80+Int(-Music[i-1]*30),i, 80+Int(-Music[i]*30)
Next
SetColor 255,255,0
For Local i%=1 Until size/2
DrawRect  i, 500,1, Int(-out[i]/3)
Next
SetColor 111,111,255
Local last%
DrawRect 0,500,800,1
For Local i%=1 Until size/2
If i*SAMPLE_RATE/size > last+480
last=last+500
Local t$= last/1000 +"kHz"
DrawText  "|", i, 500
If last Mod 1000=0 DrawText  t, i-TextWidth(t)/2, 520
EndIf
Next
SetColor 255,255,255
DrawText "20msec Audio:     3 Tones:   " + Hertz + "Hz   +"+ (2*Hertz) + "Hz    + "+ (4*Hertz) + "Hz", 10,10
DrawText "Length of FFT-Window = " + Int(size*1000/SAMPLE_RATE) + "msec", 50,560

DrawText "FFT found  Frequenz = " + Int(PeakAt*SAMPLE_RATE/size) + "Hz", 550,560
DrawText "Resolution = " + (SAMPLE_RATE/size) + "Hz", 550,580
End Function


What is FFT good for?

With FFT you can analyse audio files or audio signal, for involved frequencies. This does not work only with "simple" audios containing a single frequency, but also complex music recording. We can use this for finding out notes that are played. Or we can tune an instrument by displaying its frequency.

FFT breaks down one audio into many frequency bands. As there is also a IFFT (inverse FFT) you can recalculate frequency bands back to audio signal. This gives you the chance to manipulate a certain frequency band of a audio signal. Lets say: increase trebles or remove a singer's voice. Digital Filters

As a third you can move the frequency bands before you recalculte them with IFFT. This enables auto correction of tune of single notes or the whole song. Also time strechting is possible.


How does it work?

FFT is a fast Radix2 variant of the Diskrete-Fourier-Transformation. This examines a small period of time (called "window") and offers all frequencies that were involved during this period.

To calculate this in the fast way of FFT the window length needs to be a multiple of 2^n: 2, 4, 8, 16, 32, 64,....
The window length also defines the numbers of result frequecy bands. The bands are linear. In a 48kHz recording with a Window Length of 2048 you would get 2048 Bands of 23Hz width:
48000Hz / 2048 = 23.43Hz
better would be to work with 4096:
48000Hz / 4096 = 11.71Hz

This Window length would mean a audio length of 42.6msec
1000msec *2048 / 48000 = 42.6msec
better would be to work with 1024:
1000msec *1024 / 48000 = 21.3msec

This means as bigger the window length as more precise is the frequency information, but the audio length is increasing. This harbors the danger of accidentialy analysing more than one note of the melody into the same result.

The "quality" of the frequency information differs for the frequencies you are searching.

A typical female singer around 400Hz.

• 
  Difference between two halftones in this area: 20Hz
• 
  A 23Hz Bandwidth is nearly the same as 20Hz
• 
  This would mean exact enough to assign the FFT-Result to a true note
  

A typical Bass-note is around 100Hz

• 
  Difference between two halftones in this area: 5Hz
• 
  A 23Hz Bandwidth is more than 4 times bigger than 5Hz
• 
  This would mean impossible to assign the FFT-Result to one of the 4 notes
  

Nyquist Limitations

A window length of 2048 delivers 2048 frequency bands, but you can only use the half of them. A 48kHz Sample can only contain valid information of frequencies below 24kHz. As the Bandwidth is 23Hz, all result above band 1025 contain nonsense:
1024 x 23.43Hz = 24000


Window Limitations

All Fourier Transformations expect a Window with "perfect builded" audio-waves. This means (in a ideal case)  that all waves started at the beginning of the window at 0 and end extact at the end of the window.

In the real world you cutted the window from anywhere in your audio and the waves (at left and right border) were trunced. This would cause artifacts (wrong results). A solution We try to fade-in and fade-out the volume of the window, so that the audiosamples have value 0 at beginning and end of the window. There are several methods to do this from HAMM to GAUß. But this results in new artifacts.
[/code]



How to select Window Length and Sample Rate?

Both depents on the frequencies you want to examine. If your maximal frequency you expect is around 1000Hz, it makes no sense to sample with more that 3000Hz. With a Window Lenght of 3000 you would get very good 1Hz result bands, but this needs a long constant tone of 1000msec!!! With a Window Lenghth of 300 you would need a tone length of only 100msec, what corresponds to a 1/8-note. Now you would receive 150 bands each 10Hz.  This is good enough for assign the result to a note ( halftone at 1000Hz = 50Hz), but bad for a instrument tuning, where you need 1-5Hz exactness)

FFT-Class in BlitzMax NG

I transformed the FFT-knowledge into a first try of a FFT-Class. It enables to use FFT in a simple way.

New functions are:


FFT.AnalyseTAudioSample()
Function  AnalyseTAudioSample(Audio:TAudioSample, Start:Int, Length:Int, CallBackName: Byte Ptr)

This is the universal MAIN function when you have a trad. TAudioSample in SF_MONO16LE. It expects a TAudioSample and the desired window size. At the end the function reports "ready" to a callback function.
Example:
FFT:FFT_Class = New FFT_Class
Audio:TAudioSample = LoadAudioSample("...."
FFT.AnalyseTAudioSample Audio, 0, 2048 , MyDoneFunction
....
Function MyDoneFunction()
local Out:Double[] = FFT.GetResultArray()
' do whatever you want to do
End Function



FFT.MainFrequency()
Function  MainFrequency:Int()

returns the loudest frequency found in audio
Example:
FFT.Analyse Audio, 2048 , MyDoneFunction
...
local Loudest:Int = FFT.MainFrequency()



FFT.GetListOfFrequencies()
Function  GetListOfFrequencies:TList( sorted:Int=1, NewSampleRate:Int=0)

returns a list of peak frequencies found by FFT
Example:




FFT.GetListOfTones()
GetListOfTones:TList(sorted:Int=1, NewSampleRate:Int=0)

Returns a list of Tones (as MIDI NOTE NUMBERS) found by FFT

Example:




FFT.MidiFrom()
Function MidiFrom:Int(Frequency:Double)

Calculates the MIDI number related to a frequency

Example:
local Midinote:Int= FFT.MidiFrom(442)




FFT.GetResultArray()
Function  GetResultArray:Double[](mode:Int=1)

returns the raw result array of the FFT
Example:




FFT.Analyse()
Function  AnalyseTAudioSample(Audio:TAudioSample, Start:Int, Length:Int, CallBackName: Byte Ptr)

This is a deeper MAIN function for valid Audio in a 64bit-"Double" array give FFT a audio adress and a window size at the end the function reports "ready" with a callback function.

Example:
Global Audio:Double[2048] = ....
FFT.Analyse Audio, 2048 , MyDoneFunction
....
Function MyDoneFunction()
local Out:Double[] = FFT.GetResultArray()
' do whatever you want to do
End Function



FFT.Resolution()
Resolution:Int()

Returns the width of one band of the analysis in Hz

Example:




FFT.WindowTime()
WindowTime:Int()

Returns the length of the examined audio in msec

Example:



FFT.RemoveBand()
Function  RemoveBand(FromHz:Int, ToHz: Int)
You can manipulate the result of the FFT. Afterward you can recreate the audio from this.
Here: removed frequencies in an audio material (still in development).

Example:
FFT.Analyse Audio, 2048 , MyDoneFunction
...
FFT.Remove(2000,3000)
NewAudio:Double=FFT.InverseFFT()



FFT.InverseFFT()
Function  InverseFFT:Double[]()
You can manipulate the result of the FFT. Afterward you can recreate the audio from this.
Here: Calculates back the audio signal from a given frequencies array (still in development).

Example:
FFT.Analyse Audio, 2048 , MyDoneFunction
...
FFT.Remove(2000,3000)
NewAudio:Double=FFT.InverseFFT()





FFT_CLASS.bmx


Type FFT
' F F T - F A S T   F O U R I E R   T R A N S F O R M A T I O N
'
' Author Midimaster www.midimaster.de
'
' Copyright: Public Domain
' Version-Info  0.99
'
' The FFT needs an array REAL[] with the audio-samples
' The array needs a defined LENGTH
' LENGTH need to be 2^n (128, 256, 512, ... ,32768)
' A bigger LENGTH brings more results,
' but a long audio often contains already more than one tone
' so a good compromize is LENGTH=2^12, this means FFT-window is 85msec

Global ListOfTones:TList, ListOfFrequencies:TList
Global Real:Double[] , Imag:Double[], Out:Double[]
Global WindowLength:Int
Global HERTZ:Int
Global AlreadyDone:Int =False
Global CallBack()

' **************************************************************************
' PUBLIC FUNCTIONS
' **************************************************************************

Function  AnalyseTAudioSample(Audio:TAudioSample, Start:Int, Length:Int, CallBackName: Byte Ptr)
' PUBLIC
' this is the universal MAIN function when you have a trad. TAudioSample in SF_MONO16LE
' expects a TAudioSample and the desired window size
' at the end the function reports "ready" to a callback function.
' example:
'           Audio:TAudioSample =LoadAudioSample("....
'           
'           FFT.AnalyseTAudioSample Audio, 0, 2048 , MyDoneFunction
'           ....
'           Function MyDoneFunction()
'                 local Out:Double[] = FFT.GetResultArray()
'                 ' do whatever you want to do
'           End Function
'
WindowLength = CalculateWindowSize(Length)
HERTZ=Audio.Hertz
Local locAudio:Double[]=ConvertToDouble(Audio, Start, WindowLength)
Analyse locAudio, WindowLength, CallBackName
End Function


Function Analyse(Audio:Double[], Length:Int, CallBackName: Byte Ptr)
' PUBLIC
' this is a deeper MAIN function for valid Audio in a 64bit-"Double" array
' give FFT a audio adress and a window size
' at the end the function reports "ready" with a callback function.
' example:
'           Global Audio:Double[2048] = ....
'           FFT.Analyse Audio, 2048 , MyDoneFunction
'           ....
'           Function MyDoneFunction()
'                 local Out:Double[] = FFT.GetResultArray()
'                 ' do whatever you want to do
'           End Function
'
CallBack=CallBackName
WindowLength = CalculateWindowSize(Length)
PrepareWindow Audio, WindowLength
Process
CallBack()
End Function


Function GetResultArray:Double[](mode:Int=1)
' PUBLIC
' returns the raw result array of the FFT
' MODE=1 return Power of REAL and IMAG
' MODE=2 return simply REAL
Print mode
If Mode=1
Return CalculatePower()
ElseIf Mode=3
Return Imag
Else
Return OnlyReal()
EndIf
End Function


Function GetListOfFrequencies:TList( sorted:Int=1, NewSampleRate:Int=0)
' PUBLIC
' returns a list of peak frequencies found by FFT
' SampleRate= Samplerate of the TaudioSample, if not already known
' SORTED=1 Sorted by Frequency
' SORTED=2 Sorted by Signalstrength (loudest first)
If  NewSampleRate<>0 HERTZ = NewSampleRate

Out = CalculatePower()
Local List:TList =New TList
RuntimeError   "function GetListOfFrequencies not ready"
Return List
End Function


Function GetListOfTones:TList(sorted:Int=1, NewSampleRate:Int=0)
' PUBLIC
' returns a list of Tones (as MIDI NOTE NUMBERS) found by FFT
' SampleRate= Samplerate of the TaudioSample, if not already known

' SORTED=1 Sorted by Frequency
' SORTED=2 Sorted by Signalstrength (loudest first)
Out = CalculatePower()
Local PeakList:TList= GetListOfFrequencies( sorted, NewSampleRate)
Local List:TList =New TList
RuntimeError   "function GetListOfTones not ready"
Return List
End Function


Function MainFrequency:Int()
Local da:Int, Maxi:Double
For Local i%=0 Until OUT.Length
If OUT[i]>Maxi
Maxi=OUT[i]
da=i
EndIf
Next
Return GetFrequency(Da)
End Function


Function Resolution:Int()
Return HERTZ/WindowLength
End Function


Function WindowTime:Int()
Return 1000*WindowLength/HERTZ
End Function


' **************************************************************************
' PRIVATE FUNCTIONS
' **************************************************************************

Function CalculateWindowSize:Int(AnySize:Int)
' PRIVATE
    ' calculates allowed power-of-2 array length for FFT
Local ValidSize:Int=1
Repeat
anysize   = anysize/2
ValidSize = ValidSize*2
Until AnySize<2
Return ValidSize
End Function


Function GetFrequency:Int(At:Int)
If HERTZ=0 RuntimeError "Need given Sample-Rate first"
If WindowLength=0 RuntimeError "Need FFT-Analyse first"
Return Int(At*HERTZ/WindowLength)
End Function


Function GetPeak:Double(At:Int)
If AlreadyDone=False RuntimeError "first calculate results!"
Return Out[At]
End Function


Function CalculatePower:Double[]()
' PRIVATE
' combines the Real and the Imaginary Array to a result
Out = New Double[WindowLength/2]
For Local i:Int=1 Until WindowLength/2
out[i] = Sqr((Imag[i] * Imag[i]) + (Real[i] * Real[i]))
Next
AlreadyDone=True
Return out
End Function


Function OnlyReal:Double[]()
' PRIVATE
' returns the Real Array as a result
Out = New Double[WindowLength/2]
For Local i:Int=1 Until WindowLength/2
out[i] = Real[i]
Next
AlreadyDone=True
Return out
End Function


Function Process()
' PRIVATE
' STEP: BIT REVERSING:
Local Length:Int = WindowLength
Local Middle:Int = Length/2
Local j:Int = Middle

For Local I:Int = 1 To Length - 2
If I < j Then
Local loc:Double = Real[j]
Real[j] = Real[I]
Real[i] = loc
' absolet, wenn Imag[]=0
loc     = Imag[j]
Imag[j] = Imag[I]
Imag[i] = loc
End If
Local K:Int  = Middle
While K <= j
j = j - K
K = K / 2
Wend
j = j + K
Next
' STEP: FFT CALCULATION
Local Potenz:Int=1
For Local L:Int  = 1 Until Int(Log(Length) / Log(2))   
Potenz:Int = Potenz Shl 1

Local Ur:Double = 1.0
Local Ui:Double = 0
Local arc:Double=  Double(180.0)/ Double(Potenz)
Local Sr:Double =  Cos(arc)
Local Si:Double = -Sin(arc)

For Local j:Int = 0 To Potenz-1
Local i:Int = j
While I < Length
Local Ip:Int    = I + Potenz
Local TR:Double = Real[ip] * Ur - Imag[ip] * Ui
Local TI:Double = Real[ip] * Ui + Imag[ip] * Ur
Real[ip] = Real[i] - TR
Imag[ip] = Imag[i] - TI
Real[i]  = Real[i] + TR
Imag[i]  = Imag[i] + TI
I        = I + 2*Potenz
Wend
Local loc:Double = Ur
Ur = loc * Sr - Ui * Si
Ui = loc * Si + Ui * Sr
Next
Next

For Local i%=0 Until length/2
' Print i + "-->REAL=" + Int(real[i]) + "  --->IMAG="  + Int(imag[i])
Next
End Function



Function PrepareWindow(Audio:Double[], size:Int)
'PRIVATE
' fills the two FFT-arrays: one with the audio-window, the other is  reseted only
AlreadyDone=False
Real = New Double[size]
Imag = New Double[size]
For Local I%=0 Until  size
Real[i] = Audio[i]
Imag[i] = 0
Next
End Function


Function ConvertToDouble:Double[](Audio: TAudioSample, StartAt:Int, size:Int)
'PRIVATE
' brings trad. BlitzMax TAudioSample into Double-FlotingPoint-Format
Local locAudio:Double[size]
For Local I%=0 Until  size
Local s:Double=ShortToInt(Audio.Samples[StartAt + i])
locAudio[i] = S/32000.0
Next
Return locAudio
End Function

Function ShortToInt:Int( s:Int )
Return (s Shl 16) Sar 16
End Function


Function MidiFrom:Int(Frequency:Double)
Local NoteC:Double    = 16.3515625
Local NoteStep:Double =  1.0594630961707261
Local HalfStep:Double =  1.0293022375234235

Local Oktave%=1
Local tone%=0
Repeat
oktave:+1
NoteC:*2.0
Until NoteC*2>Frequency
Repeat
tone:+1
NoteC:*NoteStep
Until NoteC*HalfStep>Frequency
Return (oktave*12+tone)
End Function


Function InverseFFT:Double[]()
Local Teiler:Double = 1.0/ Double(WindowLength)

For Local i%=0 Until WindowLength
'FFT_Imag[i]: *-1
Imag[i]:*-1
'Imag[i]=0
Next

'FFTAnalyse_II WindowLength, FFT_Real, FFT_IMag, Imag, Real
'Analyse Real,  WindowLength, CallBackDone
process
For Local i:Int=0 Until WindowLength
Real[i]:*Teiler
' Imag[i]:= *-1
' Real[i]= Real[i]
Next
Return Real
End Function


Function FreqToBand:Int(Hz:Int)
For Local i%=0 Until WindowLength
If GetFrequency(i)>Hz
Return i
EndIf
Next
Return WindowLength-1
End Function


Function RemoveBand(FromHz:Int, ToHz: Int)
Local Start:Int= FreqToBand(FromHz)
Local Stop:Int= FreqToBand(ToHz)

For Local i%=start To  Stop
Real[i]=0
Imag[i]=0
Next
End Function


Function MoveBands(f:Double)
If f>1 Return
For Local da%=WindowLength To 0 Step -1
Real[da]=Real[da*f]
Imag[da]=Imag[da*f]
Next
End Function

End Type


runnable FFT-Calling example :

' F F T - F A S T   F O U R I E R   T R A N S F O R M A T I O N
'
' Author Midimaster www.midimaster.de
'
' Copyright: Public Domain
' Version 0.99
'
'
'
' This demo shows:
' how to discover all frequencies in a complex audio sample   
' The FFT take a look (window) into the TAudioSample
' Size is the lenghth of the window (in samples) and need to be 2^n (128, 256, 512, ... ,32768)
' As bigger size is as more precise is the result, but a big size already contains often more than one tone
' so a good compromize is size=2^12, this means window is 85msec


SuperStrict
Include "FFT_Class.bmx"
Graphics 800,600

Const SAMPLE_RATE:Int = 48000
Global  Size:Int=2^12

Global Music:Double[size] ' 1 second of audio material in 48Khz
Global OUT:Double[]
FFT.HERTZ=SAMPLE_RATE


'  move the MouseY up and down to change tone frequency
'  build your own audio in function CreateMusic()
' 
'
Global PeakAt:Double, Hertz:Int
Local LastMouse:Int=-1
MoveMouse 100,100
Flip

Global RetReal:Double[size]
Repeat
Cls
If MouseY()<>LastMouse
LastMouse=MouseY()
Hertz=LastMouse^1.35
CreateMusic Hertz
FFT.Analyse Music, size, CallBackDone
'FFT.removeBand 3000,9999
FFT.moveBands 0.5
'CallBackDone()
RetReal= FFT.InverseFFT()
EndIf
Paint
Flip 1
Until AppTerminate()


Function CallBackDone()
Out = FFT.GetResultArray(1)
PeakAt=FFT.MainFrequency()
End Function



Function CreateMusic(Hertz:Int)
' fills an array with sample-values (64 bit floating point, values from -1 to +1)
Local f:Double = Hertz*360.0/SAMPLE_RATE

' first tone
For Local I:Int = 0 Until Music.Length
Music[i]=Sin(f*i)
Next

'second tone: double frequency but half volume
For Local I:Int = 0 Until Music.Length
Music[i]  =Music[i]+ Sin(2*f*i)/2
Next

'third tone: 4x frequency but 1/4 volume
For Local I:Int = 0 Until Music.Length
Music[i] = Music[i]+ Sin(4*f*i)/4
Next
End Function




'  only graphics things for the demo app:

Function Paint()
SetColor 88,88,88
DrawRect 5,5,790,140
SetColor 0,255,0
For Local i%=1 Until Min(800,Out.Length)
DrawLine  i-1, 80+Int(-Music[i-1]*30),i, 80+Int(-Music[i]*30)
Next
SetColor 0,255,0
For Local i%=2 Until Min(800,Out.Length) Step 2
DrawLine  i-1, 280+Int(-retReal[i-2]*150),i, 280+Int(-retReal[i]*150)
Next
SetColor 255,255,0
For Local i%=1 Until OUT.length
DrawRect  i, 500,1, Int(-out[i])
Next
SetColor 111,111,255
Local last%
DrawRect 0,500,800,1
For Local i%=1 Until OUT.length
If i*SAMPLE_RATE/size > last+480
last=last+500
Local t$= last/1000 +"kHz"
DrawText  "|", i, 500
If last Mod 1000=0 DrawText  t, i-TextWidth(t)/2, 520
EndIf
Next
SetColor 255,255,255
DrawText "20msec Audio:     3 Tones:   " + Hertz + "Hz   +"+ (2*Hertz) + "Hz    + "+ (4*Hertz) + "Hz", 10,10
DrawText "Length of FFT-Window = " + FFT.windowtime() + "msec", 50,560
DrawText "FFT Main-Frequency = " + Int(PeakAt) + "Hz", 550,560
DrawText "         MIDI-Note = " + FFT.MidiFrom(PeakAt) + "", 550,580
DrawText "        Resolution = " + (FFT.Resolution()) + "Hz", 50,580
End Function


I played around with other Algorithm related to Frequency Anlysis and here is the...

Goertzel Algorithm

Goertzel also dectects a frequency in audio-samples. The difference to FFT is, that you can search for only one frequency band instead of receiving all audio bands as result. Another advantage is that there is no need to select a 2^n window size, but you can choose the window size as you like. So Goertzel looks like beeing faster as FFT.

But there alre also disadvantages: The result is a (single) value between 0 and some billions. And because you dont have the values of the other frequecy bands you can badly decide how "good" is this result. If you need more than a dozend frequencies to examine FFT is faster. FFT offers always thousands of bands in one run.

In my experiments I could see, that Goertzel delivers "very good" results at bass notes around 110Hz only when the window is longer than 500msec. The higher the note, the better are te results. For melody notes around 440Hz a windows size of 125msec is necessary.

Here the FFT allowes more conclusions, because you can better evaluate a single result by comparing it with all the neighbor band results you received in the same run.


The Algorithm

Function ExamineGoertzel:Double(Samples:Double[], SampleFreqency:Int, TargetFreq:Int)

• 
  Samples is an array of samples values. The length of the array defines the quality of Goertzel.
  You can convert a TAudioSample simply by moving each sample value to one cell of the array.
• 
  SampleFrequency is the HERTZ of the TAudioSample
• 
  TargetFreq is the frequency you want to detect

Function ExamineGoertzel:Double(Samples:Double[], SampleFreqency:Int, TargetFreq:Int)
Local OverLast:Double=0
Local Last:Double=0
Local N:Double = Double(Samples.Length)
Local k:Int            = 0.5 + N * TargetFreq / SampleFreqency
Local W:Double         = k/N *360
Local Cosinus:Double   = Cos(w)
Local Sinus:Double     = Sin(w)
Local coeff:Double     = 2*cosinus

For Local i%=0 Until N
Local Result:Double   = coeff*Last - OverLast + Samples[i]
OverLast = Last
Last     = Result
Next

Local real:Double      = Last- OverLast*cosinus
Local imag:Double      = Last*Sinus
Local Magnitude:Double
Magnitude              = (real^2 + imag^2)
'Magnitude             = real^2   +   imag^2   - ( real * imag *coeff )
Return  Max(0, Magnitude)
End Function


Here is a first testing app. (I will update it several times during the next days...)

Testing App

The app creates an Audio-Sample with some noise and some  other frequencies:

• 
  One is one halftone below 110Hz
• 
  one is one halftone above 110Hz
• 
  one is one octave below 110Hz
• 
  one is one octave above 110Hz
  

You can add your 110Hz signal (bass note) by pressing the mouse button. The graph will show the audio signal and the result of Goertzel. When moving the mouse horizontally you can change the size of the audio window and so observe the change of quality of Goertzel.

Graphics 800,600
Local ExamineFrequence:Int =110 'Hz
Const HERTZ:Int=3000
Global AudioLength=3000
Global Samples:Double[]
Local zeit%=MilliSecs()

Repeat
Cls
' some background noise and other frequenccies:
Samples = BuildNoiseAudio(AudioLength)
Samples = AddFrequency(Samples, HERTZ, 1/1.06*ExamineFrequence, 0.5)
Samples = AddFrequency(Samples, HERTZ, 1.06*ExamineFrequence, 0.5)
Samples = AddFrequency(Samples, HERTZ, ExamineFrequence/2, 0.5)
Samples = AddFrequency(Samples, HERTZ, ExamineFrequence*2, 0.5)

' change the examined samplesize with the mouse:
AudioLength = MouseX()*3+400
If MouseDown(1)
' now add your own frequency
Samples = AddFrequency(Samples, HERTZ, ExamineFrequence,0.5)
EndIf


'paint the design:
SetColor 0,255,255
DrawText " Sample-Rate: " + HERTZ + " Hz", 100 , 200
DrawText "Audio-Length: " + AudioLength*1000/HERTZ + " msec", 100 , 260
DrawText " Audio-Noise: " + Int(ExamineFrequence/2) + "Hz + " + Int(ExamineFrequence/1.06) + "Hz + " + Int(ExamineFrequence*1.06) + "Hz + " + Int(ExamineFrequence*2) + "Hz + WHITE NOISE", 100 , 290
Local LeftFrequency:Int = ExamineFrequence-15
Local Xadd:Int=20
If MouseDown(1)
DrawText "   Test-Tone: " + ExamineFrequence + " Hz", 100 , 230
SetColor 0,255,0
DrawText ExamineFrequence, (ExamineFrequence-LeftFrequency)*Xadd-8 , 530
Else
DrawText "   Test-Tone: --- Hz", 100 , 230
EndIf
DrawText Int(ExamineFrequence*1.06), (ExamineFrequence*1.06-LeftFrequency)*Xadd-8 , 530
DrawText Int(ExamineFrequence/1.06), (ExamineFrequence/1.06-LeftFrequency)*Xadd-8 , 530
SetColor 255,255,255
DrawRect 30,500,750,1
For Local i% = 0 To 100
If i Mod 5 = 0
DrawText i+LeftFrequency,i*Xadd-8,515
DrawRect i*30+4,500,2,14
EndIf
DrawRect i*30+4,500,2,7
Next 

'paint the waveform:
SetColor 0,255,0
For Local i%=1 To AudioLength
DrawLine i*2,Samples[i-1]*30+100,i*2+2,Samples[i]*Xadd+100
Next

'paint the goertzel result:
For Local i%=0 To 30
Local v:Double  = ExamineGoertzel(Samples,HERTZ,i+ExamineFrequence-15)
SetColor 0,255,255
DrawRect i*Xadd,500,10,-v/1000 +1
SetColor 255,255,255
' DrawText (i+ExamineFrequence-15),i*30-40,520
Next
Flip
Until AppTerminate()


Function ExamineGoertzel:Double(Samples:Double[], SampleFreqency:Int, TargetFreq:Int)
Local OverLast:Double=0
Local Last:Double=0
Local N:Double = Double(Samples.Length)
Local k:Int            = 0.5 + N * TargetFreq / SampleFreqency
Local W:Double         = k/N *360
Local Cosinus:Double   = Cos(w)
Local Sinus:Double     = Sin(w)
Local coeff:Double     = 2*cosinus

For Local i%=0 Until N
Local Result:Double   = coeff*Last - OverLast + Samples[i]
OverLast = Last
Last     = Result
Next

Local real:Double      = Last- OverLast*cosinus
Local imag:Double      = Last*Sinus
Local Magnitude:Double
Magnitude              = (real^2 + imag^2)
'Magnitude = real^2   +   imag^2   - ( real * imag *coeff )
Return  Max(0, Magnitude)
Return magnitude
End Function


Function AddFrequency:Double[](samples:Double[], SampleRate:Int, Frequency:Int, Volume:Double=1)
Local f:Double = (360.0*Frequency)/SampleRate

For Local i=0 Until samples.length
samples[i] = samples[i]+ Sin(i*f)*Volume
Next
Return samples
End Function 

Function BuildNoiseAudio:Double[](Length:Int)
Local samples:Double[Length]
For Local i=0 Until samples.length
samples[i] = Rnd(-0.3, + 0.3)
Next
Return samples
End Function



The app shows that Goertzel is not as good as FFT to find out notes in an audio file. It is more usefull in cases where you already know the tone you will hear, and you only want to find out if this tone is on or off. Goertzel is extremly resistent agains noise. -noise does not influence the result, so you could us it to listen to a certain frequency in a loud surrounding. Or find out a frequency in a very weak signal. (f.e. WLAN or lang distance wires, etc..)
 

I read a lot about audio manipulation this time and try to understand what happens. Now here is a first tutorial about the IIR-filtering or Biquad-Filtering:

Biquad Filter

To filter frequency bands from audio it is not neccessary to use FFT or manipulate the complete audio as a block. With a method called Biquad Filtering manipulations are possible by only have a look on three neighbor sample values.

Let us call them :

• 
  Sample: the current moment sample (n)
• 
  Last: the sample before the current (n-1)
• 
  OverLast: the sample 2 steps prior (n-2)
  

Also we will need these results of our calculations:

• 
  Result: the current moment sample after we "transformed" it. (n)
• 
  LastResult: the result of the last transformation (n-1)
• 
  OverLastResult: the result from 2 steps prior (n-2)
  

With only have a look on these  6 sample values at the same time, you can calculate all "Filter-jobs"  with audio signals: High-Pass, Low-Pass, Notch,  EQ, Delay, Wha-wha, Chorus, etc...


Only One Code Line For Filtering

The Result is always calculated by an addition of the other 5. Each of the 5 elements needs to be multiplicated by an individual constant value.
Result =(Sample*B_0) + (Last*B_1) + (OverLast*B_2) -(LastResult*A_1) - (OverLastResult*A_2)


Defining The Filter Means Only Defining 5 Constants

The 5 constants B_0 , B_1 , B_2 , A_1 , A_2 are pre-calculated when selecting the filter type.

This is the "genius" part of filter design. Here is an example for a Low-Pass-Filter that deletes high frequencies and let through low frequencies:

Low-Pass Filter

Function DefineLowPass(LimitFrequency:Int, SampleRate:Int)
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
ResetAllGlobals
Local w:Double = 360.0*LimitFrequency/SampleRate
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

' here are now the 5 contants:
B_0 = Float ((1-COSINUS)*0.5 * INV_ALPHA)
B_1 = Float ((1-COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Function

This function only needs the SampleRate of your TAudio and the start frequency from where it will cut out the higher frequencies. It is only need once to pre-calculate the 5 constants.

Here is a runnable example:

SuperStrict
Global Sample:TAudioSample =LoadAudioSample("test.ogg")

TBiQuad.DefineLowPass (1500,Sample.Hertz)
Print "HERTZ = " + Sample.Hertz
Print "Length of audio = " + sample.Length
Global Pointer:Short Ptr = Sample.Samples

Local zeit%=MilliSecs()
For Local i%=0 Until sample.Length
Local value:Double = ShortToInt(Pointer[i])
value = TBiquad.ProcessFilter(value)
Pointer[i] = Int(Value)
Next
Print "used time = " + (MilliSecs()-zeit) + "msec"
Global Sound:TSound=LoadSound(Sample)
PlaySound Sound
WaitKey()


Type TBiquad
Global B_0:Double, B_1:Double, B_2:Double, A_1:Double, A_2:Double
Global OverLastResult:Double, LastResult:Double
Global Last:Double, OverLast:Double


Function DefineLowPass(LimitFrequency:Int, SampleRate:Int)
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
ResetAllGlobals
Local w:Double = 360.0*LimitFrequency/SampleRate
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float ((1-COSINUS)*0.5 * INV_ALPHA)
B_1 = Float ((1-COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Function


Function ProcessFilter:Double(Sample:Double)
Local Result:Double = (Sample*B_0) + (Last*B_1) + (OverLast*B_2) -(LastResult*A_1) - (OverLastResult*A_2)
OverLast       = Last
Last           = Sample
OverLastResult = LastResult
LastResult     = Result
Return Result
End Function


Function ResetAllGlobals()
B_0=0
B_1=0
B_2=0
A_1=0
A_2=0
OverLastResult=0
LastResult=0
Last=0
OverLast=0
End Function

End Type


Function ShortToInt:Int( s:Int )
    Return (s Shl 16) Sar 16
End Function


The test.ogg is here in attachment:

Next step was to create a HIGH-PASS filter.


High PassFilter

it uses nearly the same caclulations for the 5 constants:

Method DefineHighPass()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float ((1+COSINUS)*0.5 * INV_ALPHA)
B_1 = Float ((-1-COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method



The TFilter Class

Now I packed everything into a real class. This means you can now run more than one filter on the same audio material. Current version of the class:

SuperStrict
Global Sample:TAudioSample =LoadAudioSample("test.ogg")
Print "HERTZ = " + Sample.Hertz
Print "Length of audio = " + sample.Length
Global Pointer:Short Ptr = Sample.Samples


Local LowPass:TFilter  = TFilter.CreateLowPass  (1500,Sample.Hertz)
Local HighPass:TFilter = TFilter.CreateHighPass (1000,Sample.Hertz)


Local zeit%=MilliSecs()
For Local i%=0 Until sample.Length
Local value:Double = ShortToInt(Pointer[i])
value = LowPass.Process(value)
value = HighPass.Process(value)
Pointer[i] = Int(Value)
Next
Print "used time = " + (MilliSecs()-zeit) + "msec"

Global Sound:TSound=LoadSound(Sample)
PlaySound Sound
WaitKey()



Type TFilter
Const   HIGH:Int=1, LOW:Int=0
Field   Typ:Int, Limit:Int, HERTZ:Int
Field   B_0:Double, B_1:Double, B_2:Double, A_1:Double, A_2:Double
Field OverLastResult:Double, LastResult:Double
Field Last:Double, OverLast:Double

'**********************************************************************************
' PUBLIC

' LOW PASS FILTER:
Function CreateLowPass:TFilter( LimitFrequency:Int, SampleRate:Int)
Local loc:TFilter = New TFilter
loc.Typ   = LOW
loc.Limit = LimitFrequency
loc.HERTZ = SampleRate
loc.DefineLowPass
Return loc
End Function


' HIGH PASS FILTER:
Function CreateHighPass:TFilter( LimitFrequency:Int, SampleRate:Int)
Local loc:TFilter = New TFilter
loc.Typ   = LOW
loc.Limit = LimitFrequency
loc.HERTZ = SampleRate
loc.DefineHighPass
Return loc
End Function


' PROCESS FOR ALL FILTERS:
Method Process:Double(Sample:Double)
Local Result:Double = (Sample*B_0) + (Last*B_1) + (OverLast*B_2) -(LastResult*A_1) - (OverLastResult*A_2)
OverLast       = Last
Last           = Sample
OverLastResult = LastResult
LastResult     = Result
Return Result
End Method


'**********************************************************************************
' PPRIVATE

Method DefineLowPass()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float ((1-COSINUS)*0.5 * INV_ALPHA)
B_1 = Float ((1-COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method



Method DefineHighPass()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float ((1+COSINUS)*0.5 * INV_ALPHA)
B_1 = Float ((-1-COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method



'**********************************************************************************
' EXPERIMENTAL

Method OptimalProcess:Double(Sample:Double)
' 20% faster for LowPass and HighPass
Local Result:Double = ((Sample+OverLast)*B_0) + (Last*B_1)  -(LastResult*A_1) - (OverLastResult*A_2)
OverLast       = Last
Last           = Sample
OverLastResult = LastResult
LastResult     = Result
Return Result
End Method


Method ResetAllGlobals()
B_0=0
B_1=0
B_2=0
A_1=0
A_2=0
OverLastResult=0
LastResult=0
Last=0
OverLast=0
End Method

End Type


Function ShortToInt:Int( s:Int )
    Return (s Shl 16) Sar 16
End Function



The audiofile test.ogg is avaiable as attachment in post#4: https://www.syntaxbomb.com/tutorials/fft-fast-fourier-transformation-in-blitzmax-ng/?action=dlattach;attach=5167

Ok... here are two more filter:

BandPass Filter

Only the frequencies closed the band frequency keep remain

Method DefineBandPass()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float (ALPHA * INV_ALPHA)
B_1 = 0
B_2 = -B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method



Notch Filter

The frequencies around the filter frequency will be cut away from the audio

Method DefineNotch()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = INV_ALPHA
B_1 = Float((-2*COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = B_1
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method


The last version of TFilter-Class:

SuperStrict
Global Sample:TAudioSample =LoadAudioSample("test.ogg")
Print "HERTZ = " + Sample.Hertz
Print "Length of audio = " + sample.Length
Global Pointer:Short Ptr = Sample.Samples


Local LowPass:TFilter  = TFilter.CreateLowPass  (200,Sample.Hertz)
Local HighPass:TFilter = TFilter.CreateHighPass (1000,Sample.Hertz)
Local BandPass:TFilter  = TFilter.CreateBandPass  (1500,Sample.Hertz)
Local Notch:TFilter  = TFilter.CreateNotch  (1500,Sample.Hertz)


Local zeit%=MilliSecs()
For Local i%=0 Until sample.Length
Local value:Double = ShortToInt(Pointer[i])
value = Notch.Process(value)
Pointer[i] = Int(Value)
Next
Print "used time = " + (MilliSecs()-zeit) + "msec"

Global Sound:TSound=LoadSound(Sample)
PlaySound Sound
WaitKey()



Type TFilter
Const   HIGH:Int=1, LOW:Int=0, BAND:Int=2
Field   Typ:Int, Limit:Int, HERTZ:Int
Field   B_0:Double, B_1:Double, B_2:Double, A_1:Double, A_2:Double
Field OverLastResult:Double, LastResult:Double
Field Last:Double, OverLast:Double

'**********************************************************************************
' PUBLIC

' LOW PASS FILTER:
Function CreateLowPass:TFilter( LimitFrequency:Int, SampleRate:Int)
Local loc:TFilter = New TFilter
loc.Typ   = LOW
loc.Limit = LimitFrequency
loc.HERTZ = SampleRate
loc.DefineLowPass
Return loc
End Function


' HIGH PASS FILTER:
Function CreateHighPass:TFilter( LimitFrequency:Int, SampleRate:Int)
Local loc:TFilter = New TFilter
loc.Typ   = LOW
loc.Limit = LimitFrequency
loc.HERTZ = SampleRate
loc.DefineHighPass
Return loc
End Function


' BAND PASS FILTER:
Function CreateBandPass:TFilter( LimitFrequency:Int, SampleRate:Int)
Local loc:TFilter = New TFilter
loc.Typ   = BAND
loc.Limit = LimitFrequency
loc.HERTZ = SampleRate
loc.DefineBandPass
Return loc
End Function


' NOTCH FILTER:
Function CreateNotch:TFilter( LimitFrequency:Int, SampleRate:Int)
Local loc:TFilter = New TFilter
loc.Typ   = BAND
loc.Limit = LimitFrequency
loc.HERTZ = SampleRate
loc.DefineNotch
Return loc
End Function


' PROCESS FOR ALL FILTERS:
Method Process:Double(Sample:Double)
Local Result:Double = (Sample*B_0) + (Last*B_1) + (OverLast*B_2) -(LastResult*A_1) - (OverLastResult*A_2)
OverLast       = Last
Last           = Sample
OverLastResult = LastResult
LastResult     = Result
Return Result
End Method


'**********************************************************************************
' PPRIVATE

Method DefineLowPass()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float ((1-COSINUS)*0.5 * INV_ALPHA)
B_1 = Float ((1-COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method



Method DefineHighPass()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float ((1+COSINUS)*0.5 * INV_ALPHA)
B_1 = Float ((-1-COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method



Method DefineBandPass()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = Float (ALPHA * INV_ALPHA)
B_1 = 0
B_2 = -B_0
A_1 = Float((-2*COSINUS)     * INV_ALPHA)
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method



Method DefineNotch()
Local COSINUS:Double, ALPHA:Double, INV_ALPHA:Double
Local QUALITY:Double = 0.7
Local w:Double = 360.0*Limit/HERTZ
COSINUS  = Cos(w)
ALPHA    = Sin(w)/2/QUALITY
INV_ALPHA = 1/(1+ALPHA)

B_0 = INV_ALPHA
B_1 = Float((-2*COSINUS)     * INV_ALPHA)
B_2 = B_0
A_1 = B_1
A_2 = Float ((1-ALPHA)       * INV_ALPHA)
End Method


'**********************************************************************************
' EXPERIMENTAL

Method OptimalProcess:Double(Sample:Double)
' 20% faster for LowPass and HighPass
Local Result:Double = ((Sample+OverLast)*B_0) + (Last*B_1)  -(LastResult*A_1) - (OverLastResult*A_2)
OverLast       = Last
Last           = Sample
OverLastResult = LastResult
LastResult     = Result
Return Result
End Method


Method OptimalBandProcess:Double(Sample:Double)
' 20% faster for BandPass
Local Result:Double = ((Sample-OverLast)*B_0) - (LastResult*A_1) - (OverLastResult*A_2)
OverLast       = Last
Last           = Sample
OverLastResult = LastResult
LastResult     = Result
Return Result
End Method


Method OptimalNotchProcess:Double(Sample:Double)
' 20% faster for Notch
Local Result:Double = (Sample+OverLast)*B_0 + (Last-LastResult)*B_1  - (OverLastResult*A_2)
OverLast       = Last
Last           = Sample
OverLastResult = LastResult
LastResult     = Result
Return Result
End Method

Method ResetAllGlobals()
B_0=0
B_1=0
B_2=0
A_1=0
A_2=0
OverLastResult=0
LastResult=0
Last=0
OverLast=0
End Method

End Type


Function ShortToInt:Int( s:Int )
    Return (s Shl 16) Sar 16
End Function


The audiofile test.ogg is avaiable as attachment in post#4: https://www.syntaxbomb.com/tutorials/fft-fast-fourier-transformation-in-blitzmax-ng/?action=dlattach;attach=5167

I see this is ancient, but I've just discovered it. I play a lot with synthesis and recording nowadays, so it's really cool to see the code behind this and to have nice working examples. Many thanks!

@Midimaster
That's a nice piece of work there.