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.
(https://www.syntaxbomb.com/tutorials/fft-fast-fourier-transformation-in-blitzmax-ng/?action=dlattach;attach=5088)
' 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 LimitationsA 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 LimitationsAll 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 AlgorithmFunction 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 AppThe 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 FilterTo 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 FilteringThe
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 ConstantsThe 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!
@MidimasterThat's a nice piece of work there.