fftlib.bas
'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
|
|
搜索
论坛
Home
回顶部
g.yi.org Hits:0 Last modified:2010-12-07 22:25:00
