Dark Bit Factory & Gravity
PROGRAMMING => General coding questions => Topic started by: va!n on October 19, 2007
-
hi guys! me and a friend are trying to code something, where we need a fast FFT algo... we searched the web and found a lot of algos but sadly all are written cryted math signs (like on wikipedia), we dont understand anymore... (to long ago, when we had some of this math signs in scool ^^)
So, does someone of you know or have a working and easy to read and undestanding codesnip that works? atm we are trying to understand and convert a VB source with bad variablenames and without comments ^^
First attached image (scope) looks like something, we still managed to work atm...
Second attached image (fft) is what we are looking for (in realtime !)
many thanks in advance!!!
-
Morning dudes,
@va!n: I'd recommend having a look at this url: http://www.musicdsp.org/archive.php?classid=2#79
I'm pretty sure you'll find something useful there as this site is a real cool knowledgebase ... :)
SLiPPY
-
Wow fantastic link! Karma++ , thanks slippy.
-
Also http://www.nrbook.com/a/bookcpdf.php (http://www.nrbook.com/a/bookcpdf.php)
Chapter 12. Watch out, code works but coding style is shocking!
Jim
-
Jim -> numerical correct & fast algo is always horrible to read :P.
Btw if you needed basic FFT explanation try cormen's book.
-
Yeah, I know :) The problem is the code was ported from FORTRAN to C by the authors, and that makes it far worse that if it had been re-coded cleanly.
Jim
-
thx guys... i know the url dsp... but i dont found something that is easy to understand / useable for me... i will try it again ;)
-
here is some simple FFT processing (although it can be optimized further)
It has simple code to create different classic analyzers from the basic FFT
-
The Numerical Recipes code is actually pretty straightforward, but not very optimised, and the accompanying explanation is pretty good - if you're not grokking that then it might not matter how simple the code is. Can you post some attempts or ideas, or plans that you have for using the code?
Jim
-
@Stormbringer:
thanks for the code.... atm i dont really understand it (need to convert it to PB or a working (very fast and extreme small in size compiled lib ^^) however, own PB code would be much better, even to understand the math and how it works ^^
btw, is is possible, there are 2 differend kind of sound FFT algos available? One floating and one integer version? O.O If so, i think i only need a fast integer version ;) 10x in advance... best regards
-
@va!n: you are welcome. The code to perform the 1D FFT is quite simple actually. It takes a window of samples (BTW, can be any signal, not just audio) and does the transform forwards and backwards. The floating point version is slower but propably simpler to understand as the basic idea is to transform real numbers into complex numbers (therefore sin() and cosin() functions are used). I suggest you have a look at some theory here: http://zone.ni.com/devzone/cda/tut/p/id/4278
The explanation is quite good in my opinion.
Try the floating point version with some basic signal (a sine function for which you know the frequency) first. Let me know how it goes.
-
@stormbringer:
thanks a lot... i will try it out very soon... first i have to send my notebook to the manufractor to repair the keyboard and wait to get it hopefully very soon back :P
-
Due fact we need a very very fast and small size FFT algo, i found that intel Math Lib (comercial) has a very fast FFT inbuild... O.O Does someone has a very fast optimized version like intels or something like this? ^^ thx (source and compiled lib would be nice :P)
-
va!n ... I just found another solution which might work for you as well ... as Hitchhikr updated his site, he released the sourcecode of a music disk (Scoopex'n'Chips #3 Final) ... there's an fft algo built in ... perhaps it might help if you study these sources ...
http://pagesperso-orange.fr/franck.charlet/demoscene.html
Cheers,
SLiPPY
-
@slippy:
thanks, i will try to take a look to his source! K++
Edit:
Ah as i see, he is using bass.dll and i think there is a inbuild FFT like in fmod which is not acceptable to use any external DLLs like Bass, Fmod ^^ Its for a tiny exe size projec. anyway thx
-
You really should try the floating point code from Numerical Recipes. I used it in 1992 on a 20MHz Sun Sparc and was able to FFT 1MB or so in 20seconds. Today's CPUs are way faster, and you are only getting around 800samples/frame or 44100/second with audio. You don't need the Intel one, or anything near it - it's total overkill.
Jim
-
Are you sure you actually need an FFT ? I guess you want to detect the beats and that sort of things. Reading the patterns of your track, or synth, should be enough.
-
Va!n,
how about this: http://local.wasp.uwa.edu.au/~pbourke/other/dft/ ? Both dft and fft are in the source half way down the page...
Chris
-
Oh, nice link!
Jim
-
sorry but i am back again with this topic... first i want to thank you all for some nice links and your help... my main problem is to read/understand and convert this crypted-looked math symbols. (sadly my best friend and me cant read them like following graphic
(http://upload.wikimedia.org/math/9/e/a/9ea4fa6829333c2a4cba78e474bdb6a4.png))
However, until now we did not managed to code an own FFT nor to understand how does it work, because readable basic sources like VB-Sources are using silly variable-names where we dont know what they stand for and whats going on in the procedure.
I have to thanks JIM, who showed me a C/CPP source for DFT, i tried to convert to PB but we still have even problems with this DFT source. here is what i have (DFT instead needed (fast and small FFT) ^^
Type : fourier transform
References : Posted by Andy Mucho
Code :
AnalyseWaveform(float *waveform, int framesize)
{
float aa[MaxPartials];
float bb[MaxPartials];
for(int i=0;i<partials;i++)
{
aa[i]=0;
bb[i]=0;
}
int hfs=framesize/2;
float pd=pi/hfs;
for (i=0;i<framesize;i++)
{
float w=waveform[i];
int im = i-hfs;
for(int h=0;h<partials;h++)
{
float th=(pd*(h+1))*im;
aa[h]+=w*cos(th);
bb[h]+=w*sin(th);
}
}
for (int h=0;h<partials;h++)
amp[h]= sqrt(aa[h]*aa[h]+bb[h]*bb[h])/hfs;
}
Type : fourier transform
References : Posted by Andy Mucho
PB conversion: Mr.Vain / Secretly!
Code :
Procedure AnalyseWaveform( *waveform.f, framesize.l )
aa.f ( MaxPartials )
bb.f ( MaxPartials )
For i = 0 To partials
aa(i) = 0
bb(i) = 0
Next
hfs.l = framesize/2
pd.f = #Pi/hfs
For i = 0 to framzesize
w.f = waveform( i )
im.l = i-hfs
For h = 0 To partials
th.f = (pd*(h+1))*im
aa( h ) = aa( h ) + w*Cos( th )
bb( h ) = aa( h ) + w*Sin( th )
Next
Next
For h = 0 to partials
amp( h ) = Sqr( aa(h)*aa(h)+bb(h)*bb(h) ) /hfs
Next
EndProcedure
I hope i converted the source correctly to PB. Now lets have a look to the sources and try to explain following *confusing* things, we dont understand pls ^^
a)
There is the variable named MaxPartials, but what is it for and what value does it must have?
aa.f ( MaxPartials )
bb.f ( MaxPartials )
b)
Why does we must fill aa(i) and bb(i) with 0? If they will be initialized the first time, they will be always 0. Seems to me, its useless code that can be removed...?
For i = 0 To partials
aa(i) = 0
bb(i) = 0
Next
c)
What does the procedure should return, as it can only return one value each call...
Does someone have an own or know a very fast and very small FFT *.lib? Credits will be given! Thanks!
-
All those squiggles are trying to tell you that the frequency at a given point in time is the sum of all the sub-frequencies in your signal. Don't let it bother you too much...
In answer to your questions.
a) MaxPartials/partials is how much detail you want to analyse - the number of sub-frequencies you want. Set it to something like 128. For every incoming chunk of signal (say, 1/50th second's worth of audio), you get an ARRAY of MaxPartials as the output.
b) yes, you need to zero those arrays. if pb does that for you then you don't need to do it too.
c) the output is an array. it's stored in 'amp'. If you chose 128 partials, then effectively you are dividing the incoming signal into 128 frequencies. The values in 'amp' are the magnitudes of each of those sub-frequencies.
For i = 0 to framzesize
Typo?
If you plot the array vertically against time, (so, 128 dark to light pixels based on the values in 'amp' vertically per 1/50th of a second) you will get a pretty picture - something you might well recognise.
If it's too slow to start with, drop the partials to 64 or 32.
You don't need an FFT library, you have everything you need!
Jim
-
Since yesterday i am sitting by my best friend in duisburg and we are trying to code and finish some missing routines for our actually project like a good and fast FFT routine. From yesterday until this morning we have tried to get convert and working a FFT version, which seems to work very nice. *yeah*
So i will post the source(s) to help hopefully other peoples and possible some other guys here may be help to optimize the code in speed =)
Here is the original VB sourcecode we have tried to convert:
Dim REX(512) 'REX[ ] holds the real part of the frequency domain
Dim IMX(512) 'IMX[ ] holds the imaginary part of the frequency domain
Const N = 512
Public Sub fft()
pi = 3.14159265 'Set constants
1000 'THE FAST FOURIER TRANSFORM
'copyright © 1997-1999 by California Technical Publishing
'published with permission from Steven W Smith, www.dspguide.com
'GUI by logix4u , www.logix4u.net
'modified by logix4u, www.logix4.net
1010 'Upon entry, N% contains the number of points in the DFT, REX[ ] and
1020 'IMX[ ] contain the real and imaginary parts of the input. Upon return,
1030 'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1.
1060 NM1% = N% - 1
1070 ND2% = N% / 2
1080 M% = CInt(Log(N%) / Log(2))
1090 J% = ND2%
1100 '
1110 For i% = 1 To N% - 2 'Bit reversal sorting
1120 If i% >= J% Then GoTo 1190
1130 TR = REX(J%)
1140 TI = IMX(J%)
1150 REX(J%) = REX(i%)
1160 IMX(J%) = IMX(i%)
1170 REX(i%) = TR
1180 IMX(i%) = TI
1190 K% = ND2%
1200 If K% > J% Then GoTo 1240
1210 J% = J% - K%
1220 K% = K% / 2
1230 GoTo 1200
1240 J% = J% + K%
1250 Next i%
1260 '
1270 For L% = 1 To M% 'Loop for each stage
1280 LE% = CInt(2 ^ L%)
1290 LE2% = LE% / 2
1300 UR = 1
1310 UI = 0
1320 SR = Cos(pi / LE2%) 'Calculate sine & cosine values
1330 SI = -Sin(pi / LE2%)
1340 For J% = 1 To LE2% 'Loop for each sub DFT
1350 JM1% = J% - 1
1360 For i% = JM1% To NM1% Step LE% 'Loop for each butterfly
1370 IP% = i% + LE2%
1380 TR = REX(IP%) * UR - IMX(IP%) * UI 'Butterfly calculation
1390 TI = REX(IP%) * UI + IMX(IP%) * UR
1400 REX(IP%) = REX(i%) - TR
1410 IMX(IP%) = IMX(i%) - TI
1420 REX(i%) = REX(i%) + TR
1430 IMX(i%) = IMX(i%) + TI
1440 Next i%
1450 TR = UR
1460 UR = TR * SR - UI * SI
1470 UI = TR * SI + UI * SR
1480 Next J%
1490 Next L%
1500 '
End Sub
Here is the converted ( /2 instead >> 1 unoptimized) PB version, which seems to work fine:
pi.f = 3.14159265
N = 8 ; Number of samples
Dim rex.f( N - 1 )
Dim imx.f( N - 1 )
; Upon entry, N contains the number of points in the DFT, REX() and
; IMX() contain the real and imaginary parts of the input. Upon return,
; REX() and IMX() contain the DFT output. All signals run from 0 to N-1.
NM1.l = N - 1
ND2.l = N / 2
M.l = Int( Log( N ) / Log( 2 ) )
J.l = ND2
For i.l = 1 To N - 2 ; Bit reversal sorting
If i < J
TR.f = REX( J )
TI.f = IMX( J )
REX( J ) = REX( i )
IMX( J ) = IMX( i )
REX( i ) = TR
IMX( i ) = TI
EndIf
K = ND2
While K <= J
J = J - K
K = K / 2
Wend
J = J + K
Next
For L = 1 To M ; Loop for each stage
LE.l = Int( Pow( 2, L ) )
LE2.l = LE / 2
UR.f = 1
UI.f = 0
SR.f = Cos( #PI / LE2 ) ; Calculate sine & cosine values
SI.f = - Sin( #PI / LE2 )
For J.l = 1 To LE2 ; Loop for each sub DFT
JM1.l = J - 1
For i = JM1 To NM1 ; Loop for each butterfly
IP.l = i + LE2
TR = REX( IP ) * UR - IMX( IP ) * UI ; Butterfly calculation
TI = REX( IP ) * UI + IMX( IP ) * UR
REX( IP ) = REX( i ) - TR
IMX( IP ) = IMX( i ) - TI
REX( i ) = REX( i ) + TR
IMX( i ) = IMX( i ) + TI
i + LE - 1
Next
TR = UR
UR = TR * SR - UI * SI
UI = TR * SI + UI * SR
Next
Next
And at least here is my latest try to convert a small modified version from C to PB:
; / ***************************************************************************/
; / * fft.pb */
; / ***************************************************************************/
; / * Forward And Reverse Fast Fourier Transform routines. */
; / * ==================================================== */
; / * Original Source by Tobias Dammers (c) 2005. */
; / * PB version by Thorsten Will aka va!n in 2008 */
; / * */
; / * Reference: http://wiki.allegro.cc/FFT */
; / * */
; / * Feel free To use these routines As you see fit. */
; / ***************************************************************************/
; // Fourier Transform is the process of transforming signals from the time domain
; // (levels over time) into the frequency domain (levels over frequency). The dis-
; // crete-time version of this is called the Discrete Fourier Transform (DFT).
; // The standard approach is a rather costly operation, but For blocks of Data
; // whose length is a power of 2, a dramatic optimization is possible, resulting
; // in an algorithm called Fast Fourier Transform, FFT. This is a C implementation
; // of the algorithm, And of its reverse process.
; // The forward FFT. rex holds the real part of the signal, imx the imaginary part.
; // Both arrays must be at least 2^m samples long. The result is returned in the
; // same two arrays, ie the transform is in-place. After transforming, rex holds
; // the real part of the frequency domain, And imx the imaginary part.
; // To calculate the frequency spectrum of an arbitrary block of wave Data, put the
; // input in rex[], And set imx[] To all-zero.
Global Dim rex.f( 512 )
Global Dim imx.f( 512 )
Procedure FFT( m.l )
Protected nm1.l, nd2.l
Protected i.l, j.l, l.l, n.l
Protected tr.f, ti.f
n.l = 1 << m
nm1.l = n - 1
nd2.l = n >> 1
j.l = nd2
tr.f
ti.f
i.l
l.l
i = 1
While i < n
k.l
If ( i < j )
tr = rex( j )
ti = imx( j )
rex( j ) = rex( i )
imx( j ) = imx( i )
rex( i ) = tr
imx( i ) = ti
EndIf
k = nd2
While ( ( k <= j ) And ( k > 0 ) ) ; &&
j = j - k
k = k >> 1
Wend
j = j + k
i +1
Wend
l = 1
While l <= m
le.l = 1 << l
le2.l = le >> 1
ur.f = 1.0
ui.f = 0.0
sr.f = Cos( #PI / le2 )
si.f = -Sin( #PI / le2 )
j = 1
While j <= le2
i = j-1
While i <= nm1
ip.l = i + le2
tr = rex( ip ) * ur - imx( ip ) * ui
ti = rex( ip ) * ui + imx( ip ) * ur
rex( ip ) = rex( i ) - tr
imx( ip ) = imx( i ) - ti
rex( i ) = rex( i ) + tr
imx( i ) = imx( i ) + ti
i = i + le
Wend
tr = ur
ur = tr * sr - ui * si
ui = tr * si + ui * sr
j + 1
Wend
l + 1
Wend
EndProcedure
; // The inverse FFT. This function takes advantage of the special symmetry of
; // Fourier pairs; all we need to do is scale the input, apply a forward FFT
; // then scale the output.
; // As a test, you can feed a signal through the forward fft, then feed the result
; // back through the inverse fft, giving you the original signal, give Or take
; // some roundoff errors.
Procedure Inv_FFT( m.l )
n.l = 1 << m
i.l
i = 0
While i < n
imx( i ) = -imx( i )
i + 1
Wend
FFT( m )
i = 0
While i < n
rex( i ) = rex( i ) / n ; rex[i] /= (float)n;
imx( i ) = imx( i ) / -n ; imx[i] /= -(float)n
i + 1
Wend
EndProcedure
If some of you have a benchmark for this both routines (which is the fastest) or how and where to optimize the code (basic/ x86 or MMX or SSEx) it would be really cool. Thanks for any feedback. Hope this source(s) may be helpfull for some guys here.
-
Thats a karma++ moment.
-
Va!n - there are lots of places in your code where you divide by something. Often that division could be replaced by a multiply by the reciprocal (ie. x = a/n, y=b/n could be re-written t = 1/n, x=a*t, y=b*t). And often, that reciprocal calculation can be hauled outside the loops.
Also, you should replace the separate calls to Sin and Cos library functions to use the inline assembly instruction fsincos which calculates both values at the same time.
For instance, I have this code
static void sincos(double a, double *s, double *c)
{
double zin,coz;
//a *= 3.141592653589/180.0;
_asm
{
fld [a]
fsincos
fstp [coz]
fstp [zin]
}
*s = zin;
*c = coz;
}
in most of my demos. You need the multiply if you're working in degrees (which you're not).
Jim
-
@taj:
Thanks for k++ and feedback =)
@JIM:
Thanks for feedback and the great SIN/COS inline asm hint, i will try to add this later. *thumbs up* K++ for the tipp! =)
@all:
I managed it to optimize the FFT code a little bit more. I will try to optimize the SIN/COS stuff with inline asm later.
; // *************************************************************************************
; // *
; // * Project: FFT (THE FAST FOURIER TRANSFORM) v0.1b
; // *
; // * Original version: logix4u, www.logix4.net
; // * Coverted version: Tranquil and Mr.Vain of Secretly! (Thorsten Will) in 2008
; // * Optimisations by: va!n aka Thorsten Will
; // *
; // * Reference:
; // * http://logix4u.net/DSP/Fast_Fourier_Transform/Visual_Basic_program_for_Fast_Fourier_Transform.html
; // *
; // * Upon entry, N.w contains the number of points in the DFT, REX() and
; // * IMX() contain the real and imaginary parts of the input. Upon return,
; // * REX() and IMX() contain the DFT output. All signals run from 0 to N.w-1
; // *
; // *************************************************************************************
Procedure record_doFFT( *scope.SCOPE )
; // -------- Init some values for FFT analysing --------
N.w = 1024 ; // Number of samples
NM1.l = N - 1
ND2.l = N >> 1 ; // Optmimized, instead N / 2
M.l = 10 ; // If N.w = 1024 == 10 == Same as: Int(Log(N) / 0.69314718055994529)
J.l = ND2
If Config\buffer = 0 : ProcedureReturn : EndIf
; // -------- Clear and Fill array values for analysing in just only one loop --------
For i = 0 To Config\size Step 2 ; // Optimized by merging clear and fill array in one loop
rex( i >> 1 ) = 0 ; // 0 to 512
imx( i >> 1 ) = 0 ; // 0 to 512
rex( i >> 1 + N >> 1 ) = 0 ; // 513 to 1024
imx( i >> 1 + N >> 1 ) = 0 ; // 513 to 1024
value.w = PeekW( Config\buffer + i )
; value.w = PeekW( Config\buffer + i + *scope\channel * 2 ) ; // Enable this For Stereo Inpus
rex( i >> 1 ) = value / 32767 ; // Optimized by doing i >> 1
; imx( i >> 1 ) = 0 ; // Can be removed now
Next
; // -------- Start FFT --------
For i.l = 1 To N - 2 ; // Bit reversal sorting
If i < J
TR.f = REX( J )
TI.f = IMX( J )
REX( J ) = REX( i )
IMX( J ) = IMX( i )
REX( i ) = TR
IMX( i ) = TI
EndIf
K = ND2
While K <= J
J = J - K
K = K >> 1 ; // Optmimized, instead N / 2
Wend
J = J + K
Next
For L = 1 To M ; // Loop for each stage
LE.l = Int( Pow( 2, L ) )
LE2.l = LE >> 1 ; // Optmimized, instead N / 2
UR.f = 1
UI.f = 0
SR.f = Cos( #PI / LE2 ) ; // Calculate sine & cosine values
SI.f = - Sin( #PI / LE2 )
For J.l = 1 To LE2 ; // Loop for each sub DFT
JM1.l = J - 1
For i = JM1 To NM1 ; // Loop for each butterfly
IP.l = i + LE2
TR = REX( IP ) * UR - IMX( IP ) * UI ; // Butterfly calculation
TI = REX( IP ) * UI + IMX( IP ) * UR
REX( IP ) = REX( i ) - TR
IMX( IP ) = IMX( i ) - TI
REX( i ) = REX( i ) + TR
IMX( i ) = IMX( i ) + TI
i + LE - 1
Next i
TR = UR
UR = TR * SR - UI * SI
UI = TR * SI + UI * SR
Next
Next
; // -------- Calculate Outputarray --------
For cnt = 0 To N.w ; fixed to N.w instead wrong fixed value
outputarray( cnt ) = ( IMX( cnt ) * IMX( cnt ) ) + ( REX( cnt ) * REX( cnt ) )
Next
; // -------- Search for MaxValue of the Paket --------
maxvalue = 0
For cnt = 0 To N.w ; // fixed to N.w instead wrong fixed value
If maxvalue < outputarray( cnt )
maxvalue = outputarray( cnt )
EndIf
Next
; // -------- Draw FFT --------
StartDrawing( WindowOutput( FFTWnd ) )
Box( 0, 0, N.w, 500, $0 )
MaxPeak = 0
For cnt = 0 To 500 ; ; // Change fixed value to N.w !?
DiffY = Outputarray( cnt ) / MaxValue * 400
Box( cnt, 400, 1, -DiffY, $FFFFFF)
If DiffY > MaxPeak
MaxPeak = DiffY
MaxPos = cnt
EndIf
Next
StopDrawing()
SetWindowTitle( FFTWnd, ShowNote_Get( record_FindNote( MaxPos ) ) )
EndProcedure
-
Okay guys... i am sure some of you are interested in FFT coding, so here we did a small, fast test application (sound recorder) to test our FFT routine and its results. Here the full PB sourcecode. An executeable version will be added for testing too.
Sourcecode using our first "FFT v0.1" version:
Global Dim rex.f(512*2+1)
Global Dim imx.f(512*2+1)
Global Dim OutPutArray.f(512*2+1)
Global FFTWnd
Structure NoteRange
Note.l
FromPos.l
ToPos.l
EndStructure
Procedure.l ShowNote_Init()
Global Dim NoteRange.NoteRange(53)
For Note=0 To 53
Read FromPos.w
Read ToPos.w
NoteRange(Note)\FromPos = FromPos
NoteRange(Note)\ToPos = ToPos
NoteRange(Note)\Note = Note
Next
Global Dim g_RealNote.s( 53 )
For i = 0 To 42+12-1
Read sRealNote.s
g_RealNote( i ) = sRealNote.s
Next
EndProcedure
Procedure.s ShowNote_Get( lValue)
ProcedureReturn g_RealNote.s( lValue )
EndProcedure
ShowNote_Init()
Structure SCOPE
channel.b
left.l
top.l
width.l
height.l
middleY.l
quarterY.l
EndStructure
Structure CONFIG
hWindow.l ; Window handle
size.l ; Wave buffer size
buffer.l ; Wave buffer pointer
output.l ; WindowOutput()
wave.l ; Address of waveform-audio input device
format.WAVEFORMATEX ; Capturing WaveFormatEx
lBuf.l ; Capturing Buffer size
nBuf.l ; Capturing Buffer number
nDev.l ; Capturing Device identifier
nBit.l ; Capturing Resolution (8/16)
nHertz.l ; Capturing Frequency (Hertz)
nChannel.l ; Capturing Channels number (Mono/Stereo)
LScope.SCOPE ; Wave form display
RScope.SCOPE ; Wave form display
EndStructure
Global Config.CONFIG
Global Dim inHdr.WAVEHDR(16)
Config\format\wFormatTag = #WAVE_FORMAT_PCM
Procedure Record_Start()
Config\format\nChannels = 1
Config\format\wBitsPerSample = 16
Config\format\nSamplesPerSec = 8000
Config\nDev = 0 ; (0 default MS Sound Mapper device)
Config\lBuf = 1024
Config\nBuf = 8
Config\nBit = 1
Config\format\nBlockAlign = (Config\format\nChannels*Config\format\wBitsPerSample)/8
Config\format\nAvgBytesPerSec = Config\format\nSamplesPerSec*Config\format\nBlockAlign
If #MMSYSERR_NOERROR = waveInOpen_(@Config\wave,#WAVE_MAPPER+Config\nDev,@Config\format,Config\hWindow,#Null,#CALLBACK_WINDOW|#WAVE_FORMAT_DIRECT)
For i=0 To Config\nBuf-1
inHdr(i)\lpData=AllocateMemory(Config\lBuf)
inHdr(i)\dwBufferLength=Config\lBuf
waveInPrepareHeader_(Config\wave,inHdr(i),SizeOf(WAVEHDR))
waveInAddBuffer_(Config\wave,inHdr(i),SizeOf(WAVEHDR))
Next
If #MMSYSERR_NOERROR = waveInStart_(Config\wave)
SetTimer_(Config\hWindow,0,1,0)
EndIf
EndIf
EndProcedure
Procedure Record_Read(hWaveIn.l,lpWaveHdr.l)
*hWave.WAVEHDR=lpWaveHdr
Config\buffer=*hWave\lpData
Config\size=*hWave\dwBytesRecorded
waveInAddBuffer_(hWaveIn,lpWaveHdr,SizeOf(WAVEHDR))
EndProcedure
Procedure record_FindNote(Value)
For Note = 0 To 53
If Value=>NoteRange(Note)\FromPos And Value<=NoteRange(Note)\ToPos
ProcedureReturn note
EndIf
Next
EndProcedure
Procedure record_doFFT(*scope.SCOPE)
If Config\buffer = 0 : ProcedureReturn : EndIf
For pos=0 To 1024:rex(pos)=0:imx(pos)=0:Next
pos = 0
For i=0 To Config\size Step 2
value.w=PeekW(Config\buffer+i)
;value.w=PeekW(Config\buffer+i+*scope\channel*2) Enable this for Stereo Inpus
rex(pos) = value/32767
imx(pos) = 0 :
pos + 1
Next
N.w = 1024 ; Num Samples
;N.w = 512
m-w = 102
NM1.l = N - 1
ND2.l = N / 2
M.l = Int(Log(N) / 0.69314718055994529)
J.l = ND2
For i.l = 1 To N - 2 ; Bit reversal sorting
If i < J
TR.f = REX(J)
TI.f = IMX(J)
REX(J) = REX(i)
IMX(J) = IMX(i)
REX(i) = TR
IMX(i) = TI
EndIf
K = ND2
While K <= J
J = J - K
K = K / 2
Wend
J = J + K
Next i
For L = 1 To M ; Loop for each stage
LE.l = Int(Pow(2, L))
;LE2.l = LE / 2
LE2.l = LE >> 1
UR.f = 1
UI.f = 0
SR.f = Cos(#PI / LE2) ; Calculate sine & cosine values
SI.f = - Sin(#PI / LE2)
For J.l = 1 To LE2 ; Loop for each sub DFT
JM1.l = J - 1
For i = JM1 To NM1 ; Loop for each butterfly
IP.l = i + LE2
TR = REX(IP) * UR - IMX(IP) * UI ; Butterfly calculation
TI = REX(IP) * UI + IMX(IP) * UR
REX(IP) = REX(i) - TR
IMX(IP) = IMX(i) - TI
REX(i) = REX(i) + TR
IMX(i) = IMX(i) + TI
i + LE - 1
Next i
TR = UR
UR = TR * SR - UI * SI
UI = TR * SI + UI * SR
Next J
Next L
; Outputarray berechnen
For cnt=0 To 512
outputarray(cnt) = (IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt))
Next cnt
;Search for MaxValue of the Paket
maxvalue = 1
;
For cnt = 0 To 1024
If (maxvalue < outputarray(cnt))
maxvalue = outputarray(cnt)
EndIf
Next cnt
StartDrawing(WindowOutput(FFTWnd))
Box(0,0,500,500,$0)
MaxPeak=0
For cnt = 5 To 512
DiffY=Outputarray(cnt)/MaxValue*400
LineXY(cnt,400,cnt,400-DiffY,$FFFFFF)
If DiffY>MaxPeak
MaxPeak=DiffY
AkNote=cnt
EndIf
Next cnt
SetWindowTitle(FFTWnd,Str(record_FindNote(aknote))+" on:"+Str(aknote)+" note:"+ShowNote_Get(record_FindNote(aknote)))
StopDrawing()
EndProcedure
Procedure record_CallBack(hWnd.l,Msg.l,wParam.l,lParam.l)
Result.l=#PB_ProcessPureBasicEvents
Select Msg
Case #WM_TIMER : record_doFFT(Config\LScope)
Case #MM_WIM_DATA : record_Read(wParam,lParam)
EndSelect
ProcedureReturn Result
EndProcedure
FFTWnd = OpenWindow(#PB_Any,0,0,500,500,"",#PB_Window_ScreenCentered | #PB_Window_SystemMenu)
Config\hWindow=WindowID(FFTWnd)
Config\output=WindowOutput(FFTWnd)
SetWindowCallback(@record_CallBack())
Record_Start()
Repeat:Until WaitWindowEvent() = #PB_Event_CloseWindow
End
DataSection
Notes:
Data.w 14, 17
Data.w 18, 18
Data.w 19, 19
Data.w 20, 20
Data.w 21, 21
Data.w 22, 23
Data.w 24, 24
Data.w 25, 26
Data.w 27, 27
Data.w 28, 29
Data.w 30, 31
Data.w 32, 32
Data.w 33, 34
Data.w 35, 37
Data.w 38, 39
Data.w 40, 41
Data.w 42, 44
Data.w 45, 46
Data.w 47, 49
Data.w 50, 52
Data.w 53, 55
Data.w 56, 59
Data.w 60, 62
Data.w 63, 66
Data.w 67, 70
Data.w 71, 74
Data.w 75, 79
Data.w 80, 83
Data.w 84, 88
Data.w 89, 94
Data.w 95, 99
Data.w 100, 105
Data.w 106, 112
Data.w 113, 118
Data.w 119, 125
Data.w 126, 133
Data.w 134, 141
Data.w 142, 149
Data.w 150, 158
Data.w 159, 168
Data.w 169, 178
Data.w 179, 188
Data.w 189, 200
Data.w 201, 212
Data.w 213, 224
Data.w 225, 238
Data.w 239, 252
Data.w 253, 267
Data.w 268, 283
Data.w 284, 300
Data.w 301, 318
Data.w 319, 337
Data.w 338, 357
Data.w 358, 375
RealNotes:
Data.s "C0", "C#0", "D0", "D#0", "E0", "F0", "F#0", "G0", "G#0"
Data.s "A0", "A#0", "B0", "C1", "C#1", "D1", "D#1", "E1", "F1", "F#1", "G1", "G#1"
Data.s "A2", "A#2", "B2", "C2", "C#2", "D2", "D#2", "E2", "F2", "F#2", "G2", "G#2"
Data.s "A3", "A#3", "B3", "C3", "C#3", "D3", "D#3", "E3", "F3", "F#3", "G3", "G#3"
Data.s "A4", "A#4", "B4", "C4", "C#4", "D4", "D#4", "E4", "F4"
EndDataSection
"FFT v0.1c" is the latest speed optimized version procedure i got managed to work. ASM optimisations are not impented yet.
; // *************************************************************************************
; // *
; // * Project: FFT (THE FAST FOURIER TRANSFORM) v0.1c
; // *
; // * Original version: logix4u, www.logix4.net
; // * Coverted version: Tranquil and Mr.Vain (Thorsten Will) in 2008
; // * Optimisations by: va!n aka Thorsten Will
; // *
; // * Special thanks 2: MrMat, Jim, Hellhound
; // *
; // * Reference:
; // * http://logix4u.net/DSP/Fast_Fourier_Transform/Visual_Basic_program_for_Fast_Fourier_Transform.html
; // *
; // * Upon entry, N.w contains the number of points in the DFT, REX() and
; // * IMX() contain the real and imaginary parts of the input. Upon return,
; // * REX() and IMX() contain the DFT output. All signals run from 0 to N.w-1
; // *
; // *************************************************************************************
Procedure record_doFFT( *scope.SCOPE )
Define.d TR, TI, SR, SI, UR, UI
Define.l J, K, L, M, N, NM1, ND2, cnt
Define.L MaxPeak, MaxValue, DiffY
Define.w value
; // -------- Init some values for FFT analysing --------
N = 1024 ; // Number of samples
NM1 = N - 1
ND2 = N >> 1 ; // Optmimized, instead N / 2
M = 10 ; // If N.w = 1024 == 10 == Same as: Int(Log(N) / 0.69314718055994529)
J = ND2
If Config\buffer = 0 : ProcedureReturn : EndIf
; // -------- Clear and Fill array values for analysing in just only one loop --------
For i = 0 To Config\size Step 2 ; // Optimized by merging clear and fill array in one loop
rex( i >> 1 ) = 0 ; // 0 to 512
imx( i >> 1 ) = 0 ; // 0 to 512
rex( i >> 1 + N >> 1 ) = 0 ; // 513 to 1024
imx( i >> 1 + N >> 1 ) = 0 ; // 513 to 1024
value = PeekW( Config\buffer + i )
; value.w = PeekW( Config\buffer + i + *scope\channel * 2 ) ; // Enable this For Stereo Inpus
rex( i >> 1 ) = value / 32767 ; // Optimized by doing i >> 1
; imx( i >> 1 ) = 0 ; // Can be removed now
Next
; // -------- Start FFT --------
For i.l = 1 To N - 2 ; // Bit reversal sorting
If i < J
TR = REX( J )
TI = IMX( J )
REX( J ) = REX( i )
IMX( J ) = IMX( i )
REX( i ) = TR
IMX( i ) = TI
EndIf
K = ND2
While K <= J
J = J - K
K = K >> 1 ; // Optmimized, instead N / 2
Wend
J = J + K
Next
For L = 1 To M ; // Loop for each stage
LE.l = 1 << L ; // Optimized, instead LE.l = Int( Pow( 2, L ) )
LE2.l = LE >> 1 ; // Optimized, instead N / 2
UR = 1
UI = 0
SR = Cos( #PI / LE2 ) ; // Calculate sine & cosine values
SI = - Sin( #PI / LE2 )
For J.l = 1 To LE2 ; // Loop for each sub DFT
JM1.l = J - 1
For i = JM1 To NM1 ; // Loop for each butterfly
IP.l = i + LE2
TR = REX( IP ) * UR - IMX( IP ) * UI ; // Butterfly calculation
TI = REX( IP ) * UI + IMX( IP ) * UR
REX( IP ) = REX( i ) - TR
IMX( IP ) = IMX( i ) - TI
REX( i ) = REX( i ) + TR
IMX( i ) = IMX( i ) + TI
i + LE - 1
Next i
TR = UR
UR = TR * SR - UI * SI
UI = TR * SI + UI * SR
Next
Next
; // -------- Calculate Outputarray and search for MaxValue of the Paket --------
maxvalue = 0 ; // Optimized by merging calculate Outputarray
; // and search MaxValue into just one loop.
For cnt = 0 To N ; // fixed to N.w instead wrong fixed value
outputarray( cnt ) = ( IMX( cnt ) * IMX( cnt ) ) + ( REX( cnt ) * REX( cnt ) )
If maxvalue < outputarray( cnt )
maxvalue = outputarray( cnt )
EndIf
Next
; // -------- Draw FFT --------
StartDrawing( WindowOutput( FFTWnd ) )
Box( 0, 0, N, 500, $0 )
MaxPeak = 0
For cnt = 0 To 500 ; ; // Change fixed value to N.w !?
DiffY = Outputarray( cnt ) / MaxValue * 400
Box( cnt, 400, 1, -DiffY, $FFFFFF)
If DiffY > MaxPeak
MaxPeak = DiffY
MaxPos = cnt
EndIf
Next
StopDrawing()
SetWindowTitle( FFTWnd, ShowNote_Get( record_FindNote( MaxPos ) ) )
EndProcedure
ZIP attachement contains "FFT v0.1" (topic source == first version) and my latest "FFT v0.1c" versions with some speed optimisations (still not ASM optimized, yet)
-
great! :clap:
it's what I want to do with a module tracker, but in C. is it simple to adapt it in C for other player using wavout ?
-
@[Y]
Yes it should be easy to port to any other language... You just only have to put your output-sound-datas to the FFT input array to have the graph. good luck ;)
-
rex( i >> 1 + N >> 1 ) = 0 ; // 513 to 1024
imx( i >> 1 + N >> 1 ) = 0 ; // 513 to 1024
[code]
In C, this would be equivalent to
(i>>(1+N))>>1
where I think you want
(i>>1)+(N>>1)
Might be worth putting the brackets in :)
Your routine will now be spending nearly all its time in Sin and Cos library routines...time to get the asm in ;)
Jim
[/code]