Hello everyone,
Although I have taken my leave from these forums a while ago, I still wanted to share this plugin with you.
Sometimes I use DBP for prototyping, and a week ago I had the need for an FFT algorithm in DBP (engineers like to work in the frequency domain). DBP doesn't have this (FFT's usually aren't needed for game development), so I searched the internet and found FFTW3 (
http://www.fftw.org/) and made a wrapper for it.
I realize that the audience for this is quite small, but I still think it would be a good plugin for DBP to have.
Here are the commands
memblock fft Memblock_from, Memblock_to, Points
memblock ifft Memblock_from, Memblock_to, Points
memblock fftr Memblock_from, Memblock_to, Points
memblock ifftr Memblock_from, Memblock_to, Points
fft [complex_ptr_in, complex_ptr_out, int_N, int_Sign, int_Flags][complex_ptr_in, complex_ptr_out, int_N]
fft_2d complex_ptr_in, complex_ptr_out, int_N0, int_N1, int_Sign, int_Flags
fft_3d complex_ptr_in, complex_ptr_out, int_N0, int_N1, int_N2, int_Sign, int_Flags
fft_c2r complex_ptr_in, float_ptr_out, int_N, int_Flags
fft_c2r_2d complex_ptr_in, float_ptr_out, int_N0, int_N1, int_Flags
fft_c2r_3d complex_ptr_in, float_ptr_out, int_N0, int_N1, int_N2, int_Flags
fft_r2c float_ptr_in, complex_ptr_out, int_N, int_Flags
fft_r2c_2d float_ptr_in, complex_ptr_out, int_N0, int_N1, int_Flags
fft_r2c_3d float_ptr_in, complex_ptr_out, int_N0, int_N1, int_N2, int_Flags
ifft complex_ptr_in, complex_ptr_out, int_N
fftr float_ptr_in, complex_ptr_out, int_N
ifftr complex_ptr_in, float_ptr_out, int_N
Notes
Attached is a zip file containing the DLL, keywords, constants and a PDF with information that comes with FFTW3. All commands use this prototype of C++ code:
fftwf_plan p = [create plan]
fftwf_execute(p);
fftwf_destroy_plan(p);
where [create plan] corresponds to the command name in question:
fft -> fftwf_plan_dft_1d()
fft_2d -> fftwf_plan_dft_2d()
fft_3d -> fftwf_plan_dft_3d()
fft_r2c -> fftwf_plan_dft_r2c_1d()
...
with a 1 on 1 mapping of the parameters.
fft/ifft are the same as fft with Flags and Sign equal to FFTW_ESTIMATE and FFTW_FORWARD(fft)/FFT_BACKWARD(ifft).
fftr/ifftr are the same as fft_r2c/fft_c2r with the flag FFTW_ESTIMATE.
To provide easier access, I also included some commands to work with memblocks.
memblock fft/ifft/...
These commands can check the memblock size and see if there is enough memory for the FFT. There is also no need for pointers this way. Downside is that it is a little slower, and that you can only compute the FFT from the start of the memblock.
Caution when using
When using the pointer versions, you must always double check that the size is large enough to contain the FFT result. If the allocated memory (memblock or 'make memory') is not large enough, then the program will CRASH. Unfortunately, this crash will often happen randomly and can become hard to trace because of that.
Examples
Record sound and display FFT (like an equalizer)
sync on
sync rate 0
#constant SAMPLE_DURATION 50
#constant FFT_SAMPLES 512
global SoundRecording = 0
global LastTimer = 0
global SampleRate as float
global RecordedSamples as integer
`Make a memblock with sufficient space for the samples + DC component
make memblock 2, FFT_SAMPLES * 4 + 8
do
cls
`Clear the screen
KeepRecording()
f# = DrawSound(2, 0, screen height() / 4 * 3, 10.0 / FFT_SAMPLES, FFT_SAMPLES / 2)
text 0, 0, str$(f#) + " Hz"
text 0, 20, "Sample rate: " + str$(SampleRate) + " Hz"
text 0, 40, "Samples: " + str$(RecordedSamples)
sync
loop
function KeepRecording()
if timer() - LastTimer > SAMPLE_DURATION
`Stop recording the sound
stop recording sound
`Convert to a memblock
if memblock exist(1) then delete memblock 1
if sound exist(1) then make memblock from sound 1, 1
`Start recording
record sound 1, SAMPLE_DURATION
LastTimer = timer()
if memblock exist(1)
`Read sampling rate
SampleRate = memblock dword(1, 8)
RecordedSamples = (get memblock size(1) - 28)
if RecordedSamples > FFT_SAMPLES
`Use a window to reduce leakage
for i = 0 to FFT_SAMPLES - 1
a# = 0.5 * (1.0 - cos(180.0 * i / (1.0 * FFT_SAMPLES)))
write memblock float 2, i*4, (memblock byte(1, 28 + i) - 128) * a#
next i
memblock fftr 2, 2, FFT_SAMPLES
endif
endif
endif
endfunction
function DrawSound(mem, x, y, scale#, Samples)
maxA# = 0.0
maxFreq# = 0.0
if memblock exist(mem)
ix = -1
r = make vector2(1)
for i = 0 to Samples - 1
`Get average of sound samples
set vector2 1, memblock float(mem, i*8), memblock float(mem, i*8 + 4)
SoundH# = length vector2(1) * scale#
if SoundH# > maxA#
maxA# = SoundH#
maxFreq# = 0.5 * SampleRate / Samples * i
endif
`Drawing
if ix >= 0 then line x + ix, y - oldSoundH#, x + ix + 2, y - SoundH#
inc ix, 2
oldSoundH# = soundH#
next i
endif
endfunction MaxFreq#
Run the program and try whistling or singing in your microphone. You can see how the frequency spectrum changes when you change your pitch height.
FFT example with time and frequency domain
N = 1024 : `Number of points
Fs# = 2.5e3 : `Sample frequency
Bin# = Fs# / N : `Frequency resolution
F1# = Fs# / 16.0 : `Excited sine wave
F2# = Fs# / 4.0
`Make a sampled sinewave (and pad some zeros for the fft)
make memblock 1, N*4 + 8
for i = 0 to N - 1
write memblock float 1, i*4, sin(360.0 * F1# * i / Fs#) + sin(360.0 * F2# * i / Fs# + 121.0)
next i
`TIME DOMAIN
cls
`Grid line
ink rgb(0, 255, 0), 0
line 0, 128, 1024, 128
ink rgb(255, 255, 255), 0
for i = 1 to N - 1
if i mod 100 = 0
ink rgb(0, 255, 0), 0
line i, 123, i, 138
center text i, 140, str$(i / Fs# * 1000, 0) + " ms"
ink rgb(255, 255, 255), 0
endif
line i - 1, 128.0 - (memblock float(1, (i-1)*4) * 64.0), i, 128.0 - (memblock float(1, i*4) * 64.0)
next i
text 0, 256, "Sample frequency: " + str$(Fs#, 1) + "kHz"
text 0, 276, "Excited sinewaves: " + str$(F1#, 1) + "Hz and " + str$(F2#, 1) + "Hz"
`FFT
memblock fftr 1, 1, N
`FREQUENCY DOMAIN
`Grid line
ink rgb(0, 255, 0), 0
line 0, 700, 1024, 700
ink rgb(255, 255, 255), 0
prevA# = sqrt(memblock float(1, 0)^2.0 + memblock float(1, 4)^2.0) / 2.0
for i = 1 to N/2
`Draw spectrum
A# = sqrt(memblock float(1, i*8)^2.0 + memblock float(1, i*8 + 4)^2.0) / 2.0
line i - 1, 700 - prevA#, i, 700 - A#
`Spectrum is symmetrical, so draw symmetric from other side
line 1024 - i + 1, 700 - prevA#, 1024 - i, 700 - A#
prevA# = A#
next i
`Show excited lines
ink rgb(255, 0, 0), 0
x = int(F1# / Bin# + 0.5)
line x, 695, x, 705
center text x, 710, str$(F1#, 1) + "Hz"
x = int(F2# / Bin# + 0.5)
line x, 695, x, 705
center text x, 710, str$(F2#, 1) + "Hz"
wait key
end
A demonstration on using FFT's. The total signal is a sum of two sine waves, and by using the FFT, you can see which frequency and phase these sine waves have.
Low pass filter on an image
load image "image.png", 1, 1
W = image width(1)
H = image height(1)
make memblock from image 1, 1
`Convert to float map
make memblock 2, W*H*4
ConvertToFloatMap(1, 2, 12, W, H)
ConvertToByteMap(2, 1, 12, W, H, 1.0)
make image from memblock 1, 1
`FFT
make memblock 3, (W+2)*H*4
fft_r2c_2d get memblock ptr(2), get memblock ptr(3), H, W, 1 << 6
`Convert to a heightmap
nW = W/2 + 1
nH = H/2
ConvertComplexToAmplitude(3, 2, nW, nH, 0.5 / sqrt(W * H))
`Make a new image for the fft
nW = W/2 + 1
nH = H/2
make memblock 4, 12 + (nW*nH*4)
write memblock dword 4, 0, nW
write memblock dword 4, 4, nH
write memblock dword 4, 8, 32
ConvertToByteMap(2, 4, 12, nW, nH, 1.0)
make image from memblock 2, 4
paste image 2, 0, 0
`Filter out the high frequency components (right and bottom area)
for y = 0 to nH - 1
for x = 0 to nW - 1
pos = (y*nW + x)*8
if x > nW/8 or y > nH/8
write memblock float 3, pos, 0.0
write memblock float 3, pos + 4.0, 0.0
endif
next x
next y
ConvertComplexToAmplitude(3, 2, nW, nH, 0.5 / sqrt(W * H))
ConvertToByteMap(2, 4, 12, nW, nH, 1.0)
make image from memblock 3, 4
paste image 3, 512, 0, 1
`Take inverse fft
fft_c2r_2d get memblock ptr(3), get memblock ptr(2), H, W, 1 << 6
`Convert to byte map
ConvertToByteMap(2, 1, 12, W, H, 1.0 / (H * W))
make image from memblock 4, 1
paste image 1, 0, 300
paste image 4, 512, 300
wait key
end
function ConvertToFloatMap(mem, tomem, offset, width, height)
`go off all pixels
for y = 0 to height - 1
for x = 0 to width - 1
pos = (y * width + x) * 4
write memblock float tomem, pos, memblock byte(mem, pos + offset)
next x
next y
endfunction
function ConvertToByteMap(mem, tomem, offset, width, height, factor#)
`Go off all floats
for y = 0 to height - 1
for x = 0 to width - 1
pos = (y * width + x) * 4
v = memblock float(mem, pos) * factor#
if v < 0 then v = 0
if v > 255 then v = 255
write memblock dword tomem, pos + offset, rgb(v, v, v)
next x
next y
endfunction
`width and height in complex samples
function ConvertComplexToAmplitude(mem, tomem, width, height, factor#)
for y = 0 to height - 1
for x = 0 to width - 1
pos = (y*width + x) * 8
a# = memblock float(mem, pos)^2 + memblock float(mem, pos + 4)^2
pos = (y*width + x) * 4
write memblock float tomem, pos, sqrt(a#) * factor#
next x
next y
endfunction
An example of a 2D FFT. By filtering the high frequency components, we can "blur" the image. However, due to the nature of the filter (a rectangular filter), we have artifacts near sharp edges of the image.
That was all. Have fun everyone!
Sven B
[edit]Forgot attachment.
[edit2] Updated with missing libfftw3f-3.dll file and updated instructions.