'DLL code for the FreeBasic compiler by V1ctor ' to compile --> fbc -dll fftlib.bas '-------------------------------------------------------------------- 'code from several sources: ' VB FFT Release 2-B by Murphy McCauley (MurphyMc@Concentric.NET) (removed that fft for errors) ' Stanford Meterology dept -- I think they derived their code from SETOOLS? option explicit option dynamic '$INCLUDE: 'crt.bi' '$INCLUDE: 'windows.bi' 'second protots DECLARE FUNCTION CheckPowerOfTwo lib "fftlib" alias "CheckPowerOfTwo"(X As Long) AS INTEGER DECLARE FUNCTION ExactBlackman lib "fftlib" alias "ExactBlackman"(n AS INTEGER, j AS INTEGER) AS SINGLE DECLARE SUB FFTMagnitude lib "fftlib" alias "FFTMagnitude"(lpxr AS SINGLE PTR, lpyi AS SINGLE PTR, BYVAL n AS INTEGER, lpOutVal AS SINGLE PTR) DECLARE SUB FFTPhase lib "fftlib" alias "FFTPhase"(lpxr AS SINGLE PTR, lpyi AS SINGLE PTR, BYVAL n AS INTEGER, lpPhase AS SINGLE PTR) DECLARE FUNCTION Hamming lib "fftlib" alias "Hamming"(n AS INTEGER, j AS INTEGER) AS SINGLE DECLARE FUNCTION Hanning lib "fftlib" alias "Hanning"(n AS INTEGER, j AS INTEGER) AS SINGLE DECLARE FUNCTION Parzen lib "fftlib" alias "Parzen"(n AS INTEGER, j AS INTEGER) AS SINGLE DECLARE FUNCTION SigmaGauss lib "fftlib" alias "SigmaGauss"(n AS INTEGER, j AS INTEGER) AS SINGLE DECLARE FUNCTION SquareAndSum lib "fftlib" alias "SquareAndSum"(BYVAL a AS SINGLE, BYVAL b AS SINGLE) AS SINGLE DECLARE FUNCTION Welch lib "fftlib" alias "Welch"(n AS INTEGER, j AS INTEGER) AS SINGLE DECLARE SUB FFTFrequency lib "fftlib" alias "FFTFrequency"(lpHz AS SINGLE PTR, BYVAL Sample_Hz AS SINGLE, BYVAL n AS INTEGER) DECLARE SUB Convolve lib "fftlib" alias "Convolve"(lpfiltcoef AS SINGLE PTR, lpIn AS SINGLE PTR, lpOut AS SINGLE PTR, BYVAL filtLen AS INTEGER, BYVAL aLen AS INTEGER) DECLARE SUB DFT lib "fftlib" alias "DFT"(lpReal AS SINGLE PTR, lpImag AS SINGLE PTR, BYVAL Sample_N AS INTEGER, BYVAL inverse AS INTEGER) DECLARE SUB DFTMagnitude lib "fftlib" alias "DFTMagnitude"(lpReal AS SINGLE PTR, lpImag AS SINGLE PTR, BYVAL Sample_N AS INTEGER, lpAmp AS SINGLE PTR) DECLARE SUB DFTPhase lib "fftlib" alias "DFTPhase"(lpReal AS SINGLE PTR, lpImag AS SINGLE PTR, BYVAL Sample_N AS INTEGER, lpPhase AS SINGLE PTR) DECLARE SUB FFT2DCalc lib "fftlib" alias "FFT2DCalc"(lpxreal AS SINGLE PTR, lpyImag AS SINGLE PTR, BYVAL c AS INTEGER, BYVAL r AS INTEGER, BYVAL flag AS INTEGER) DECLARE SUB FFT2DCRCalc (xreal() AS SINGLE, yimag() AS SINGLE, c AS INTEGER, r AS INTEGER, flag AS INTEGER) DECLARE SUB FFT2DRCCalc (xreal() AS SINGLE, yimag() AS SINGLE, c AS INTEGER, r AS INTEGER, flag AS INTEGER) DECLARE SUB PowerSpectrumCalc lib "fftlib" alias "PowerSpectrumCalc"(lpxreal AS SINGLE PTR, lpyimag AS SINGLE PTR, BYVAL numdat AS INTEGER, BYVAL delta AS SINGLE) DECLARE SUB RealFFT2 (xr() AS SINGLE, sinc() AS SINGLE, n AS INTEGER, inverse AS INTEGER) DECLARE SUB WindowData lib "fftlib" alias "WindowData"(lpx AS SINGLE PTR, numdat AS INTEGER, wind AS INTEGER) DECLARE SUB fft lib "fftlib" alias "fft"(lpxr AS SINGLE PTR, lpyi AS SINGLE PTR, BYVAL n AS INTEGER, BYVAL inverse AS INTEGER) DECLARE SUB fft2 (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, m AS INTEGER, ks AS INTEGER) DECLARE SUB reorder (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, m AS INTEGER, ks AS INTEGER, reel AS INTEGER) DECLARE SUB rfft2 (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, m AS INTEGER, ks AS INTEGER) DECLARE SUB trans (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, eval AS INTEGER) DECLARE Sub CalcFrequency lib "fftlib" alias "CalcFrequency"(BYVAL NumberOfSamples As Long, BYVAL FrequencyIndex As Long, lpRealIn AS SINGLE PTR, lpImagIn AS SINGLE PTR, BYREF RealOut AS SINGLE, BYREF ImagOut AS SINGLE) DECLARE Function NumberOfBitsNeeded(PowerOfTwo As Long) As Byte DECLARE Function ReverseBits(ByVal Index As Long, NumBits As Byte) As Long CONST fft.pi = 3.141592653589793# ' *** convolution of filtcoef (length filtLen) with input x() ' *** input x() of aLen => output y() SUB Convolve (lpfiltcoef AS SINGLE PTR, lpIn AS SINGLE PTR, lpOut AS SINGLE PTR, BYVAL filtLen AS INTEGER, BYVAL aLen AS INTEGER) EXPORT DIM n AS INTEGER DIM m AS INTEGER DIM summ AS INTEGER DIM fp As Single Ptr: fp = @lpfiltcoef : DIM filtcoef(filtLen) AS SINGLE : MEMCPY(@filtcoef(0), fp, filtLen * SIZEOF(SINGLE)) DIM ip As Single Ptr: ip = @lpIn: DIM Inx(aLen) AS SINGLE : MEMCPY(@Inx(0), ip, aLen * SIZEOF(SINGLE)) DIM op As Single Ptr: op = @lpOut: DIM Outy(aLen) AS SINGLE : MEMCPY(@Outy(0), op, aLen * SIZEOF(SINGLE)) FOR n = 0 TO (aLen + filtLen - 2) summ = 0! FOR m = 0 TO filtLen - 1 IF ((n - m) >= 0) AND ((n - m) <= (aLen - 1)) THEN summ += filtcoef(m) * Inx(n - m) END IF NEXT Outy(n) = summ NEXT MEMCPY(fp, @filtcoef(0), filtLen * SIZEOF(SINGLE)) MEMCPY(ip, @Inx(0), aLen * SIZEOF(SINGLE)) MEMCPY(op, @Outy(0), aLen * SIZEOF(SINGLE)) END SUB '********************************************************* 'sub section of calculations for forward and inverse FFT '********************************************************* SUB fft2 (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, m AS INTEGER, ks AS INTEGER) DIM cc(0 TO m) AS INTEGER DIM indx AS INTEGER, k0 AS INTEGER, k1 AS INTEGER, k2 AS INTEGER, k3 AS INTEGER DIM span AS INTEGER, j AS INTEGER, k AS INTEGER, kb AS INTEGER DIM kn AS INTEGER, fp AS INTEGER DIM mm AS INTEGER, mk AS INTEGER, breakf AS INTEGER, bf2 AS INTEGER DIM rad AS SINGLE, c1 AS SINGLE, c2 AS SINGLE DIM c3 AS SINGLE, s1 AS SINGLE, s2 AS SINGLE, s3 AS SINGLE, ck AS SINGLE, sk AS SINGLE, sq DIM x0 AS SINGLE, x1 AS SINGLE, x2 AS SINGLE, x3 AS SINGLE DIM y0 AS SINGLE, y1 AS SINGLE, y2 AS SINGLE, y3 sq = .707106781187# sk = .382683432366# ck = .92387953251# cc(m) = ks mm = (m \ 2) * 2 kn = 0 FOR k = m - 1 TO 0 STEP -1 cc(k) = cc(k + 1) \ 2 NEXT rad = 6.28318530718# / ((cc(0)) * (ks)) mk = m - 5 breakf = 0 WHILE ((kn < n) AND (breakf = 0)) 999 breakf = 0 kb = kn kn = kn + ks IF (mm <> m) THEN k2 = kn k0 = cc(mm) + kb WHILE (k0 > kb) k2 = k2 - 1 k0 = k0 - 1 x0 = xr(k2) y0 = y(k2) xr(k2) = xr(k0) - x0 xr(k0) = xr(k0) + x0 y(k2) = y(k0) - y0 y(k0) = y(k0) + y0 WEND END IF c1 = 1! s1 = 0! indx = 0 k = mm - 2 j = 3 IF (k >= 0) THEN GOTO 1002 ELSE IF (kn < n) THEN GOTO 999 ELSE breakf = 1 GOTO 1008 END IF END IF 1001 bf2 = 1 WHILE (bf2 = 1) IF (cc(j) <= indx) THEN indx = indx - cc(j) j = j - 1 IF (cc(j) <= indx) THEN indx = indx - cc(j) j = j - 1 k = k + 2 END IF ELSE bf2 = 0 END IF WEND indx = cc(j) + indx j = 3 1002 span = cc(k) fp = 0 IF (indx <> 0) THEN c2 = (indx) * (span) * rad c1 = COS(c2) s1 = SIN(c2) END IF 1005 IF ((fp = 1) OR (indx <> 0)) THEN c2 = c1 * c1 - s1 * s1 s2 = 2! * s1 * c1 c3 = c2 * c1 - s2 * s1 s3 = c2 * s1 + s2 * c1 END IF FOR k0 = kb + span - 1 TO kb STEP -1 k1 = k0 + span k2 = k1 + span k3 = k2 + span x0 = xr(k0) y0 = y(k0) IF (s1 = 0!) THEN x1 = xr(k1) y1 = y(k1) x2 = xr(k2) y2 = y(k2) x3 = xr(k3) y3 = y(k3) ELSE x1 = xr(k1) * c1 - y(k1) * s1 y1 = xr(k1) * s1 + y(k1) * c1 x2 = xr(k2) * c2 - y(k2) * s2 y2 = xr(k2) * s2 + y(k2) * c2 x3 = xr(k3) * c3 - y(k3) * s3 y3 = xr(k3) * s3 + y(k3) * c3 END IF xr(k0) = x0 + x2 + x1 + x3 y(k0) = y0 + y2 + y1 + y3 xr(k1) = x0 + x2 - x1 - x3 y(k1) = y0 + y2 - y1 - y3 xr(k2) = x0 - x2 - y1 + y3 y(k2) = y0 - y2 + x1 - x3 xr(k3) = x0 - x2 + y1 - y3 y(k3) = y0 - y2 - x1 + x3 NEXT IF (k > 0) THEN k = k - 2 GOTO 1002 END IF kb = k3 + span IF (kb < kn) THEN IF (j = 0) THEN k = 2 j = mk GOTO 1001 END IF j = j - 1 c2 = c1 IF (j = 1) THEN c1 = c1 * ck + s1 * sk s1 = s1 * ck - c2 * sk ELSE c1 = (c1 - s1) * sq s1 = (s1 + c2) * sq END IF fp = 1 GOTO 1005 END IF 1008 ' CONTINUE WEND END SUB 'need all arrays to be a power of 2, except DFT(,,,) FUNCTION CheckPowerOfTwo (X As Long) AS INTEGER EXPORT DIM chk as integer: chk = 0 If (X And (X - 1)) = 0 Then chk = 1 if x < 2 Then chk = 0 CheckPowerOfTwo = chk End FUNCTION '============================================================================================ 'THE DISCRETE FOURIER TRANSFORM by JohnK, use this only for small data sets 'code adapted from from: "The Scientist and Engineer's Guide to Digital Signal 'Processing (www.dspguide.com), which allows programs to be copied, distributed, ' and used for any noncommercial purpose. ' ***** the signal in does not have to be a power of 2 **** ' The frequency domain signals, held in Real and Imag components, are ' calculated from a time domain signal. ' RawDat( ) holds the time domain signal ' RealCom holds the real part of the frequency domain ' Imag holds the imaginary part of the frequency domain ' SUB DFT (lpReal AS SINGLE PTR, lpImag AS SINGLE PTR, BYVAL Sample_N AS INTEGER, BYVAL inverse AS INTEGER) EXPORT DIM Omega AS SINGLE: Omega = 2.0! * fft.pi / Sample_N DIM k AS INTEGER 'loops DIM i AS INTEGER 'loops DIM Rsin AS SINGLE DIM Isin AS SINGLE 'complex sinusoids DIM RealSum(0 TO Sample_N) AS SINGLE DIM ImagSum(0 TO Sample_N) AS SINGLE DIM rp As Single Ptr: rp = @lpReal : DIM Real(Sample_N) AS SINGLE : MEMCPY(@Real(0), rp, Sample_N * SIZEOF(SINGLE)) DIM ip As Single Ptr: ip = @lpImag : DIM Imag(Sample_N) AS SINGLE : MEMCPY(@Imag(0), ip, Sample_N * SIZEOF(SINGLE)) IF inverse = 0 THEN 'THE FORWARD DISCRETE FOURIER TRANSFORM FOR k = 0 TO (Sample_N/2) 'loops through each sample in Real and Imaginary component array FOR i = 0 TO Sample_N 'for all samples in output ( ) Rsin = COS(Omega * k * i) 'speed issues Isin = -SIN(Omega * k * i) RealSum(k) += Real(i) * Rsin ImagSum(k) += Real(i) * Isin IF Imag(i) <> 0! THEN 'don't waste computation for null array RealSum(k) += Imag(i) * Isin 'complex sinusoid ImagSum(k) += Imag(i) * Rsin END IF NEXT i NEXT k ELSE 'THE INVERSE DISCRETE FOURIER TRANSFORM 'The time domain signal, held in SigOut[ ] of size N, is calculated from the frequency 'domain signals, held in the real (real[ ]) and imaginary frequency domains (imag[ ]) FOR k = 0 TO Sample_N/2 'Find the cosine and sine wave amplitudes Real(k) = Real(k) / (Sample_N/2) Imag(k) = -Imag(k) / (Sample_N/2) NEXT k Real(0) = Real(0) / 2 Real(Sample_N/2) = Real(Sample_N/2) / 2 'SYNTHESIS METHOD #1 Loop through each frequency generating the entire length of the sine & cosine 'waves, and add them to the accumulator OUTPUT signal FOR k = 0 TO Sample_N\2 'loops through each sample in Real( ) and Imag( ) FOR i = 0 TO (Sample_N-1) 'loops through each sample in Output -- RealSum( ) RealSum(i) += Real(k) * COS(Omega * k * i) RealSum(i) += Imag(k) * SIN(Omega * k * i) NEXT i NEXT k END IF 'FORWARD or INVERSE MEMCPY(rp, @RealSum(0), Sample_N/2 * SIZEOF(SINGLE)) MEMCPY(ip, @ImagSum(0), Sample_N/2 * SIZEOF(SINGLE)) END SUB SUB DFTMagnitude (lpReal AS SINGLE PTR, lpImag AS SINGLE PTR, BYVAL Sample_N AS INTEGER, lpAmp AS SINGLE PTR) EXPORT DIM k AS INTEGER 'loops DIM nm AS SINGLE: nm = 2/Sample_N DIM rp As Single Ptr: rp = @lpReal : DIM Real(Sample_N) AS SINGLE : MEMCPY(@Real(0), rp, Sample_N * SIZEOF(SINGLE)) DIM ip As Single Ptr: ip = @lpImag : DIM Imag(Sample_N) AS SINGLE : MEMCPY(@Imag(0), ip, Sample_N * SIZEOF(SINGLE)) DIM ap As Single Ptr: ap = @lpAmp : DIM Amp(Sample_N) AS SINGLE : MEMCPY(@Amp(0), ap, Sample_N * SIZEOF(SINGLE)) FOR k = 0 TO (Sample_N/2) ' Amp(k) = 2 * SQR(Imag(k) * Imag(k) + Real(k) * Real(k))/ Sample_N Amp(k) = nm * SQR(Imag(k) * Imag(k) + Real(k) * Real(k)) NEXT k MEMCPY(ap, @Amp(0), Sample_N * SIZEOF(SINGLE)) END SUB SUB DFTPhase (lpReal AS SINGLE PTR, lpImag AS SINGLE PTR, BYVAL Sample_N AS INTEGER, lpPhase AS SINGLE PTR) EXPORT DIM k AS INTEGER 'loops DIM rp As Single Ptr: rp = @lpReal : DIM Real(Sample_N) AS SINGLE : MEMCPY(@Real(0), rp, Sample_N * SIZEOF(SINGLE)) DIM ip As Single Ptr: ip = @lpImag : DIM Imag(Sample_N) AS SINGLE : MEMCPY(@Imag(0), ip, Sample_N * SIZEOF(SINGLE)) DIM pp As Single Ptr: pp = @lpPhase :DIM Phase(Sample_N) AS SINGLE : MEMCPY(@Phase(0), pp, Sample_N * SIZEOF(SINGLE)) DIM RealComp AS SINGLE 'store results for convenience of phase calc DIM ImagComp AS SINGLE FOR k = 0 TO (Sample_N/2) RealComp = Real(k) 'store results for convenience of phase calc ImagComp = Imag(k) 'Compute the phase coefficient IF (ImagComp > 0) AND (RealComp > 0) THEN Phase(k) = ATN(ImagComp/RealComp) IF (ImagComp > 0) AND (RealComp < 0) THEN Phase(k) = fft.pi - ATN(-ImagComp/RealComp) IF (ImagComp < 0) AND (RealComp > 0) THEN Phase(k) = -1 * ATN(-ImagComp/RealComp) IF (ImagComp < 0) AND (RealComp < 0) THEN Phase(k) = -fft.pi + ATN(ImagComp/RealComp) IF (ImagComp = 0) AND (RealComp >= 0) THEN Phase(k) = 0.00 IF (ImagComp = 0) AND (RealComp < 0) THEN Phase(k) = fft.pi IF (RealComp = 0) AND (ImagComp > 0) THEN Phase(k) = fft.pi / 2 IF (RealComp = 0) AND (ImagComp < 0) THEN Phase(k) = -fft.pi / 2 Phase(k) = Phase(k) * 180.0 / fft.pi 'convert to degrees NEXT k MEMCPY(pp, @Phase(0), Sample_N * SIZEOF(SINGLE)) END SUB SUB FFT2DCalc (lpxreal AS SINGLE PTR, lpyImag AS SINGLE PTR, BYVAL c AS INTEGER, BYVAL r AS INTEGER, BYVAL flag AS INTEGER) EXPORT DIM rp As Single Ptr: rp = @lpxReal : DIM xReal(r, c) AS SINGLE : MEMCPY(@xReal(0, 0), rp, r * c * SIZEOF(SINGLE)) DIM ip As Single Ptr: ip = @lpyImag : DIM yImag(r, c) AS SINGLE : MEMCPY(@yImag(0, 0), ip, r * c * SIZEOF(SINGLE)) IF CheckPowerOfTwo(r) AND CheckPowerOfTwo(c) THEN 'make sure Radix-2 IF (flag = 0) THEN FFT2DRCCalc(xReal(), yImag(), c, r, flag) 'real to complex ELSE FFT2DCRCalc(xReal(), yImag(), c, r, flag) 'complex to real END IF MEMCPY(rp, @xReal(0, 0), r * c * SIZEOF(SINGLE)) MEMCPY(ip, @yImag(0, 0), r * c * SIZEOF(SINGLE)) END IF END SUB SUB FFT2DCRCalc (xreal() AS SINGLE, yimag() AS SINGLE, c AS INTEGER, r AS INTEGER, flag AS INTEGER) DIM i AS INTEGER, a AS INTEGER, j AS INTEGER DIM xVector(c - 1) AS SINGLE DIM yVector(c - 1) AS SINGLE FOR j = 0 TO r - 1 FOR a = 0 TO c - 1 xVector(a) = xreal(j, a) yVector(a) = yimag(j, a) NEXT IF (flag = 0) THEN 'forward FFT fft(@xVector(0), @yVector(0), c, 0) ELSE 'inverse FFT fft(@xVector(0), @yVector(0), c, 1) END IF FOR a = 0 TO c - 1 xreal(j, a) = xVector(a) yimag(j, a) = yVector(a) NEXT NEXT FOR i = 0 TO c - 1 FOR a = 0 TO r - 1 xVector(a) = xreal(a, i) yVector(a) = yimag(a, i) NEXT IF (flag = 0) THEN 'forward fft fft(@xVector(0), @yVector(0), r, 0) ELSE 'inverse fft fft(@xVector(0), @yVector(0), r, 1) END IF FOR a = 0 TO r - 1 xreal(a, i) = xVector(a) yimag(a, i) = yVector(a) NEXT NEXT END SUB '2d fast fourier transform from real to complex ' xreal() is a real 2D array of size(r,c) Radix-2 ' yreal() is an imaginary 2D array of size(r,c) Radix-2 SUB FFT2DRCCalc (xreal() AS SINGLE, yimag() AS SINGLE, c AS INTEGER, r AS INTEGER, flag AS INTEGER) DIM i AS INTEGER, a AS INTEGER, j AS INTEGER DIM xVector(r - 1) AS SINGLE, yVector(r - 1) AS SINGLE FOR j = 0 TO c - 1 FOR a = 0 TO r - 1 xVector(a) = xreal(a, j) yVector(a) = yimag(a, j) NEXT a IF (flag = 0) THEN 'forward fft fft(@xVector(0), @yVector(0), c, 0) ELSE 'inverse fft(@xVector(0), @yVector(0), c, 1) END IF FOR a = 0 TO r - 1 xreal(a, j) = xVector(a) yimag(a, j) = yVector(a) NEXT a NEXT j FOR i = 0 TO r - 1 FOR a = 0 TO c - 1 xVector(a) = xreal(i, a) yVector(a) = yimag(i, a) NEXT a IF (flag = 0) THEN 'forward fft fft(@xVector(0), @yVector(0), r, 0) ELSE 'inverse fft fft(@xVector(0), @yVector(0), r, 1) END IF FOR a = 0 TO c - 1 xreal(i, a) = xVector(a) yimag(i, a) = yVector(a) NEXT a NEXT i END SUB SUB FFTFrequency (lpHz AS SINGLE PTR, BYVAL Sample_Hz AS SINGLE, BYVAL n AS INTEGER) EXPORT DIM k AS INTEGER DIM hzp As Single Ptr: hzp = @lpHz : DIM Hz(n) AS SINGLE : MEMCPY(@Hz(0), hzp, n * SIZEOF(SINGLE)) FOR k = 0 TO (n/2) Hz(k) = k * Sample_Hz / n NEXT k FOR k = (n/2 + 1) TO n Hz(k) = -Sample_Hz * (n-k) / n NEXT k MEMCPY(hzp, @Hz(0), n * SIZEOF(SINGLE)) END SUB ' CalcFrequency() transforms a single frequency. Sub CalcFrequency (BYVAL NumberOfSamples As Long, BYVAL FrequencyIndex As Long, lpRealIn AS SINGLE PTR, lpImagIn AS SINGLE PTR, BYREF RealOut AS SINGLE, BYREF ImagOut AS SINGLE) EXPORT Dim K As Long Dim Cos1 AS SINGLE, Cos2 AS SINGLE, Cos3 AS SINGLE, Theta AS SINGLE, Beta AS SINGLE Dim Sin1 AS SINGLE, Sin2 AS SINGLE, Sin3 AS SINGLE DIM rp As Single Ptr: rp = @lpRealIn : DIM RealIn(NumberOfSamples) AS SINGLE : MEMCPY(@RealIn(0), rp, NumberOfSamples * SIZEOF(SINGLE)) DIM ip As Single Ptr: ip = @lpImagIn : DIM ImagIn(NumberOfSamples) AS SINGLE : MEMCPY(@ImagIn(0), ip, NumberOfSamples * SIZEOF(SINGLE)) Theta = 2 * fft.pi * FrequencyIndex / NumberOfSamples Sin1 = Sin(-2 * Theta) Sin2 = Sin(-Theta) Cos1 = Cos(-2 * Theta) Cos2 = Cos(-Theta) Beta = 2 * Cos2 For K = 0 To NumberOfSamples - 2 'Update trig values Sin3 = Beta * Sin2 - Sin1 Sin1 = Sin2 Sin2 = Sin3 Cos3 = Beta * Cos2 - Cos1 Cos1 = Cos2 Cos2 = Cos3 RealOut = RealOut + RealIn(K) * Cos3 - ImagIn(K) * Sin3 ImagOut = ImagOut + ImagIn(K) * Cos3 + RealIn(K) * Sin3 Next k MEMCPY(rp, @RealIn(0), NumberOfSamples * SIZEOF(SINGLE)) MEMCPY(ip, @ImagIn(0), NumberOfSamples * SIZEOF(SINGLE)) End Sub SUB FFTMagnitude (lpxr AS SINGLE PTR, lpyi AS SINGLE PTR, BYVAL n AS INTEGER, lpOutVal AS SINGLE PTR) EXPORT DIM i AS INTEGER DIM rr AS SINGLE, ii AS SINGLE DIM xp As Single Ptr: xp = @lpxr : DIM xr(n) AS SINGLE : MEMCPY(@xr(0), xp, n * SIZEOF(SINGLE)) DIM yp As Single Ptr: yp = @lpyi : DIM yi(n) AS SINGLE : MEMCPY(@yi(0), yp, n * SIZEOF(SINGLE)) DIM op As Single Ptr: op = @lpOutVal : DIM OutVal(n) AS SINGLE : MEMCPY(@OutVal(0), op, n * SIZEOF(SINGLE)) OutVal(0) = SQR(xr(0) * xr(0) + yi(0) * yi(0)) / n FOR i = 1 TO n rr = ABS(xr(i)) + ABS(xr(n - i)) ii = ABS(yi(i)) + ABS(yi(n - i)) OutVal(i) = SQR(rr * rr + ii * ii) / n NEXT i MEMCPY(op, @OutVal(0), n * SIZEOF(SINGLE)) END SUB SUB FFTPhase (lpxr AS SINGLE PTR, lpyi AS SINGLE PTR, BYVAL n AS INTEGER, lpPhase AS SINGLE PTR) EXPORT DIM xp As Single Ptr: xp = @lpxr : DIM xr(n) AS SINGLE : MEMCPY(@xr(0), xp, n * SIZEOF(SINGLE)) DIM yp As Single Ptr: yp = @lpyi : DIM yi(n) AS SINGLE : MEMCPY(@yi(0), yp, n * SIZEOF(SINGLE)) DIM pp As Single Ptr: pp = @lpPhase : DIM Phase(n) AS SINGLE : MEMCPY(@Phase(0), pp, n * SIZEOF(SINGLE)) DIM NearZero AS SINGLE: NearZero = 1E-20 DIM TempAngle AS SINGLE DIM i AS INTEGER FOR i = 0 TO n IF ABS(xr(i)) < NearZero THEN IF ABS(yi(i)) < NearZero THEN TempAngle = 0! IF yi(i) < -(NearZero) THEN TempAngle = 3! * fft.pi / 2! IF yi(i) > NearZero THEN TempAngle = fft.pi / 2! Phase(i) = TempAngle ' EXIT FOR ELSE IF ABS(yi(i)) < NearZero THEN IF (xr(i) < -(NearZero)) THEN TempAngle = fft.pi IF xr(i) > NearZero THEN TempAngle = 0 Phase(i) = TempAngle ' EXIT FOR ELSE IF (xr(i) > NearZero) AND (yi(i) > NearZero) THEN TempAngle = ATN(yi(i) / xr(i)) IF (xr(i) < -NearZero) AND (yi(i) > NearZero) THEN TempAngle = fft.pi - ATN(-yi(i) / xr(i)) IF (xr(i) < -NearZero) AND (yi(i) < -NearZero) THEN TempAngle = fft.pi + ATN(yi(i) / xr(i)) IF (xr(i) > NearZero) AND (yi(i) < -NearZero) THEN TempAngle = 2! * fft.pi - ATN(-yi(i) / xr(i)) END IF END IF Phase(i) = TempAngle NEXT i MEMCPY(pp, @Phase(0), n * SIZEOF(SINGLE)) END SUB '*************** power spectrum ************** ' input real data (xreal) ==> output power ' input imag data (yimag) ==> output each Frequency, or Hz SUB PowerSpectrumCalc (lpxreal AS SINGLE PTR, lpyimag AS SINGLE PTR, BYVAL numdat AS INTEGER, BYVAL delta AS SINGLE) EXPORT DIM numdatr AS SINGLE, normal AS SINGLE, timespan AS SINGLE DIM i AS INTEGER DIM rp As Single Ptr: rp = @lpxReal : DIM xReal(numdat) AS SINGLE : MEMCPY(@xReal(0), rp, numdat * SIZEOF(SINGLE)) DIM ip As Single Ptr: ip = @lpyImag : DIM yImag(numdat) AS SINGLE : MEMCPY(@yImag(0), ip, numdat * SIZEOF(SINGLE)) IF CheckPowerOfTwo(numdat) = 0 THEN EXIT SUB numdatr = numdat * 1! normal = 1! / numdatr ^ 2 fft(@xreal(0), @yimag(0), numdat, 0) 'perform the fft 'now set up the arrays for power and frequency output xreal(0) = SquareAndSum(xreal(0), yimag(0)) * normal FOR i = 1 TO (numdat \ 2) - 1 xreal(i) = normal * (SquareAndSum(xreal(i), yimag(i)) + SquareAndSum(xreal(numdat - i), yimag(numdat - i))) NEXT i = numdat \ 2 xreal(i) = SquareAndSum(xreal(i), yimag(i)) * normal timespan = numdat * delta FOR i = 0 TO (numdat \ 2) yimag(i) = INT(i) / timespan NEXT MEMCPY(rp, @xReal(0), numdat * SIZEOF(SINGLE)) MEMCPY(ip, @yImag(0), numdat * SIZEOF(SINGLE)) END SUB SUB RealFFT2 (xr() AS SINGLE, sinc() AS SINGLE, n AS INTEGER, inverse AS INTEGER) DIM nn AS INTEGER, j AS INTEGER, m AS INTEGER DIM pp IF CheckPowerOfTwo(n) = 0 THEN EXIT SUB m = -1 nn = n \ 2 ' determine power of 2 WHILE (nn > 0) nn = nn \ 2 m = m + 1 WEND nn = n \ 2 IF (inverse = 1) THEN pp = .5 / (nn) trans(xr(), sinc(), nn, 1) fft2(xr(), sinc(), nn, m, nn) FOR j = 0 TO nn - 1 xr(j) = xr(j) * pp sinc(j) = sinc(j) * (-pp) NEXT reorder(xr(), sinc(), nn, m, nn, 1) FOR j = 0 TO nn - 1 ' combine xr(j + nn) = sinc(j) NEXT ELSE FOR j = 0 TO nn - 1 ' split array in 2 sinc(j) = xr(j + nn) NEXT reorder(xr(), sinc(), nn, m, nn, 1) rfft2(xr(), sinc(), nn, m, 1) FOR j = 0 TO nn - 1 xr(j) = xr(j) * .5 sinc(j) = sinc(j) * .5 NEXT trans(xr(), sinc(), nn, 0) FOR j = 0 TO nn - 1 sinc(j) = -sinc(j) NEXT FOR j = 1 TO nn - 1 xr(j + nn) = xr(nn - j) sinc(j + nn) = -sinc(nn - j) NEXT END IF END SUB SUB reorder (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, m AS INTEGER, ks AS INTEGER, reel AS INTEGER) DIM i AS INTEGER, j AS INTEGER, indx AS INTEGER, k AS INTEGER, kk AS INTEGER, kb AS INTEGER, k2 AS INTEGER, KU AS INTEGER, lim AS INTEGER, p AS INTEGER DIM breakf AS INTEGER DIM temp DIM cc(m) AS INTEGER DIM lst(m) cc(m) = ks FOR k = m TO 1 STEP -1 cc(k - 1) = cc(k) \ 2 NEXT p = m - 1 j = m - 1 i = 0 kb = 0 IF (reel <> 0) THEN KU = n - 2 FOR k = 0 TO KU STEP 2 temp = xr(k + 1) xr(k + 1) = y(k) y(k) = temp NEXT ELSE m = m - 1 END IF lim = (m + 2) \ 2 breakf = 0 IF (p > 0) THEN WHILE (breakf = 0) 3000 KU = cc(j) + kb k2 = KU indx = cc(m - j) kk = kb + indx WHILE (kk < KU) k = kk + indx WHILE (kk < k) temp = xr(kk) xr(kk) = xr(k2) xr(k2) = temp temp = y(kk) y(kk) = y(k2) y(k2) = temp kk = kk + 1 k2 = k2 + 1 WEND 'kk < k kk = kk + indx k2 = k2 + indx WEND ' KK < KU IF (j > lim) THEN j = j - 1 i = i + 1 lst(i) = j GOTO 3000 END IF kb = k2 IF (i > 0) THEN j = lst(i) i = i - 1 GOTO 3000 END IF IF (kb < n) THEN j = p ELSE breakf = 1 GOTO 3010 END IF 3010 ' CONTINUE WEND ' True END IF END SUB SUB rfft2 (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, m AS INTEGER, ks AS INTEGER) DIM k0 AS INTEGER, k1 AS INTEGER, k2 AS INTEGER, k3 AS INTEGER, k4 AS INTEGER, span AS INTEGER, j AS INTEGER, indx AS INTEGER, k AS INTEGER, kb AS INTEGER, nt AS INTEGER, kn AS INTEGER, mk AS INTEGER DIM Breakf1 AS INTEGER, fp AS INTEGER DIM rad AS SINGLE, c1 AS SINGLE, c2 AS SINGLE, c3 AS SINGLE, s1 AS SINGLE DIM s2 AS SINGLE, s3 AS SINGLE, ck AS SINGLE, sk AS SINGLE, sq AS SINGLE DIM x0 AS SINGLE, x1 AS SINGLE, x2 AS SINGLE, x3 AS SINGLE DIM y0 AS SINGLE, y1 AS SINGLE, y2 AS SINGLE, y3 AS SINGLE, re AS SINGLE, im AS SINGLE DIM cc(m) AS INTEGER sq = .707106781187# sk = .382683432366# ck = .92387953251# cc(0) = ks kn = 0 k4 = 4 * ks mk = m - 4 FOR k = 1 TO m ks = ks + ks cc(k) = ks NEXT rad = 3.14159265359# / ((cc(0)) * (ks)) WHILE (kn < n) kb = kn + 4 kn = kn + ks IF (m <> 1) THEN k = 0 indx = 0 j = mk nt = 3 c1 = 1! s1 = 0! Breakf1 = 0 WHILE (Breakf1 = 0) 2000 span = cc(k) fp = 0 IF (indx <> 0) THEN c2 = (indx) * (span) * rad c1 = COS(c2) s1 = SIN(c2) END IF 2001 IF ((fp = 1) OR (indx <> 0)) THEN c2 = c1 * c1 - s1 * s1 s2 = 2! * s1 * c1 c3 = c2 * c1 - s2 * s1 s3 = c2 * s1 + s2 * c1 ELSE s1 = 0! END IF k3 = kb - span WHILE (k3 < kb) k2 = k3 - span k1 = k2 - span k0 = k1 - span x0 = xr(k0) y0 = y(k0) x1 = xr(k1) y1 = y(k1) x2 = xr(k2) y2 = y(k2) x3 = xr(k3) y3 = y(k3) xr(k0) = x0 + x2 + x1 + x3 y(k0) = y0 + y2 + y1 + y3 IF (s1 = 0!) THEN xr(k1) = x0 - x1 - y2 + y3 y(k1) = y0 - y1 + x2 - x3 xr(k2) = x0 + x1 - x2 - x3 y(k2) = y0 + y1 - y2 - y3 xr(k3) = x0 - x1 + y2 - y3 y(k3) = y0 - y1 - x2 + x3 ELSE re = x0 - x1 - y2 + y3 im = y0 - y1 + x2 - x3 xr(k1) = re * c1 - im * s1 y(k1) = re * s1 + im * c1 re = x0 + x1 - x2 - x3 im = y0 + y1 - y2 - y3 xr(k2) = re * c2 - im * s2 y(k2) = re * s2 + im * c2 re = x0 - x1 + y2 - y3 im = y0 - y1 - x2 + x3 xr(k3) = re * c3 - im * s3 y(k3) = re * s3 + im * c3 END IF 's1 = 0 k3 = k3 + 1 WEND 'while (k3 < kb) nt = nt - 1 IF (nt >= 0) THEN c2 = c1 IF (nt = 1) THEN c1 = c1 * ck + s1 * sk s1 = s1 * ck - c2 * sk ELSE c1 = (c1 - s1) * sq s1 = (s1 + c2) * sq END IF kb = kb + k4 IF (kb <= kn) THEN fp = 1 GOTO 2001 ELSE Breakf1 = 1 GOTO 2005 END IF END IF 'nt >= 0 IF (nt = -1) THEN k = 2 GOTO 2000 END IF IF (cc(j) <= indx) THEN indx = indx - cc(j) j = j - 1 IF (cc(j) <= indx) THEN indx = indx - cc(j) j = j - 1 k = k + 2 ELSE indx = cc(j) + indx j = mk END IF ELSE indx = cc(j) + indx j = mk END IF 'cc[j] <= indx IF (j < mk) THEN GOTO 2000 k = 0 nt = 3 kb = kb + k4 IF (kb > kn) THEN Breakf1 = 1 2005 ' CONTINUE WEND ' while true END IF 'm <> 1 k = (m \ 2) * 2 IF (k <> m) THEN k2 = kn k0 = kn - cc(k) j = k0 WHILE (k2 > j) k2 = k2 - 1 k0 = k0 - 1 x0 = xr(k2) y0 = y(k2) xr(k2) = xr(k0) - x0 xr(k0) = xr(k0) + x0 y(k2) = y(k0) - y0 y(k0) = y(k0) + y0 WEND 'while (k2 > j) END IF ' k <> m WEND 'while (kn < n) END SUB FUNCTION SquareAndSum (BYVAL a AS SINGLE, BYVAL b AS SINGLE) AS SINGLE EXPORT SquareAndSum = a ^ 2 + b ^ 2 END FUNCTION SUB trans (xr() AS SINGLE, y() AS SINGLE, n AS INTEGER, eval AS INTEGER) DIM k AS INTEGER, nk AS INTEGER, nn AS INTEGER DIM x1 AS SINGLE, ab AS SINGLE, ba AS SINGLE, y1 AS SINGLE DIM re AS SINGLE, im AS SINGLE, ck AS SINGLE, sk AS SINGLE, dc AS SINGLE, ds AS SINGLE, r AS SINGLE nn = n \ 2 r = 3.14159265359# / (n) ds = SIN(r) r = 2! * SIN(.5 * r) r = -r * r dc = -.5 * r ck = 1! sk = 0! IF (eval = 0) THEN xr(n) = xr(0) y(n) = y(0) END IF FOR k = 0 TO nn nk = n - k x1 = xr(k) + xr(nk) ab = xr(k) - xr(nk) ba = y(k) + y(nk) y1 = y(k) - y(nk) re = ck * ba + sk * ab im = sk * ba - ck * ab y(nk) = im - y1 y(k) = im + y1 xr(nk) = x1 - re xr(k) = x1 + re dc = r * ck + dc ck = ck + dc ds = r * sk + ds sk = sk + ds NEXT END SUB SUB WindowData (lpx AS SINGLE PTR, numdat AS INTEGER, wind AS INTEGER) EXPORT DIM multiplier AS SINGLE, i AS INTEGER DIM xp As Single Ptr: xp = @lpx : DIM x(numdat)AS SINGLE : MEMCPY(@x(0), xp, numdat * SIZEOF(SINGLE)) FOR i = 0 TO numdat - 1 SELECT CASE wind CASE 0 multiplier = 1 CASE 1 multiplier = Parzen(numdat, i) CASE 2 multiplier = Hanning(numdat, i) CASE 3 multiplier = Welch(numdat, i) CASE 4 multiplier = Hamming(numdat, i) CASE 5 multiplier = ExactBlackman(numdat, i) CASE ELSE multiplier = 1! END SELECT x(i) = multiplier * x(i) NEXT i MEMCPY(xp, @x(0), numdat * SIZEOF(SINGLE)) END SUB ' **************** window filters *********************** FUNCTION ExactBlackman (n AS INTEGER, j AS INTEGER) AS SINGLE EXPORT ExactBlackman = .42 - .5 * COS(2! * fft.pi * j / (n - 1)) + .08 * COS(4! * fft.pi * j / (n - 1)) END FUNCTION FUNCTION Hamming (n AS INTEGER, j AS INTEGER) AS SINGLE EXPORT Hamming = .54 - .46 * COS(2! * fft.pi * j / (n - 1)) END FUNCTION FUNCTION Hanning (n AS INTEGER, j AS INTEGER) AS SINGLE EXPORT Hanning = .5 * (1 - COS(2! * fft.pi * j / (n - 1!))) END FUNCTION FUNCTION Parzen (n AS INTEGER, j AS INTEGER) AS SINGLE EXPORT Parzen = 1 - ABS((j - .5 * (n - 1)) / (.5 * (n - 1))) END FUNCTION FUNCTION SigmaGauss (n AS INTEGER, j AS INTEGER) AS SINGLE EXPORT SigmaGauss = Exp(-4.5 * (2 * j / n - 1) ^ 2) END FUNCTION FUNCTION Welch (n AS INTEGER, j AS INTEGER) AS SINGLE EXPORT Welch = 1 - ((j - .5 * (n - 1)) / (.5 * (n + 1))) ^ 2 END FUNCTION '-------------------------------------------------------------------- ' VB FFT Release 2-B ' by Murphy McCauley (MurphyMc@Concentric.NET) ' 10/01/99 '-------------------------------------------------------------------- ' About: ' This code is very, very heavily based on Don Cross's fourier.pas ' Turbo Pascal Unit for calculating the Fast Fourier Transform. ' I've not implemented all of his functions, though I may well do ' so in the future. ' For more info, you can contact me by email, check my website at: ' http://www.fullspectrum.com/deeth/ ' or check Don Cross's FFT web page at: ' http://www.intersrv.com/~dcross/fft.html ' You also may be intrested in the FFT.DLL that I put together based ' on Don Cross's FFT C code. It's callable with Visual Basic and ' includes VB declares. You can get it from either website. '-------------------------------------------------------------------- ' History of Release 2-B: ' Fixed a couple of errors that resulted from me mucking about with ' variable names after implementation and not re-checking. BAD ME. ' -------- ' History of Release 2: ' Added FrequencyOfIndex() which is Don Cross's Index_to_frequency(). ' FourierTransform() can now do inverse transforms. ' Added CalcFrequency() which can do a transform for a single ' frequency. '-------------------------------------------------------------------- ' Usage: ' The useful functions are: ' FourierTransform() performs a Fast Fourier Transform on an pair of ' Double arrays -- one real, one imaginary. Don't want/need ' imaginary numbers? Just use an array of 0s. This function can ' also do inverse FFTs. ' FrequencyOfIndex() can tell you what actual frequency a given index ' corresponds to. ' CalcFrequency() transforms a single frequency. '-------------------------------------------------------------------- ' Notes: ' All arrays must be 0 based (i.e. Dim TheArray(0 To 1023) or ' Dim TheArray(1023)). ' The number of samples must be a power of two (i.e. 2^x). ' FrequencyOfIndex() and CalcFrequency() haven't been tested much. ' Use this ENTIRELY AT YOUR OWN RISK. '-------------------------------------------------------------------- Function NumberOfBitsNeeded(PowerOfTwo As Long) As Byte Dim I As Byte For I = 0 To 16 If (PowerOfTwo And (2 ^ I)) <> 0 Then NumberOfBitsNeeded = I Exit Function End If Next End Function Function ReverseBits(ByVal Index As Long, NumBits As Byte) As Long Dim I As Byte, Rev As Long For I = 0 To NumBits - 1 Rev = (Rev * 2) Or (Index And 1) Index = Index \ 2 Next ReverseBits = Rev End Function '****************************************************************** ' Fourier transform : ' RealIn() real part in --> real part out ' yi() imaginary part in --> imaginary part out ' n sample size ' inverse -- if true then perform inverse Fourier transform '****************************************************************** SUB fft (lpRealIn AS SINGLE PTR, lpImagIn AS SINGLE PTR, BYVAL NumSamples AS INTEGER, BYVAL inverse AS INTEGER) EXPORT IF CheckPowerOfTwo(NumSamples) = 0 THEN EXIT SUB 'can't do Radix-N 'must transfer the data at the pointer to an internal array ' DIM rp As Single Ptr: rp = @lpRealIn : DIM RealIn(NumSamples) AS SINGLE : MEMCPY(@RealIn(0), rp, NumSamples * SIZEOF(SINGLE)) DIM RealIn As Single PTR: RealIn = @lpRealIn DIM ip As Single Ptr: ip = @lpImagIn : DIM ImagIn(NumSamples) AS SINGLE : MEMCPY(@ImagIn(0), ip, NumSamples * SIZEOF(SINGLE)) DIM RealOut(NumSamples) as single 'these hold temp output of real / imag DIM ImagOut(NumSamples) as single Dim AngleNumerator as single Dim NumBits As Byte, I As Long, j As Long, K As Long, n As Long, BlockSize As Long, BlockEnd As Long Dim DeltaAngle as single, DeltaAr as single Dim Alpha as single, Beta as single Dim TR as single, TI as single, AR as single, AI as single IF Inverse > 0 THEN AngleNumerator = -2# * fft.pi ELSE AngleNumerator = 2# * fft.pi END IF NumBits = NumberOfBitsNeeded(NumSamples) For I = 0 To (NumSamples - 1) j = ReverseBits(I, NumBits) RealOut(j) = RealIn[I] 'RealIn(I) ' ImagOut(j) = ImagIn(I) Next BlockEnd = 1 BlockSize = 2 Do While BlockSize <= NumSamples DeltaAngle = AngleNumerator / BlockSize Alpha = Sin(0.5 * DeltaAngle) Alpha = 2# * Alpha * Alpha Beta = Sin(DeltaAngle) I = 0 Do While I < NumSamples AR = 1# AI = 0# j = I For n = 0 To BlockEnd - 1 K = j + BlockEnd TR = AR * RealOut(K) - AI * ImagOut(K) TI = AI * RealOut(K) + AR * ImagOut(K) RealOut(K) = RealOut(j) - TR ImagOut(K) = ImagOut(j) - TI RealOut(j) += TR ImagOut(j) += TI DeltaAr = Alpha * AR + Beta * AI AI -= (Alpha * AI - Beta * AR) AR -= DeltaAr j += 1 Next n I += BlockSize Loop BlockEnd = BlockSize BlockSize = BlockSize * 2 Loop If Inverse Then 'Normalize the resulting time samples... For I = 0 To NumSamples - 1 RealOut(I) /= NumSamples ImagOut(I) /= NumSamples Next I End If ' MEMCPY(rp, @RealOut(0), NumSamples * SIZEOF(SINGLE)) 'restore the arrays to calling program MEMCPY(RealIn, @RealOut(0), NumSamples * SIZEOF(SINGLE)) 'restore the arrays to calling program MEMCPY(ip, @ImagOut(0), NumSamples * SIZEOF(SINGLE)) End Sub