'11/18/2006
'This block of code performs a discrete Fourier transform (DFT).  You may use it freely
'but I would appreciate if you would leave this title block in place.  This and more free
'examples (in MathCad and Excel) can be found at my website.  Updates to this particular
‘code block can be found at www.planetsource.code
'
'Russell W. Patterson, P.E.
'rwpatterson357@aol.com
'http://webpages.charter.net/rwpatterson357/
'
'Reference: "The Scientist and Engineer's Guide to Digital Signal Processing"
'by Steven W. Smith  ISBN: 0-9660176-3-3, Website: www.DSPguide.com '

'Below I load up 16 values of sampled data in the "Sample" array.  These are from the
'original sampled waveform.  Since one complete cycle of fundamental is 16 samples in
'this case, N = 16.  This number is variable of course depending on how many samples are
'taken across the fundamental frequency.
'Sample(0) = 4.943
'Sample(1) = 4.215
'Sample(2) = 1.655
'Sample(3) = -1.144
'Sample(4) = -2.546
'Sample(5) = -4.467
'Sample(6) = -6.195
'Sample(7) = -5.658
'Sample(8) = -4.743
'Sample(9) = -4.015
'Sample(10) = -1.455
'Sample(11) = 1.344
'Sample(12) = 2.746
'Sample(13) = 4.667
'Sample(14) = 6.395
'Sample(15) = 5.858

Dim Real(16) As Single  'You can make this double if need.  This array holds the real
                        'value results of the DFT
Dim Imag(16) As Single  'Likewise, this array holds the imaginary results
Dim Mag(16) As Single   'This array holds the magnitude of the DFT results
Dim Phase(16) As Single 'This array holds the phase angles of the DFT results
Dim pi As Single
pi = 4 * Atn(1)         'Computes pi (3.14)
 'In this case the sample frequency is 960Hz and I am interested in a 60Hz fundamental.
 'N works out to 16 samples per cycle, Nyquist tells us that we can get information about
 'the harmonics up to N/2, and that if any higher frequencies are present that aliasing
 'will occur and our results will be suspect (which is why low-pass filtering may be required).

N = 960 / 60 'samples per 60Hz cycle

'calulate the real and imaginary components of each harmonic (perform DFT)
'This is the only code needed to perform the DFT.
For F = 0 To (N / 2)  'F=0 is the zero frequency (DC) component, F=1 is the fundamental
 For i = 0 To (N - 1) 'step through each of the 16 samples in 1 fundamental cycle
  Real(F) = Real(F) + Sample(i) * Cos(2 * pi * i * F / N)
  Imag(F) = Imag(F) + Sample(i) * Sin(2 * pi * i * F / N)
 Next i
Imag(F) = -Imag(F)
Next F
Real(0) = Real(0) / 2  'correct the magnitudes of the end bin components
Real(N / 2) = Real(N / 2) / 2

'This remaining code is aimed at putting the DFT rectangular complex (r + jx) results
'in polar form.

'compute harmonic magnitudes from the real and imaginary results
For F = 0 To (N / 2)
 Mag(F) = (Real(F) ^ 2 + Imag(F) ^ 2) ^ 0.5 * (2 / N)  'the 2/N is required to correct mag
Next F

'this iteration removes very small values (not required depending on the amplitudes of
'interest.  in this case it zeros harmonics that are less than 1000th of the fundamental)
'The only reason to do this is to improve ability to display results cleaner
For F = 0 To (N / 2)
 If (Mag(F) / Mag(1)) < 10 ^ -3 Then Mag(F) = 0  'get rid of bug dust
Next F

'calculate the phase angle of each harmonic (in radians, multiply by 180/pi for degrees)
For F = 0 To (N / 2)
'get rid of zeros to allow proper phase angle calculation, this is required so the Atn
'function won't blow up
 If Abs(Real(F)) < 10 ^ -10 Then Real(F) = 10 ^ -10
 Phase(F) = Atn(Imag(F) / Real(F))  'phase angles in radians
 If Abs(Phase(F)) < 10 ^ -10 Then Phase(F) = 0  'zero out small phase angles, not required
 If Mag(F) < 10 ^ -10 Then Phase(F) = 0 'zero out phase of zero amplitude harmonics
Next F

 'correct phase angle of zero frequency component (DC)
 If Sgn(Real(0)) < 0 Then
   Phase(0) = -pi
 Else
   Phase(0) = pi
 End If
 
 'correct phase angle of non zero frequency components
For i = 1 To (N / 2)
 If Sgn(Real(i)) > 0 Then
   'leave phase as-is
 Else
   If Sgn(Imag(i)) > 0 Then
     Phase(i) = Phase(i) + pi
   Else
     Phase(i) = Phase(i) - pi
   End If
 End If
Next i
‘end of code block