Attribute VB_Name = "Module2" Dim A(1024, 2) As Double Dim ip As Integer, ndata As Integer '*************************************************************** '* Program to demonstrate the Fast Fourier Transform Procedure * '* (frequency analysis of a discreet signal F(t)) * '* * '* Visual Basic 4,0 release with graph by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ----------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Input data file (tgft.dat): * '* * '* (The test signal contains 3 frequencies: 50, 250, 500 hz) * '* * '* 1024 * '* 0.00000000000000E+0000 0.00000000000000E+0000 * '* 1.95503421309917E-0004 3.53914399999776E+0001 * '* 3.91006842619834E-0004 5.95684899999760E+0001 * '* 5.86510263929974E-0004 6.54621699999552E+0001 * '* 7.82013685239669E-0004 5.24038399999845E+0001 * '* ... ... * '* 1.98826979472187E-0001 2.77372500000183E-0001 * '* 1.99022482893497E-0001 -2.43361500000174E+0000 * '* 1.99217986314807E-0001 -4.84236799999780E+0000 * '* 1.99413489736116E-0001 -6.02247899999929E+0000 * '* 1.99608993157426E-0001 -5.45615399999951E+0000 * '* 1.99804496578736E-0001 -3.22824200000105E+0000 * '* 2.00000000000045E-0001 -2.96010699999982E-0003 * '* * '* Output file (tgfft.lst): * '* * '* Frequency Value * '* ------------------------------ * '* 0.000000 0.271690 * '* 4.995117 0.546336 * '* 9.990234 0.555560 * '* 14.985352 0.572242 * '* 19.980469 0.598828 * '* 24.975586 0.640061 * '* ... ... * '* 2522.534180 0.020520 * '* 2527.529297 0.020518 * '* 2532.524414 0.020516 * '* 2537.519531 0.020513 * '* 2542.514648 0.020512 * '* 2547.509766 0.020511 * '* 2552.504883 0.020511 * '* * '* A graph of Fourier analysis is displayed. * '*************************************************************** Sub Fourier() Dim I As Integer, n1 As Integer Dim df As Double, dt As Double, f As Double, temp As Double Dim tbegin As Double, tend As Double Form1.AutoRedraw = False Form1.Cls 'open input and output file Open "c:\temp\tgfft.dat" For Input As #1 Open "c:\temp\tgfft.lst" For Output As #2 'read number of input signal points in input file Input #1, ndata 'take nearest power of two If ndata > 2048 Then ndata = 2048: ip = 11 End If If (ndata < 2048) And (ndata > 1023) Then ndata = 1024: ip = 10 End If If (ndata < 1024) And (ndata > 511) Then ndata = 512: ip = 9 End If If (ndata < 512) And (ndata > 255) Then ndata = 256: ip = 8 End If If (ndata < 256) And (ndata > 127) Then ndata = 128: ip = 7 End If If (ndata < 128) And (ndata > 63) Then ndata = 64: ip = 6 End If If (ndata < 64) And (ndata > 31) Then ndata = 32: ip = 5 End If If ndata < 32 Then Print #2, " Error: number of points too small (<32) !" Close #2 Print #2, " Results in file tfft.lst (error)." End End If 'read ndata couples T(i), Y(i) in input data file For I = 1 To ndata Input #1, temp, V(I) If I = 1 Then tbegin = temp If I = ndata Then tend = temp Next I Close #1 'Put input signal in real part of complex vector A For I = 1 To ndata A(I, 1) = V(I) A(I, 2) = 0# Next I 'call FFT procedure FFT 'get frequencies in signal real vector n1 = ndata / 2 For I = 1 To n1 temp = Sqr(A(I, 1) * A(I, 1) + A(I, 2) * A(I, 2)) If (I > 1) Then temp = 2# * temp V(I) = temp If V(I) < vmini Then vmini = V(I) If V(I) > vmaxi Then vmaxi = V(I) Next I 'calculate sampling time range dt dt = (tend - tbegin) / (ndata - 1) 'calculate frequency step df df = 1# / (ndata * dt) 'print frequency specter to output file f = 0# Print #2, " Frequency Value " Print #2, "---------------------------------" For I = 1 To n1 Print #2, f; V(I) f = f + df Next I Close #2 'draw fourier analysis CurveXY 10, n1, 0, n1 * df, 0, 1 Legends 0, 10, "Fourier Analysis", "Freq. hz", "Level m/s2" End Sub Sub FFT() '*************************************************************** '* FAST FOURIER TRANSFORM PROCEDURE * '* ----------------------------------------------------------- * '* This procedure calculates the fast Fourier transform of a * '* real function sampled in N points 0,1,....,N-1. N must be a * '* power of two (2 power M). * '* T being the sampling duration, the maximum frequency in * '* the signal can't be greater than fc = 1/(2T). The resulting * '* specter H(k) is discreet and contains the frequencies: * '* fn = k/(NT) with k = -N/2,.. -1,0,1,2,..,N/2. * '* H(0) corresponds to null fréquency. * '* H(N/2) corresponds to fc frequency. * '* ----------------------------------------------------------- * '* INPUTS: * '* * '* A(i,2) complex vector of size N, the real part of * '* which contains the N sampled points of real * '* signal to analyse (time spacing is constant). * '* * '* M integer such as N=2^M * '* * '* OUTPUTS: * '* * '* A(i,2) complex vector of size N, the vector modulus * '* contain the frequencies of input signal, the * '* vector angles contain the corresponding phases * '* (not used here). * '* * '* J-P Moreau/J-P Dumont * '*************************************************************** Dim U(2) As Double, W(2) As Double, T(2) As Double Dim N As Integer, NV2 As Integer, NM1 As Integer Dim J As Integer, I As Integer, IIP As Integer Dim K As Integer, L As Integer, LE As Integer, LE1 As Integer Dim Phi As Double, Pi As Double Pi = 3.1415926535 N = ndata NV2 = N / 2 NM1 = N - 1 J = 1 For I = 1 To NM1 If (I < J) Then T(1) = A(J, 1): T(2) = A(J, 2) A(J, 1) = A(I, 1): A(J, 2) = A(I, 2) A(I, 1) = T(1): A(I, 2) = T(2) End If K = NV2 100 If K >= J Then GoTo 200 J = J - K K = K / 2 GoTo 100 200 J = J + K Next I LE = 1 For L = 1 To ip LE1 = LE LE = LE * 2 U(1) = 1# U(2) = 0# Phi = Pi / LE1 W(1) = Cos(Phi) W(2) = Sin(Phi) For J = 1 To LE1 I = J - LE 300 If I >= N - LE Then GoTo 400 I = I + LE IIP = I + LE1 T(1) = A(IIP, 1) * U(1) - A(IIP, 2) * U(2) T(2) = A(IIP, 1) * U(2) + A(IIP, 2) * U(1) A(IIP, 1) = A(I, 1) - T(1) A(IIP, 2) = A(I, 2) - T(2) A(I, 1) = A(I, 1) + T(1) A(I, 2) = A(I, 2) + T(2) GoTo 300 400 temp = U(1) U(1) = W(1) * U(1) - W(2) * U(2) U(2) = W(1) * U(2) + W(2) * temp Next J Next L For I = 1 To N A(I, 1) = A(I, 1) / N A(I, 2) = A(I, 2) / N Next I End Sub 'End of file tgfft.bas