## FFT Fast Fourier Transformation in BlitzMax NG

Started by Midimaster, December 15, 2021, 09:10:50

#### Midimaster

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.

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.
Code (BlitzMax) Select
`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.

Code (BlitzMax) Select
`' 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 85msecSuperStrict Graphics 800,600Const SAMPLE_RATE:Int = 48000Global  Size:Int=2^12Global Music:Double[48000] ' 1 second of audio material in 48KhzGlobal Real:Double[size]   ' holds the real part of the frequency domainGlobal Imag:Double[size]   ' holds the imaginary part of the frequency domainGlobal  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:IntLocal LastMouse:Int=-1MoveMouse 100,100FlipRepeat Cls If MouseY()<>LastMouse LastMouse=MouseY() Hertz=LastMouse^1.35 CreateMusic Hertz FastFourier size, Real, Imag EndIf Paint Flip 1Until 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 NextEnd 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 FunctionFunction 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 + " " + sizeEnd 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,580End Function `
...back from Egypt

#### Midimaster

#1
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)

...back from Egypt

#### Midimaster

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:
Code (BlitzMax) Select
`FFT:FFT_Class = New FFT_ClassAudio: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:
Code (BlitzMax) Select
`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:
Code (BlitzMax) Select

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

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

Example:
Code (BlitzMax) Select

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

Calculates the MIDI number related to a frequency

Example:
Code (BlitzMax) Select
`local Midinote:Int= FFT.MidiFrom(442)`

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

returns the raw result array of the FFT
Example:
Code (BlitzMax) Select

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:
Code (BlitzMax) Select
`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:
Code (BlitzMax) Select

FFT.WindowTime()
WindowTime:Int()

Returns the length of the examined audio in msec

Example:
Code (BlitzMax) Select

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:
Code (BlitzMax) Select
`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:
Code (BlitzMax) Select
`FFT.Analyse Audio, 2048 , MyDoneFunction...FFT.Remove(2000,3000)NewAudio:Double=FFT.InverseFFT()`

FFT_CLASS.bmx
Code (BlitzMax) Select
`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 :

Code (BlitzMax) Select
`' 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 85msecSuperStrict Include "FFT_Class.bmx"Graphics 800,600Const SAMPLE_RATE:Int = 48000Global  Size:Int=2^12Global Music:Double[size] ' 1 second of audio material in 48KhzGlobal 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:IntLocal LastMouse:Int=-1MoveMouse 100,100FlipGlobal 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 1Until 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 NextEnd 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,580End Function `
...back from Egypt

#### Midimaster

#3
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

Code (BlitzMax) Select
`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.

Code (BlitzMax) Select
`Graphics 800,600Local ExamineFrequence:Int =110 'HzConst HERTZ:Int=3000Global AudioLength=3000Global 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 FlipUntil 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 magnitudeEnd 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 samplesEnd 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 samplesEnd 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..)

...back from Egypt

#### Midimaster

#4
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:

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

Code (BlitzMax) Select
` 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:

Code (BlitzMax) Select
`SuperStrictGlobal Sample:TAudioSample =LoadAudioSample("test.ogg")TBiQuad.DefineLowPass (1500,Sample.Hertz)Print "HERTZ = " + Sample.HertzPrint "Length of audio = " + sample.LengthGlobal 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 SoundWaitKey()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 TypeFunction ShortToInt:Int( s:Int )    Return (s Shl 16) Sar 16End Function`

The test.ogg is here in attachment:
...back from Egypt

#### Midimaster

#5
Next step was to create a HIGH-PASS filter.

High PassFilter

it uses nearly the same caclulations for the 5 constants:

Code (BlitzMax) Select
` 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:

Code (BlitzMax) Select
`SuperStrictGlobal Sample:TAudioSample =LoadAudioSample("test.ogg")Print "HERTZ = " + Sample.HertzPrint "Length of audio = " + sample.LengthGlobal 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) NextPrint "used time = " + (MilliSecs()-zeit) + "msec"Global Sound:TSound=LoadSound(Sample)PlaySound SoundWaitKey()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 TypeFunction ShortToInt:Int( s:Int )    Return (s Shl 16) Sar 16End 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
...back from Egypt

#### Midimaster

Ok... here are two more filter:

BandPass Filter

Only the frequencies closed the band frequency keep remain

Code (BlitzMax) Select
` 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

Code (BlitzMax) Select
` 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:

Code (BlitzMax) Select
`SuperStrictGlobal Sample:TAudioSample =LoadAudioSample("test.ogg")Print "HERTZ = " + Sample.HertzPrint "Length of audio = " + sample.LengthGlobal 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) NextPrint "used time = " + (MilliSecs()-zeit) + "msec"Global Sound:TSound=LoadSound(Sample)PlaySound SoundWaitKey()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 TypeFunction ShortToInt:Int( s:Int )    Return (s Shl 16) Sar 16End 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
...back from Egypt

#### DruggedBunny

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!

#### col

@Midimaster
That's a nice piece of work there.
https://github.com/davecamp

"When you observe the world through social media, you lose your faith in it."