Hello
I have found this 2 FFT algorithm.
The first seem to work good.
The second seem to NOT work good but should be more fast than the first.
There is some brave man that could verify if they work good and if really
the second is more fast?
Thanks
'-----------------------------------------------------------------------
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%) / 0.693147180559945)
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%
1270 For L% = 1 To M% 'Loop for each stage
1280 LE% = 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
420 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 '
'-------------------------------------------
'For z = A To B
' Debug.Print "REX ["; z; "] = "; REX(z)
' Debug.Print "IMX ["; z; "] = "; IMX(z)
' Debug.Print "---------------------------"
'Next
'--------------------------------------------
End Sub
'--------------------------------------------------------------------------
'--------------------------------------------------------------------------
Private Sub Real_FFT()
NH% = N% / 2 - 1
For I% = 0 To NH%
REX(I%) = REX(2 * I%)
IMX(I%) = REX(2 * I% + 1)
Next I%
N% = N% /2
Call FFT ' The previous subroutine
N% = N% *2
NM1% = N% - 1
ND2% = N% / 2
N4% = N% / 4 - 1
For I% = 1 To N4%
IM% = ND2% - I%
IP2% = I% + ND2%
IPM% = IM% + ND2%
REX(IP2%) = (IMX(I%) + IMX(IM%)) / 2
REX(IPM%) = REX(IP2%)
IMX(IP2%) = -(REX(I%) - REX(IM%)) / 2
IMX(IPM%) = -IMX(IP2%)
REX(I%) = (REX(I%) + REX(IM%)) / 2
REX(IM%) = REX(I%)
IMX(I%) = (IMX(I%) - IMX(IM%)) / 2
IMX(IM%) = -IMX(I%)
Next I%
REX(N% * 3 / 4) = IMX(N% / 4)
REX(ND2%) = IMX(0)
IMX(N% * 3 / 4) = 0
IMX(ND2%) = 0
IMX(N% / 4) = 0
IMX(0) = 0
PI = 3.14159265
L% = CInt(Log(N%) / Log(2))
LE% = CInt(2 ^ L%)
le2% = LE% / 2
UR = 1
UI = 0
SR = Cos(PI / le2%)
SI = -Sin(PI / le2%)
For J% = 1 To le2%
jm1% = J% - 1
For I% = jm1% To NM1% Step LE%
IP% = 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
Next I
TR = UR
UR = TR * SR - UI * SI
UI = TR * SI + UI * SR
Next J
'-------------------------------------------
'For z = A To B
' Debug.Print "REX ["; z; "] = "; REX(z)
' Debug.Print "IMX ["; z; "] = "; IMX(z)
' Debug.Print "---------------------------"
'Next
'--------------------------------------------
End Sub
'--------------------------------------------------------------------------
FFT algoritm who is best
Started by ●November 22, 2010
Reply by ●November 22, 20102010-11-22
On Nov 22, 4:35=A0pm, "gpezzella" <gpezzella@n_o_s_p_a_m.yahoo.com> wrote:> Hello > > I have found this 2 FFT algorithm. > > The first seem to work good. > The second seem to NOT work good but should be more fast than the first. > > There is some brave man that could verify if they work good and if really > the second is more fast? > > Thanks =A0 > > '----------------------------------------------------------------------- > Public Sub fft() > > PI =3D 3.14159265 'Set constants > > 1000 'THE FAST FOURIER TRANSFORM > 'copyright =A9 1997-1999 by California Technical Publishing > 'published with =A0permission 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% =3D N% - 1 > 1070 ND2% =3D N% / 2 > 1080 M% =3D CInt(Log(N%) / 0.693147180559945) > > 1090 J% =3D ND2% > 1100 ' > > 1110 For I% =3D 1 To N% - 2 'Bit reversal sorting > 1120 If I% >=3D J% Then GoTo 1190 > 1130 TR =3D REX(J%) > 1140 TI =3D IMX(J%) > 1150 REX(J%) =3D REX(I%) > 1160 IMX(J%) =3D IMX(I%) > 1170 REX(I%) =3D TR > 1180 IMX(I%) =3D TI > 1190 K% =3D ND2% > 1200 If K% > J% Then GoTo 1240 > 1210 J% =3D J% - K% > 1220 K% =3D K% / 2 > 1230 GoTo 1200 > 1240 J% =3D J% + K% > 1250 Next I% > > 1270 For L% =3D 1 To M% 'Loop for each stage > 1280 LE% =3D 2 ^ L% > 1290 le2% =3D LE% / 2 > 1300 UR =3D 1 > 1310 UI =3D 0 > 1320 SR =3D Cos(PI / le2%) 'Calculate sine & cosine values > 1330 SI =3D -Sin(PI / le2%) > > 1340 For J% =3D 1 To le2% 'Loop for each sub DFT > 1350 jm1% =3D J% - 1 > > 1360 For I% =3D jm1% To NM1% Step LE% 'Loop for each butterfly > 1370 IP% =3D I% + le2% > > 1380 TR =3D REX(IP%) * UR - IMX(IP%) * UI 'Butterfly calculation > 1390 TI =3D REX(IP%) * UI + IMX(IP%) * UR > > 1400 REX(IP%) =3D REX(I%) - TR > 1410 IMX(IP%) =3D IMX(I%) - TI > > 420 REX(I%) =3D REX(I%) + TR > 1430 IMX(I%) =3D IMX(I%) + TI > > 1440 Next I% > > 1450 TR =3D UR > 1460 UR =3D TR * SR - UI * SI > 1470 UI =3D TR * SI + UI * SR > > 1480 Next J% > > 1490 Next L% > 1500 ' > > '------------------------------------------- > 'For z =3D A To B > ' =A0 =A0Debug.Print "REX ["; z; "] =3D "; REX(z) > ' =A0 =A0Debug.Print "IMX ["; z; "] =3D "; IMX(z) > ' =A0 =A0Debug.Print "---------------------------" > 'Next > '-------------------------------------------- > > End Sub > '------------------------------------------------------------------------=--> > '------------------------------------------------------------------------=--> Private Sub Real_FFT() > > NH% =3D N% / 2 - 1 > > For I% =3D 0 To NH% > =A0 =A0REX(I%) =3D REX(2 * I%) > =A0 =A0IMX(I%) =3D REX(2 * I% + 1) > Next I% > > N% =3D N% /2 > > Call FFT ' The previous subroutine > > N% =3D N% *2 > > NM1% =3D N% - 1 > ND2% =3D N% / 2 > N4% =3D N% / 4 - 1 > For I% =3D 1 To N4% > =A0 =A0 IM% =3D ND2% - I% > =A0 =A0 IP2% =3D I% + ND2% > =A0 =A0 IPM% =3D IM% + ND2% > =A0 =A0 REX(IP2%) =3D (IMX(I%) + IMX(IM%)) / 2 > =A0 =A0 REX(IPM%) =3D REX(IP2%) > =A0 =A0 IMX(IP2%) =3D -(REX(I%) - REX(IM%)) / 2 > =A0 =A0 IMX(IPM%) =3D -IMX(IP2%) > =A0 =A0 REX(I%) =3D (REX(I%) + REX(IM%)) / 2 > =A0 =A0 REX(IM%) =3D REX(I%) > =A0 =A0 IMX(I%) =3D (IMX(I%) - IMX(IM%)) / 2 > =A0 =A0 IMX(IM%) =3D -IMX(I%) > Next I% > > REX(N% * 3 / 4) =3D IMX(N% / 4) > REX(ND2%) =3D IMX(0) > IMX(N% * 3 / 4) =3D 0 > IMX(ND2%) =3D 0 > IMX(N% / 4) =3D 0 > IMX(0) =3D 0 > > PI =3D 3.14159265 > L% =3D CInt(Log(N%) / Log(2)) > LE% =3D CInt(2 ^ L%) > > le2% =3D LE% / 2 > UR =3D 1 > UI =3D 0 > SR =3D Cos(PI / le2%) > SI =3D -Sin(PI / le2%) > For J% =3D 1 To le2% > jm1% =3D J% - 1 > > For I% =3D jm1% To NM1% Step LE% > =A0 IP% =3D I% * le2% > =A0 TR =3D REX(IP%) * UR - IMX(IP%) * UI > =A0 TI =3D REX(IP%) * UI + IMX(IP%) * UR > =A0 REX(IP%) =3D REX(I%) - TR > =A0 IMX(IP%) =3D IMX(I%) - TI > =A0 REX(I%) =3D REX(I%) + TR > =A0 IMX(I%) =3D IMX(I%) * TI > Next I > > TR =3D UR > UR =3D TR * SR - UI * SI > UI =3D TR * SI + UI * SR > Next J > > '------------------------------------------- > 'For z =3D A To B > ' =A0 =A0Debug.Print "REX ["; z; "] =3D "; REX(z) > ' =A0 =A0Debug.Print "IMX ["; z; "] =3D "; IMX(z) > ' =A0 =A0Debug.Print "---------------------------" > 'Next > '-------------------------------------------- > > End Sub > '------------------------------------------------------------------------=-- The first one is better because it has line numbers. The second is written in Klingon and is likely to be very badly behaved.
Reply by ●November 22, 20102010-11-22
On Nov 22, 11:35=A0am, "gpezzella" <gpezzella@n_o_s_p_a_m.yahoo.com> wrote:> Hello > > I have found this 2 FFT algorithm. > > The first seem to work good. > The second seem to NOT work good but should be more fast than the first. > > There is some brave man that could verify if they work good and if really > the second is more fast? > > Thanks =A0 > > '----------------------------------------------------------------------- > Public Sub fft() > > PI =3D 3.14159265 'Set constants > > 1000 'THE FAST FOURIER TRANSFORM > 'copyright =A9 1997-1999 by California Technical Publishing > 'published with =A0permission 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% =3D N% - 1 > 1070 ND2% =3D N% / 2 > 1080 M% =3D CInt(Log(N%) / 0.693147180559945) > > 1090 J% =3D ND2% > 1100 ' > > 1110 For I% =3D 1 To N% - 2 'Bit reversal sorting > 1120 If I% >=3D J% Then GoTo 1190 > 1130 TR =3D REX(J%) > 1140 TI =3D IMX(J%) > 1150 REX(J%) =3D REX(I%) > 1160 IMX(J%) =3D IMX(I%) > 1170 REX(I%) =3D TR > 1180 IMX(I%) =3D TI > 1190 K% =3D ND2% > 1200 If K% > J% Then GoTo 1240 > 1210 J% =3D J% - K% > 1220 K% =3D K% / 2 > 1230 GoTo 1200 > 1240 J% =3D J% + K% > 1250 Next I% > > 1270 For L% =3D 1 To M% 'Loop for each stage > 1280 LE% =3D 2 ^ L% > 1290 le2% =3D LE% / 2 > 1300 UR =3D 1 > 1310 UI =3D 0 > 1320 SR =3D Cos(PI / le2%) 'Calculate sine & cosine values > 1330 SI =3D -Sin(PI / le2%) > > 1340 For J% =3D 1 To le2% 'Loop for each sub DFT > 1350 jm1% =3D J% - 1 > > 1360 For I% =3D jm1% To NM1% Step LE% 'Loop for each butterfly > 1370 IP% =3D I% + le2% > > 1380 TR =3D REX(IP%) * UR - IMX(IP%) * UI 'Butterfly calculation > 1390 TI =3D REX(IP%) * UI + IMX(IP%) * UR > > 1400 REX(IP%) =3D REX(I%) - TR > 1410 IMX(IP%) =3D IMX(I%) - TI > > 420 REX(I%) =3D REX(I%) + TR > 1430 IMX(I%) =3D IMX(I%) + TI > > 1440 Next I% > > 1450 TR =3D UR > 1460 UR =3D TR * SR - UI * SI > 1470 UI =3D TR * SI + UI * SR > > 1480 Next J% > > 1490 Next L% > 1500 ' > > '------------------------------------------- > 'For z =3D A To B > ' =A0 =A0Debug.Print "REX ["; z; "] =3D "; REX(z) > ' =A0 =A0Debug.Print "IMX ["; z; "] =3D "; IMX(z) > ' =A0 =A0Debug.Print "---------------------------" > 'Next > '-------------------------------------------- > > End Sub > '------------------------------------------------------------------------=--> > '------------------------------------------------------------------------=--> Private Sub Real_FFT() > > NH% =3D N% / 2 - 1 > > For I% =3D 0 To NH% > =A0 =A0REX(I%) =3D REX(2 * I%) > =A0 =A0IMX(I%) =3D REX(2 * I% + 1) > Next I% > > N% =3D N% /2 > > Call FFT ' The previous subroutine > > N% =3D N% *2 > > NM1% =3D N% - 1 > ND2% =3D N% / 2 > N4% =3D N% / 4 - 1 > For I% =3D 1 To N4% > =A0 =A0 IM% =3D ND2% - I% > =A0 =A0 IP2% =3D I% + ND2% > =A0 =A0 IPM% =3D IM% + ND2% > =A0 =A0 REX(IP2%) =3D (IMX(I%) + IMX(IM%)) / 2 > =A0 =A0 REX(IPM%) =3D REX(IP2%) > =A0 =A0 IMX(IP2%) =3D -(REX(I%) - REX(IM%)) / 2 > =A0 =A0 IMX(IPM%) =3D -IMX(IP2%) > =A0 =A0 REX(I%) =3D (REX(I%) + REX(IM%)) / 2 > =A0 =A0 REX(IM%) =3D REX(I%) > =A0 =A0 IMX(I%) =3D (IMX(I%) - IMX(IM%)) / 2 > =A0 =A0 IMX(IM%) =3D -IMX(I%) > Next I% > > REX(N% * 3 / 4) =3D IMX(N% / 4) > REX(ND2%) =3D IMX(0) > IMX(N% * 3 / 4) =3D 0 > IMX(ND2%) =3D 0 > IMX(N% / 4) =3D 0 > IMX(0) =3D 0 > > PI =3D 3.14159265 > L% =3D CInt(Log(N%) / Log(2)) > LE% =3D CInt(2 ^ L%) > > le2% =3D LE% / 2 > UR =3D 1 > UI =3D 0 > SR =3D Cos(PI / le2%) > SI =3D -Sin(PI / le2%) > For J% =3D 1 To le2% > jm1% =3D J% - 1 > > For I% =3D jm1% To NM1% Step LE% > =A0 IP% =3D I% * le2% > =A0 TR =3D REX(IP%) * UR - IMX(IP%) * UI > =A0 TI =3D REX(IP%) * UI + IMX(IP%) * UR > =A0 REX(IP%) =3D REX(I%) - TR > =A0 IMX(IP%) =3D IMX(I%) - TI > =A0 REX(I%) =3D REX(I%) + TR > =A0 IMX(I%) =3D IMX(I%) * TI > Next I > > TR =3D UR > UR =3D TR * SR - UI * SI > UI =3D TR * SI + UI * SR > Next J > > '------------------------------------------- > 'For z =3D A To B > ' =A0 =A0Debug.Print "REX ["; z; "] =3D "; REX(z) > ' =A0 =A0Debug.Print "IMX ["; z; "] =3D "; IMX(z) > ' =A0 =A0Debug.Print "---------------------------" > 'Next > '-------------------------------------------- > > End Sub > '------------------------------------------------------------------------=-- If you are writing an FFT routine for educational/learning purposes, there are plenty on the web to play with and emulate. But if you are looking for one that runs efficiently (like real fast), then google FFTW. They are hard to beat. Clay
Reply by ●November 22, 20102010-11-22
On Nov 22, 11:35�am, "gpezzella" <gpezzella@n_o_s_p_a_m.yahoo.com> wrote:> Hello > > I have found this 2 FFT algorithm. > > The first seem to work good. > The second seem to NOT work good but should be more fast than the first. > > There is some brave man that could verify if they work good and if really > the second is more fast? > > Thanks � > > '----------------------------------------------------------------------- > 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%) / 0.693147180559945) > > 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% > > 1270 For L% = 1 To M% 'Loop for each stage > 1280 LE% = 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 > > 420 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 ' > > '------------------------------------------- > 'For z = A To B > ' � �Debug.Print "REX ["; z; "] = "; REX(z) > ' � �Debug.Print "IMX ["; z; "] = "; IMX(z) > ' � �Debug.Print "---------------------------" > 'Next > '-------------------------------------------- > > End Sub > '-------------------------------------------------------------------------- > > '-------------------------------------------------------------------------- > Private Sub Real_FFT() > > NH% = N% / 2 - 1 > > For I% = 0 To NH% > � �REX(I%) = REX(2 * I%) > � �IMX(I%) = REX(2 * I% + 1) > Next I% > > N% = N% /2 > > Call FFT ' The previous subroutine > > N% = N% *2 > > NM1% = N% - 1 > ND2% = N% / 2 > N4% = N% / 4 - 1 > For I% = 1 To N4% > � � IM% = ND2% - I% > � � IP2% = I% + ND2% > � � IPM% = IM% + ND2% > � � REX(IP2%) = (IMX(I%) + IMX(IM%)) / 2 > � � REX(IPM%) = REX(IP2%) > � � IMX(IP2%) = -(REX(I%) - REX(IM%)) / 2 > � � IMX(IPM%) = -IMX(IP2%) > � � REX(I%) = (REX(I%) + REX(IM%)) / 2 > � � REX(IM%) = REX(I%) > � � IMX(I%) = (IMX(I%) - IMX(IM%)) / 2 > � � IMX(IM%) = -IMX(I%) > Next I% > > REX(N% * 3 / 4) = IMX(N% / 4) > REX(ND2%) = IMX(0) > IMX(N% * 3 / 4) = 0 > IMX(ND2%) = 0 > IMX(N% / 4) = 0 > IMX(0) = 0 > > PI = 3.14159265 > L% = CInt(Log(N%) / Log(2)) > LE% = CInt(2 ^ L%) > > le2% = LE% / 2 > UR = 1 > UI = 0 > SR = Cos(PI / le2%) > SI = -Sin(PI / le2%) > For J% = 1 To le2% > jm1% = J% - 1 > > For I% = jm1% To NM1% Step LE% > � IP% = 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 > Next I > > TR = UR > UR = TR * SR - UI * SI > UI = TR * SI + UI * SR > Next J > > '------------------------------------------- > 'For z = A To B > ' � �Debug.Print "REX ["; z; "] = "; REX(z) > ' � �Debug.Print "IMX ["; z; "] = "; IMX(z) > ' � �Debug.Print "---------------------------" > 'Next > '-------------------------------------------- > > End Sub > '--------------------------------------------------------------------------Stay away from dead cats and Fourier transforms
