' ******************************************************************
' Fast Fourier Transform (modified from Pascal program by Don Cross)
' ******************************************************************
'
' The program generates a time signal consisting of a large 200 Hz
' sine wave added to a small 2000 Hz cosine wave, which is graphed
' on the screen (Press a key after you are done viewing each graph)
'
' Next, it performs the FFT and graphs the resulting complex
' frequency samples.
'
' Then, it filters out all frequency components above 1000 Hz in
' the transformed data.
'
' Finally, it performs the inverse transform to get a filtered
' time signal back, and graphs the result.
'
' In addition, the program compares the FFT with the direct
' computation, on a set of random data.
'
' Results are stored in the output file fftout.txt
' ******************************************************************

#INCLUDE "fbmath.bi"
#INCLUDE "fbgfx.bi"

CONST NumSamples   =  512                       ' Buffer size must be power of 2
CONST SamplingRate = 22050                      ' Sampling rate in Hz
CONST MaxIndex     = NumSamples - 1             ' Max. array index
CONST MidIndex     = NumSamples \ 2
CONST DT           = 1 / SamplingRate           ' Time unit
CONST DF           = SamplingRate / NumSamples  ' Frequency unit

FUNCTION F(T AS DOUBLE) AS DOUBLE
' Original signal
  F = SIN(200 * TwoPi * T) + 0.2 * COS(2000 * TwoPi * T)
END FUNCTION

SUB Test_CalcFrequency(InArray() AS Complex, OutArray() AS Complex)

  DIM AS Complex Z
  DIM AS INTEGER I

  ' Fill input buffers with random data
  FOR I = 0 TO MaxIndex
    InArray(I).X = 10000 * RND
    InArray(I).Y = 10000 * RND
  NEXT I

  PRINT #1, ""
  PRINT #1, " *** Testing procedure CalcFrequency *** "
  PRINT #1, ""

  FFT InArray(), OutArray()

  FOR I = 0 TO MaxIndex
    Z = CalcFrequency(I, InArray())
    PRINT #1, USING "##########"; I;
    PRINT #1, USING "##########.######"; OutArray(I).X; Z.X; OutArray(I).Y; Z.Y
  NEXT I
END SUB

SUB ListData(DataArray() AS Complex, Comment AS STRING)

  DIM AS INTEGER I

  PRINT #1, " *** "; Comment; " *** "
  PRINT #1, ""
  PRINT #1, "     index             real             imag"

  FOR I = 0 TO MaxIndex
    PRINT #1, USING "##########"; I;
    PRINT #1, USING "##########.######"; DataArray(I).X; DataArray(I).Y
  NEXT I

  PRINT #1, ""
  PRINT #1, STRING$(44, "-")
  PRINT #1, ""
END SUB

SUB PlotData(T() AS DOUBLE, Z() AS Complex, Title AS STRING)

  DIM AS DOUBLE  Xmin, Xmax, Xstep  ' Ox scale
  DIM AS DOUBLE  Ymin, Ymax, Ystep  ' Oy scale
  DIM AS INTEGER I                  ' Loop variable

  DIM AS DOUBLE X(0 TO MaxIndex)    ' Real part of Z

  IF NOT InitScreen(12, 15, 0) THEN
    PRINT "Unable to set graphic mode!"
    EXIT SUB
  END IF

  InitGraph

  FOR I = 0 TO MaxIndex : X(I) = Z(I).X : NEXT I

  AutoScale T(), LinScale, Xmin, Xmax, Xstep
  AutoScale X(), LinScale, Ymin, Ymax, Ystep

  SetOxScale LinScale, Xmin, Xmax, Xstep
  SetOyScale LinScale, Ymin, Ymax, Ystep

  PlotOxAxis Ymin, 5
  PlotOyAxis Xmin, 2

  PlotGrid

  WriteTitles Title, "Time (s)", "Amplitude"

  COLOR 10 : PlotCurve T(), X(), 0, 0, -1

  SLEEP
  SCREEN 0
END SUB

SUB PlotFFT(Freq() AS DOUBLE, Z() AS Complex, Title AS STRING)

  DIM AS DOUBLE Fq(0 TO MidIndex)  ' Frequency
  DIM AS DOUBLE X(0 TO MidIndex)   ' Real part of FFT
  DIM AS DOUBLE Y(0 TO MidIndex)   ' Imag. part of FFT

  ' Graphic parameters for real and imaginary parts of FFT
  DIM AS TCurvParam CurvParam(1 TO 2) = _
    {(12, 0, 0, -1, "Real"),            _
     (14, 0, 0, -1, "Imag")}

  DIM AS DOUBLE  Xmin, Xmax, Xstep        ' Ox scale
  DIM AS DOUBLE  Yr_min, Yr_max, Yr_step  ' Oy scale (real part)
  DIM AS DOUBLE  Yi_min, Yi_max, Yi_step  ' Oy scale (imag. part)
  DIM AS DOUBLE  Ymin, Ymax, Ystep        ' Oy scale (global)

  DIM AS INTEGER I                        ' Loop variable

  IF NOT InitScreen(12, 15, 0) THEN
    PRINT "Unable to set graphic mode!"
    EXIT SUB
  END IF

  InitGraph

  FOR I = 0 TO MidIndex
    Fq(I) = Freq(I)
    X(I) = Z(I).X
    Y(I) = Z(I).Y
  NEXT I

  AutoScale Fq(), LinScale, Xmin, Xmax, Xstep
  AutoScale X(), LinScale, Yr_min, Yr_max, Yr_step
  AutoScale Y(), LinScale, Yi_min, Yi_max, Yi_step

  Ymin  = Min(Yr_min, Yi_min)
  Ymax  = Max(Yr_max, Yi_max)
  Ystep = Min(Yr_step, Yi_step)

  SetOxScale LinScale, Xmin, Xmax, Xstep
  SetOyScale LinScale, Ymin, Ymax, Ystep

  PlotOxAxis Ymin, 5
  PlotOyAxis Xmin, 2

  PlotGrid

  WriteTitles Title, "Frequency (Hz)", "FFT"

  WITH CurvParam(1)
    COLOR .Col
    PlotCurve Fq(), X(), .Symbol, .Size, .Connect
  END WITH

  WITH CurvParam(2)
    COLOR .Col
    PlotCurve Fq(), Y(), .Symbol, .Size, .Connect
  END WITH

  WriteLegend CurvParam()

  SLEEP
  SCREEN 0
END SUB

DIM AS DOUBLE  T(0 TO MaxIndex)         ' Time
DIM AS DOUBLE  Freq(0 TO MaxIndex)      ' Frequencies
DIM AS Complex InArray(0 TO MaxIndex)   ' FFT input
DIM AS Complex OutArray(0 TO MaxIndex)  ' FFT output

DIM AS INTEGER I, FreqIndex, SymIndex

FOR I = 0 TO MaxIndex
  T(I) = I * DT
  Freq(I) = I * DF
  InArray(I).X = F(T(I))
  InArray(I).Y = 0
NEXT I

OPEN "fftout.txt" FOR OUTPUT AS #1

ListData InArray(), "Time domain data before transform"
PlotData T(), InArray(), "Original signal"

FFT InArray(), OutArray()
PlotFFT Freq(), OutArray(), "Fourier Transform"

ListData OutArray(), "Frequency domain data after transform"

' Filter out everything above 1000 Hz (low-pass)
FreqIndex = 1000 / DF
SymIndex = NumSamples - FreqIndex

FOR I = 0 TO MaxIndex
  IF (I > FreqIndex AND I < MidIndex) OR _
     (I >= MidIndex AND I < SymIndex) THEN
       OutArray(I).X = 0
       OutArray(I).Y = 0
  END IF
NEXT I

IFFT OutArray(), InArray()

ListData InArray(), "Time domain data after inverse transform"
PlotData T(), InArray(), "Filtered signal"

Test_CalcFrequency InArray(), OutArray()

CLOSE #1

