' ====================================================
' FFT_LOGIX4_SMALL_BASIC
'
' Version 4.0  - Date 5 January 2023
'
' In version 2.0 - Added INVERSE FFT with plot of the resulting waveform
' In version 3.0 -  Added BAND STOP FILTERS
' In version 3.1 -  revised the range of FOR loops, not from 1 to 512 but from 0 to 511
'
' This version - harmonics survive the windowing functions
'
' https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
' ====================================================
' THE FAST FOURIER TRANSFORM
' copyright © 1997-1999 by California Technical Publishing
' published with  permission from Steven W Smith, www.dspguide.com
' modified by logix4u (www.logix4.net)
' modified and ported in Small Basic by Cesare Brizio
' (www.cesarebrizio.it)
' ---------------------------------------------------------
' COOLEY - TOKEY ALGORITHM
' ---------------------------------------------------------
' Bit reversal is most important for radix-2 Cooley–Tukey FFT algorithms, 
' where the recursive stages of the algorithm, operating in-place, imply a bit 
' reversal of the inputs or outputs. Similarly, mixed-radix digit reversals arise 
' in mixed-radix Cooley–Tukey FFTs.
' ====================================================
 
'REX[ ] holds the real part of the frequency domain
'IMX[ ] holds the imaginary part of the frequency domain
 
GraphicsWindow.Top = 10
GraphicsWindow.Left = 10
GraphicsWindow.Show()
GraphicsWindow.Height=1000
GraphicsWindow.Width=1000
GraphicsWindow.BackgroundColor = "White"
GraphicsWindow.BrushColor = "Black"
 
CRLF = Text.GetCharacter(13) + Text.GetCharacter(10)
Message =          "This program is a rudimentary, non-parametric, fully hard-coded,"+CRLF
Message = Message +"very partial Small Basic porting of the Logix4 version of an FFT"+CRLF
Message = Message +"demonstrator, one of the very few available with full, explicit"+CRLF
Message = Message +"and crystal-clear source code."+CRLF
Message = Message +"A few options can be activated by selectively enabling or"+CRLF
Message = Message +"disabling commented lines in the source code include three"+CRLF 
Message = Message +"optional windowing functions (Hamming, Hanning, Blackman) and"+CRLF
Message = Message +"four alternatives on which FFT values should be plotted "+CRLF
Message = Message +"(Real, Imaginary, Magnitude / Power, Phase - the last option"+CRLF
Message = Message +"is the more messy and less convincing)."+CRLF
Message = Message +"I do admit that I'm a shoddy programmer, and that I stopped"+CRLF
Message = Message +"as soon as the program output approximatively matched the"+CRLF
Message = Message +"results of the original. I have no excuses."+CRLF
Message = Message +"Yet, I encountered some well-noted limitations of Small Basic"+CRLF
Message = Message +"that contributed to the poor final result."+CRLF
Message = Message +"    ====> BY NO MEANS THIS SHOULD BE CONSIDERED AN  <====    "+CRLF
Message = Message +"    ====> EXHAUSTIVE EXPLANATION OR IMPLEMENTATION  <====    "+CRLF
Message = Message +"    ====> OF THE ORIGINAL ALGORITHM!ANYBODY NEEDING <====    "+CRLF
Message = Message +"    ====> ADVICE ABOUT CODING THE FFT MAY REFER TO: <====    "+CRLF
Message = Message +"https://mathonweb.com/resources/book2/FFT-Of-A-Wavepacket.pdf"+CRLF
Message = Message +"     https://rosettacode.org/wiki/Fast_Fourier_transform     "+CRLF
Message = Message +"                                                             "+CRLF
Message = Message +"This program is a starting point that may surely be improved "+CRLF
Message = Message +"and made more interactive with limited effort."
 
GraphicsWindow.ShowMessage(Message,"«FFT LOGIX4 SMALL BASIC» - Version 3.0")
 
Message =          "For the most effective rendering you should pay attention"+CRLF
Message = Message +"to both the text window asking your input, and to the"+CRLF
Message = Message +"source code where specific options can be activated"+CRLF
Message = Message +"or inactivated by placing pr removing comment marks."+CRLF
Message = Message +"1) Choose the desired windowing functions by suitably"+CRLF
Message = Message +"   placing or removing comment marks in the source code"+CRLF
Message = Message +"   around the line number 185. If in doubt, activate NONE."+CRLF
Message = Message +"2) Choose the MAGNITUDE FFT display by suitably placing"+CRLF
Message = Message +"   or removing comment marks in the source code of the"+CRLF
Message = Message +"   Prepare_FFT_for_Plotting subroutine"+CRLF
Message = Message +"3) Launch the program"+CRLF
Message = Message +"4) When asked, opt for including harmonic frequencies"+CRLF
Message = Message +"5) If you opt for the inclusion of harmonic frequencies,"+CRLF
Message = Message +"   you'll be asked about the BAND STOP FILTER: you should"+CRLF
Message = Message +"   apply it - it makes more sense with"+CRLF
Message = Message +"   the MAGNITUDE FFT display."+CRLF
Message = Message +"6) When asked, you can choose whether or not"+CRLF
Message = Message +"   you like to perform an INVERSE FOURIER TRANSFORM"+CRLF
Message = Message +"   and re-create the waveform from the data in the"+CRLF
Message = Message +"   frequency domain."
 
GraphicsWindow.ShowMessage(Message,"MOST ADVISABLE SETTINGS")
 
' ====================================================
' Generate a sine wave 512 samples long at 20 Hz
' Prepare optional harmonic frequencies to be added
' ====================================================
 
ampl = 5 ' see comments more under - maximum amplitude of the base sine wave
Frequency = 20 ' natural fundamental frequency of the sine wave
H1_Frequency = 80 ' optional harmonic 1
H2_Frequency = 160 ' optional harmonic 2
sample_rate = 1000 
nsamples = 512 ' samples that will be generated
 
Message =  "Non-normalized parameters for fundamental sine wave "
Message = Message +"generation are as follows: Sample Rate 1000, "
Message = Message +"  Number of Samples 512, Frequency 10, Amplitude 5."
GraphicsWindow.DrawBoundText(600,20,350,Message)
 
TextWindow.Top=100
TextWindow.Left=1300
ASK_STYLE_AGAIN:
TextWindow.WriteLine("Would you like to generate a simple fundamental sine wave (F) or a fundamental with two harmonic frequencies (H) ?")
STYLE = Text.ConvertToUpperCase(TextWindow.Read())
If STYLE<>"F" And STYLE<>"H" Then
  TextWindow.WriteLine("You should enter «F» or «H»!")
  Goto ASK_STYLE_AGAIN
EndIf
 
' ----------------------------------------------------------------------------
' The bare fundamental may have an amplitude up to 
' the ampl value set above
' I decide that the max amplitude engages 50 pixel on screen
' ----------------------------------------------------------------------------
maximum_pixels = 50
 
If STYLE = "F" Then
  ' ----------------------------------------------------------------------------
  ' Vertical scale of a fundamental-only view
  ' ----------------------------------------------------------------------------  
  pixel_per_unit = math.Round(maximum_pixels / ampl)    
  For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
    samp_buf[cnt] = ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))
    ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
    ' It's apparently impossible to multiply an array element with other variables or constants, 
    ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
    ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
    Y_pos[cnt] = maximum_pixels + 10 - ((ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * pixel_per_unit)
  EndFor
Else
  ' ----------------------------------------------------------------------------
  ' By adding more sine waves, I will exceed the maximum
  ' ampl value set above
  ' I will scale the pixel_per_unit to reduce the vertical size
  ' of the  wave
  ' ----------------------------------------------------------------------------
  pixel_per_unit = math.Round((maximum_pixels / ampl) / 2)
 
  For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
    fund_buf[cnt] = ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))
  EndFor
  For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
    H1_buf[cnt] = (ampl * Math.Sin(2 * Math.Pi * H1_Frequency * (cnt / sample_rate))) / 2
  EndFor
  For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
    H2_buf[cnt] = (ampl * Math.Sin(2 * Math.Pi * H2_Frequency * (cnt / sample_rate))) / 4
  EndFor
  'Combination
  For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
    samp_buf[cnt] = fund_buf[cnt] + H1_buf[cnt] + H2_buf[cnt]
    ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
    ' It's apparently impossible to multiply an array element with other variables or constants, 
    ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
    ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
    Y_pos[cnt] = maximum_pixels + 10 - ((fund_buf[cnt] + H1_buf[cnt] + H2_buf[cnt]) * pixel_per_unit)
  EndFor
 
EndIf
 
' ------------------------------------------
' PLOT THE GENERATED WAVE
' ------------------------------------------
Plot_Input_Wave()
 
' ====================================================
'                                           WINDOWING FUNCTIONS
' -----------------------------------------------------------------------------------------
' A WINDOWING function may be chosen, particularly useful in case the 
' samples in the signal do not amount to entire wave cycles.
' In fact, the actual FFT transform assumes that the signal provided is a finite
'  data set, a continuous spectrum that is one period of a periodic signal.
' Windowing functions pre-process the sample as described in 
' "Understanding FFTs and windowing functions" by National Instruments
' https://download.ni.com/evaluation/pxi/Understanding%20FFTs%20and%20Windowing.pdf
' ====================================================
 
Win_Funct = "NONE"
'Win_Funct = "HAMM"
'Win_Funct = "HANN"
'Win_Funct = "BLKM"
 
 
If Win_Funct = "HAMM" Then
  ApplyHamming()
  GraphicsWindow.Clear()
  Plot_input_Wave()
EndIf
 
If Win_Funct = "HANN" Then
  ApplyHanning()
  GraphicsWindow.Clear()
  Plot_input_Wave()
EndIf
 
If Win_Funct = "BLKM" Then
  ApplyBlackman()
  GraphicsWindow.Clear()
  Plot_input_Wave()
EndIf
 
' =================================
' FILL ARRAYS FOR THE  FFT SUBROUTINE
' =================================
' Reminder: Small basic des not support the
' COMPLEX data type.
' ===============================
' Store the data in two separate arrays,
' one for the real part and one
' for the imaginary part
' ===============================
 
For cnt = 0 To 511 'For cnt = 1 To 512
  ' -------------------------------------------------------
  ' The real part is taken from samp_buf
  ' that may have been modified by
  ' the application of a window
  ' -------------------------------------------------------
  REX[cnt] = samp_buf[cnt]
  ' -------------------------------------------------------
  ' imaginary part will be calculated by the FFT
  ' it's not used for time-based renderings, but
  ' it allows any frequency-based rendering
  ' -------------------------------------------------------
  IMX[cnt] = 0
EndFor
 
' ========================
' CALL FFT SUBROUTINE
' ========================
N = nsamples
 
TextWindow.Top=100
TextWindow.Left=1300
TextWindow.WriteLine("FFT in progress - wait patiently!")
 
' ---------------------------------------------------------------
' DEBUG - Save pre-FFT values on disk for check
' ---------------------------------------------------------------
'FilePath = "C:\FFT Python 3.10\RealValBeforeFFT.txt"'
'For zz = 1 to nsamples
'File.InsertLine(FilePath,zz,zz+") "+ REX[zz])
'EndFor
'FilePath = "C:\FFT Python 3.10\ImagValBeforeFFT.txt"
'For zz = 1 to nsamples
'File.InsertLine(FilePath,zz,zz+") "+ IMX[zz])
'EndFor
 
FFT()
 
TextWindow.WriteLine("FFT completed!")
 
' ---------------------------------------------------------------
' DEBUG - Save post-FFT values on disk for check
' ---------------------------------------------------------------
'FilePath = "C:\FFT Python 3.10\RealValAfterFFT.txt"'
'For zz = 1 to nsamples
'  File.InsertLine(FilePath,zz,zz+") "+ REX[zz])
'EndFor
'FilePath = "C:\FFT Python 3.10\ImagValAfterFFT.txt"
'For zz = 1 to nsamples
'  File.InsertLine(FilePath,zz,zz+") "+ IMX[zz])
'EndFor
 
TextWindow.WriteLine("Output array preparation in progress - wait patiently!")
 
' ================================
' PREPARE FFT  FOR PLOTTING
' -------------------------------------------------------
' Several values can be obtained from the FFT,
' including:
' REAL PART
' IMAGINARY PART
' MAGNITUDE
' PHASE (this is the part that I've understood 
' more poorly: I cannot explain the plot
' resulting when the PHASE option is
' selected in Sub Prepare_FFT_for_Plotting() )
' ================================
Prepare_FFT_for_Plotting()
 
TextWindow.WriteLine("Output array preparation complete! Displaying FFT Results.")
 
' ================================
' PLOT FFT
' ================================
Plot_FFT()
 
GraphicsWindow.ShowMessage("FFT plotted!","FFT plotted!")
 
' =============================================================
' Ask whether BANDSTOP should be executed on the POWER SPECTRUM FFT
' =============================================================
If STYLE = "H" Then
 
  If OPTION <> "MAGNITUDE" Then
    Message = "            R E M E M B E R   !!!  !!!"+CRLF  
    Message = "The effect of band stopping will be evident in the"+CRLF
    Message = Message +"MAGNITUDE view of the FFT. If another view is selected,"+CRLF
    Message = Message +"the result may not be apparent."+CRLF
    GraphicsWindow.ShowMessage(Message,"Friendly reminder")
  EndIf
 
  ASK_BANDSTOP_AGAIN:
  TextWindow.WriteLine("Please understand that the band filtering metod used (zeroing")
  TextWindow.WriteLine("the IMX[] and REX[] values around the filtered frequencies")  
  TextWindow.WriteLine("is brutal and noise-inducing. This option is provided")  
  TextWindow.WriteLine("only as a way to introduce the concept of how an")    
  TextWindow.WriteLine("intervention on the frequency domain may impact")      
  TextWindow.WriteLine("the original waveshape in the time domain.")      
  TextWindow.WriteLine("Following the suppression, the FFT will be plotted again.")
  TextWindow.WriteLine(" ")  
  TextWindow.WriteLine("Would you like to execute a BAND STOP FILTER to suppress")
  TextWindow.WriteLine("the two harmonic frequencies?")
  OK_BANDSTOP = Text.ConvertToUpperCase(TextWindow.Read())
  If OK_BANDSTOP<>"Y" And OK_BANDSTOP<>"N" Then
    TextWindow.WriteLine("You should enter «Y» or «N»!")
    Goto ASK_BANDSTOP_AGAIN
  EndIf
 
  If OK_BANDSTOP = "Y" Then
     ' =============================================================
     ' I know that the two harmonic frequencies to suppress amount to H1_Frequency = 80 
     ' and H2_Frequency = 160.
     ' Considering that the two vectors  IMX[] and REX[] filled by the FFT have as many
     ' elements and that there is a negative-positive symmetry as side effect of putting 
     ' real-valued input into the Fourier transform, it can safely be concluded that, to affect  
     ' the two harmonic frequencies, I should zero the positions around the indexes [40] and [80]
     ' e.g. from 37 to 43 and from 77 to 83. 
     ' Symmetrically, I should zero the two frequencies at (512- 40) and (512-80),
     ' otherwise the harmonics will re-emerge after the inverse FFT 
     ' =============================================================
     First_Freq_To_Stop = Math.round(H1_Frequency/2)
     For Stop_it = First_Freq_To_Stop-3 to First_Freq_To_Stop+3
       IMX[Stop_it] = 0
       REX[Stop_it] = 0
     EndFor
     Second_Freq_To_Stop = Math.round(H2_Frequency/2)
     For Stop_it = Second_Freq_To_Stop-3 to Second_Freq_To_Stop+3
       IMX[Stop_it] = 0
       REX[Stop_it] = 0
     EndFor
     First_Symm_Freq_To_Stop = 512 - Math.round(H1_Frequency/2)
     For Stop_it = First_Symm_Freq_To_Stop-3 to First_Symm_Freq_To_Stop+3
       IMX[Stop_it] = 0
       REX[Stop_it] = 0
     EndFor
     Second_Symm_Freq_To_Stop = 512 - Math.round(H2_Frequency/2)
     For Stop_it = Second_Symm_Freq_To_Stop-3 to Second_Symm_Freq_To_Stop+3
       IMX[Stop_it] = 0
       REX[Stop_it] = 0
     EndFor     
 
 
     ' ==================================================
     ' Let's plot the power spectrum again
     ' ==================================================
     GraphicsWindow.BrushColor="White"
     GraphicsWindow.FillRectangle(0,100,990,990)
 
     TextWindow.WriteLine("Output array preparation in progress - wait patiently!")
 
     Prepare_FFT_for_Plotting()
 
     TextWindow.WriteLine("Output array preparation complete! Displaying FFT Results.")
 
     Plot_FFT()
 
     GraphicsWindow.ShowMessage("FFT plotted!","FFT plotted after harmonic frequencies suppression!")
  EndIf
EndIf
 
' =============================================================
' Ask whether INVERSE FFT should be executed
' =============================================================
ASK_INVERSE_AGAIN:
TextWindow.WriteLine("Would you like to execute an INVERSE FFT and plot the resulting wave above the original?")
TextWindow.WriteLine("(if you answer N the program will terminate)")
OK_INVERSE = Text.ConvertToUpperCase(TextWindow.Read())
If OK_INVERSE<>"Y" And OK_INVERSE<>"N" Then
  TextWindow.WriteLine("You should enter «Y» or «N»!")
  Goto ASK_INVERSE_AGAIN
EndIf
 
If OK_INVERSE = "Y" Then
 
  ' ---------------------------------------------------------------
  ' Execute INVERSE fast fourier transform
  ' ---------------------------------------------------------------  
 
  TextWindow.WriteLine("Inverse FFT in progress - wait patiently!")
 
  ifft()
 
  TextWindow.WriteLine("Inverse FFT complete! Displaying Inverse FFT Results (re-created original wave).")
 
  ' ---------------------------------------------------------------
  ' DEBUG - Save post-IFFT values on disk for check
  ' ---------------------------------------------------------------
  'FilePath = "C:\FFT Python 3.10\RealValAfterIFFT.txt"'
  'For zz = 1 to nsamples
  'File.InsertLine(FilePath,zz,zz+") "+ REX[zz])
  'EndFor
  'FilePath = "C:\FFT Python 3.10\ImagValAfterIFFT.txt"
  'For zz = 1 to nsamples
  'File.InsertLine(FilePath,zz,zz+") "+ IMX[zz])
  'EndFor
  'TextWindow.Read()
 
  GraphicsWindow.PenColor="Blue"
  Plot_Input_Wave()
 
  GraphicsWindow.ShowMessage("Re-created original wave plotted! Click OK to terminate the program","INVERSE FFT plotted!")
 
EndIf
 
 
Program.End()
 
' ===================================================
' ===================================================
' ===================================================
'                                     S  U  B  R  O  U  T  I  N  E  S
' ----------------------------------------------------------------------------------------
' The program can be effectively split into a few more separate subroutines. 
' I was too lazy in that respect.
' ===================================================
' ===================================================
' ===================================================
 
' ===================================================
' ===================================================
' APPLY BLACKMAN WINDOW TO SAMP_BUF[]
' ===================================================
' ===================================================
Sub ApplyBlackman
  ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
  ' For comments about pixel_per_unit see the fundamental wave generation phase in the main program
  If STYLE = "F" Then
    ' ----------------------------------------------------------------------------
    ' Vertical scale of a fundamental-only view
    ' ----------------------------------------------------------------------------  
    pixel_per_unit = math.Round(maximum_pixels / ampl)    
    For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
      samp_buf[cnt] = (ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.42 - (0.5 * Math.Cos(2 * 3.14159265 * cnt / 512)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / 512)))
      Y_pos[cnt] = maximum_pixels + 10 - ((ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.42 - (0.5 * Math.Cos(2 * 3.14159265 * cnt / 512)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / 512))) * pixel_per_unit)
    EndFor
  Else
    ' ----------------------------------------------------------------------------
    ' Vertical scale of a fundamental plus harmonics view
    ' ----------------------------------------------------------------------------
    pixel_per_unit = math.Round((maximum_pixels / ampl) / 2)
 
    For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
      component_F = (ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.42 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / 512)))
      ' first harmonic is assigned half the amplitude of the fundamental
      component_H1 = (ampl * Math.Sin(2 * Math.Pi * H1_Frequency * (cnt / sample_rate))) * (0.42 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / 512))) / 2
      ' second harmonic is assigned one quarter the amplitude of the fundamental      
      component_H2 = (ampl * Math.Sin(2 * Math.Pi * H2_Frequency * (cnt / sample_rate))) * (0.42 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / 512))) / 4
      'Combination
      samp_buf[cnt] = component_F + component_H1 + component_H2
      ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
      ' It's apparently impossible to multiply an array element with other variables or constants, 
      ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
      ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
      Y_pos[cnt] = maximum_pixels + 10 - ((component_F + component_H1 + component_H2) * pixel_per_unit)
    EndFor
  EndIf
EndSub
 
 
' ===================================================
' ===================================================
' APPLY HAMMING WINDOW TO SAMP_BUF[]
' ===================================================
' ===================================================
Sub ApplyHamming
  ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
  ' For comments about pixel_per_unit see the fundamental wave generation phase in the main program
  If STYLE = "F" Then
    ' ----------------------------------------------------------------------------
    ' Vertical scale of a fundamental-only view
    ' ----------------------------------------------------------------------------  
    pixel_per_unit = math.Round(maximum_pixels / ampl)    
    For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
      samp_buf[cnt] = (ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.54 - (0.46 * Math.Cos(2 * Math.Pi * cnt / 512)))
      Y_pos[cnt] = maximum_pixels + 10 - ((ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.54 - (0.46 * Math.Cos(2 * Math.Pi * cnt / 512))) * pixel_per_unit)
    EndFor
  Else
    ' ----------------------------------------------------------------------------
    ' Vertical scale of a fundamental plus harmonics view
    ' ----------------------------------------------------------------------------
    pixel_per_unit = math.Round((maximum_pixels / ampl) / 2)
 
    For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
      component_F = (ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.54 - (0.46 * Math.Cos(2 * Math.Pi * cnt / 512)))
      ' first harmonic is assigned half the amplitude of the fundamental
      component_H1 = (ampl * Math.Sin(2 * Math.Pi * H1_Frequency * (cnt / sample_rate))) * (0.54 - (0.46 * Math.Cos(2 * Math.Pi * cnt / 512))) / 2
      ' second harmonic is assigned one quarter the amplitude of the fundamental      
      component_H2 = (ampl * Math.Sin(2 * Math.Pi * H2_Frequency * (cnt / sample_rate))) * (0.54 - (0.46 * Math.Cos(2 * Math.Pi * cnt / 512))) / 4
      'Combination
      samp_buf[cnt] = component_F + component_H1 + component_H2
      ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
      ' It's apparently impossible to multiply an array element with other variables or constants, 
      ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
      ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
      Y_pos[cnt] = maximum_pixels + 10 - ((component_F + component_H1 + component_H2) * pixel_per_unit)
    EndFor
  EndIf  
EndSub
 
 
' ===================================================
' ===================================================
' APPLY HANNING WINDOW TO SAMP_BUF[]
' ===================================================
' ===================================================
Sub ApplyHanning
  ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
  ' For comments about pixel_per_unit see the fundamental wave generation phase in the main program
  If STYLE = "F" Then
    ' ----------------------------------------------------------------------------
    ' Vertical scale of a fundamental-only view
    ' ----------------------------------------------------------------------------  
    pixel_per_unit = math.Round(maximum_pixels / ampl)    
    For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
      samp_buf[cnt] = (ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.5 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512)))
      Y_pos[cnt] = maximum_pixels + 10 - ((ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.5 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512))) * pixel_per_unit)
    EndFor
  Else
    ' ----------------------------------------------------------------------------
    ' Vertical scale of a fundamental plus harmonics view
    ' ----------------------------------------------------------------------------
    pixel_per_unit = math.Round((maximum_pixels / ampl) / 2)
 
    For cnt = 0 To nsamples-1 'For cnt = 1 To nsamples
      component_F = (ampl * Math.Sin(2 * Math.Pi * Frequency * (cnt / sample_rate))) * (0.5 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512)))
      ' first harmonic is assigned half the amplitude of the fundamental
      component_H1 = (ampl * Math.Sin(2 * Math.Pi * H1_Frequency * (cnt / sample_rate))) * (0.5 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512))) / 2
      ' second harmonic is assigned one quarter the amplitude of the fundamental      
      component_H2 = (ampl * Math.Sin(2 * Math.Pi * H2_Frequency * (cnt / sample_rate))) * (0.5 - (0.5 * Math.Cos(2 * Math.Pi * cnt / 512))) / 4
      'Combination
      samp_buf[cnt] = component_F + component_H1 + component_H2
      ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
      ' It's apparently impossible to multiply an array element with other variables or constants, 
      ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
      ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
      Y_pos[cnt] = maximum_pixels + 10 - ((component_F + component_H1 + component_H2) * pixel_per_unit)
    EndFor
  EndIf  
EndSub
 
 
' ===================================================
' ===================================================
' PERFORM FAST FOURIER TRANSFORM
' =================================================== 
' ===================================================
Sub fft
' ===================================================
' Upon entry, 
'    N contains the number of points in the DFT, 
'    REX[ ] and
'    IMX[ ] contain the real and imaginary parts of the input. 
' Upon return,
'    REX[ ] and 
'    IMX[ ] contain the DFT output. 
' All signals run from 0 to N-1.
' ===================================================
NM1 = N - 1
ND2 = N / 2
M = Math.Round(Math.Log(N) / Math.Log(2))
J = ND2
' -------------------------
' Bit reversal sorting
' Needed to optimize
' execution time and
' memory usage.
' -------------------------
For i = 1 To N - 2 
  If J > i Then
    TR = REX[J]
    TI = IMX[J]
    REX[J] = REX[i]
    IMX[J] = IMX[i]
    REX[i] = TR
    IMX[i] = TI
  EndIf
  K = ND2
I_LINE_1200:
  If K > J Then 
    GoTo I_LINE_1240
  EndIf  
  J = J - K
  K = K / 2
  GoTo I_LINE_1200
I_LINE_1240:
  J = J + K
EndFor
 
' -------------------------------------------------------------------------------------------------    
' LOOP FOR EACH STAGE
' See https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
' -------------------------------------------------------------------------------------------------    
For L = 1 To M         
  LE = Math.Power(2,L) 
  LE2 = LE / 2         
  UR = 1  ' initial real multiplier is set at 1
  UI = 0  ' initial imaginary multiplier is set at 0  
  ' ---------------------------------------
  ' Calculate sine & cosine values
  ' ---------------------------------------
  SR = Math.Cos(Math.Pi / LE2) 
  SI = -1 * Math.Sin(Math.Pi / LE2)
  ' ---------------------------------------
  ' Loop for each sub DFT
  ' ---------------------------------------
  For J = 1 To LE2 
    JM1 = J - 1
    ' -------------------------------------------------------------------
    ' Loop for each butterfly
    ' See https://en.wikipedia.org/wiki/Butterfly_diagram
    ' -------------------------------------------------------------------
    For i = JM1 To NM1 Step LE 
      IP = i + LE2
      ' ---------------------------------------
      ' Butterfly calculation
      ' ---------------------------------------
      TR = REX[IP] * UR - IMX[IP] * UI
      TI = REX[IP] * UI + IMX[IP] * UR
      REX[IP] = REX[i] - TR
      IMX[IP] = IMX[i] - TI
      REX[i] = REX[i] + TR
      IMX[i] = IMX[i] + TI
    EndFor
    TR = UR 'save UR in TR! Its current value will be needed in a moment
    UR = TR * SR - UI * SI 'recalculate UR
    UI = TR * SI + UI * SR 'recalculate UI
  EndFor
EndFor
 
EndSub
 
 
' ===================================================
' ===================================================
' PERFORM INVERSE FAST FOURIER TRANSFORM
' =================================================== 
' =================================================== 
' The only differences with sub fft() are:
' DIRECT FFT    ===>  SR = Math.Cos(Math.Pi / LE2)
' DIRECT FFT    ===>  SI = -1 * Math.Sin(Math.Pi / LE2)
' INVERSE FFT ===>  SR = Math.Cos(Math.Pi / (LE2*-1))
' INVERSE FFT ===>  SI = -1 * Math.Sin(Math.Pi / (LE2*-1))
' INVERSE FFT ===>  A final loop multiplicates all the values in REX[]
'                                         and IMX[] by 1/nsamples
'  ---------------------------------------------------------------------------------------
' ...BUT OBSERVE THAT AFTER THE INVERSE FOURIER TRANSFORM
'    WE WILL BE INTERESTED IN THE REAL PART ONLY!!!!
' ===================================================
Sub ifft
' ===================================================
' Upon entry, 
'    N contains the number of points in the DFT, 
'    REX[ ] and
'    IMX[ ] contain the real and imaginary parts of the input. 
' Upon return,
'    REX[ ] and 
'    IMX[ ] contain the DFT output. 
' All signals run from 0 to N-1.
' ===================================================
NM1 = N - 1
ND2 = N / 2
M = Math.Round(Math.Log(N) / Math.Log(2))
J = ND2
' -------------------------
' Bit reversal sorting
' Needed to optimize
' execution time and
' memory usage.
' -------------------------
For i = 1 To N - 2 
  If J > i Then
    TR = REX[J]
    TI = IMX[J]
    REX[J] = REX[i]
    IMX[J] = IMX[i]
    REX[i] = TR
    IMX[i] = TI
  EndIf
  K = ND2
LINE_1200:
  If K > J Then 
    GoTo LINE_1240
  EndIf  
  J = J - K
  K = K / 2
  GoTo LINE_1200
LINE_1240:
  J = J + K
EndFor
 
' ----------------------------------    
' LOOP FOR EACH STAGE
' ----------------------------------    
For L = 1 To M
  LE = Math.Power(2,L) ' 2 at the Lth power
  LE2 = LE / 2         '  2 at the Lth power /2
  UR = 1  ' initial real component is set at 1
  UI = 0  ' initial imaginary component is set at 0
  ' ---------------------------------------
  ' Calculate sine & cosine values
  ' the -1 does the trick (compare 
  ' with fft() above)
  ' ---------------------------------------
  SR = Math.Cos(Math.Pi / (LE2*-1))
  SI = -1 * Math.Sin(Math.Pi / (LE2*-1))
 
  ' ---------------------------------------
  ' Loop for each sub DFT
  ' ---------------------------------------
  For J = 1 To LE2 
    JM1 = J - 1
    ' -------------------------------------------------------------------
    ' Loop for each butterfly
    ' See https://en.wikipedia.org/wiki/Butterfly_diagram
    ' -------------------------------------------------------------------
    For i = JM1 To NM1 Step LE 
      IP = i + LE2   
      ' ---------------------------------------
      ' Butterfly calculation
      ' ---------------------------------------
      TR = REX[IP] * UR - IMX[IP] * UI
      TI = REX[IP] * UI + IMX[IP] * UR
      REX[IP] = REX[i] - TR
      IMX[IP] = IMX[i] - TI
      REX[i] = REX[i] + TR
      IMX[i] = IMX[i] + TI
    EndFor
    TR = UR                
    UR = TR * SR - UI * SI  
    UI = TR * SI + UI * SR
  EndFor
EndFor
 
' ------------------------------------------------------------------
' I should multiply the value by 1 / N
' I do this both for the real part and the imaginary part
' but after the inverse FFT I will be interested  just 
' in the real part
' ------------------------------------------------------------------
For i = 1 To N
  REX[i] = REX[i] * (1/N)
  IMX[i] = IMX[i] * (1/N)
  ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
  ' It's apparently impossible to multiply an array element with other variables or constants, 
  ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
  ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
  Curr_REX = REX[i]
  Y_pos[I] = maximum_pixels + 10 - (Curr_REX * pixel_per_unit)  
EndFor
 
EndSub
 
 
' ===================================================
' ===================================================
' PLOT INPUT WAVE
' ===================================================
' ===================================================
Sub Plot_Input_Wave
 
  GraphicsWindow.PenWidth=1
  GraphicsWindow.DrawLine(0, maximum_pixels+10, 520,maximum_pixels+10)
 
  last_Y = maximum_pixels+10
  last_X = 0
  For smpl = 0 To 511 'For smpl = 1 To 512
    ' ------------------------------------------------------------------------------------------------------------------------------------
    ' The Pos_Y array
    ' -----------------------------------------------------------------------------------------------------------------------------------
    ' Apparently, it is not possible to multiply successfully values from an array and current memory variables
    ' such as in     
    ' plot_Y = 240 - (samp_buf[smpl] * pixel_per_unit)
    ' or in 
    ' curr_val = samp_buf[smpl] 
    ' plot_Y = 240 - (curr_val * pixel_per_unit)
    ' because such a multiplication returns ZERO.
    ' As a consequence, the Y coordinate is stored in the array Y_pos during the wave generation process.
    ' -----------------------------------------------------------------------------------------------------------------------------------   
    plot_Y = Y_pos[smpl]
 
    ' ----------------------------------------------------------------
    ' Write a dot (dismissed: the look is much better
    ' when drawing lines between two successive points)
    ' ----------------------------------------------------------------
    ' GraphicsWindow.SetPixel(smpl, plot_Y,"Red")
    ' TextWindow.WriteLine(" Pixel per unit " + pixel_per_unit + " - current samp_buf value " + curr_val + " - X plot coordinate = "+smpl+" - Y plot coordinate = "+plot_Y)
    ' TextWindow.Read()
    ' -------------------------------------------------------------
    ' Write a line
    ' -------------------------------------------------------------
    GraphicsWindow.DrawLine(last_X, last_Y, smpl, plot_Y)
    last_X = smpl
    last_Y = plot_Y    
 
  EndFor
 
EndSub
 
 
' ===================================================
' ===================================================
' PREPARE FFT FOR PLOTTING
' ===================================================
' ===================================================
Sub Prepare_FFT_for_Plotting
 
' ========================
' PREPARE OUTPUT ARRAY TO
' DISPLAY FFT RESULTS
' ========================
 
'OPTION = "REAL"
'OPTION = "IMAGINARY"
OPTION = "MAGNITUDE"
'OPTION = "POWER"
'OPTION="PHASE"
 
' ====================
' Output REAL PART
' ===================
If OPTION = "REAL" Then
  GraphicsWindow.BrushColor = "Red"
 
  ' ===================================
  ' SCALING FOR VIDEO PLOTTING
  ' ===================================
  'Empirically, it was observed that - if compared to the Visual basic original - some output options result 
  'in a much smaller vertical scale. It's necessary, at the price of doubling the time of execution, to know in advance 
  'the biggest and the smallest values encountered, so that the Y range could be suitably scaled o fit on screen
  'Lowest_Value=0
  'Highest_Value=0
  'For cnt = 1 To 256
  '  Curr_Val = REX[cnt] 
  '  If  Curr_Val > Highest_Value Then
  '    Highest_Value = Curr_Val
  '  EndIf  
  '  If  Curr_Val < Lowest_Value Then
  '    Lowest_Value = Curr_Val
  '  EndIf     
  'EndFor
  'Value_Range = Highest_Value-Lowest_Value
  'TextWindow.WriteLine("Highest_Value "+ Highest_Value + " - Lowest_Value  " + Lowest_Value + " - Value_Range " + Value_Range)
  'TextWindow.Read()     
  '
  ' Based on the code above, with around 500 pixel above and 500 pixel under the reference line, a scaling factor of 50% is applied 
  ' to the REAL option 
  ' -----------------------------------------------------------------------------------------------------------------------------------------------------------
 
  For cnt = 1 To 256
    outputarray[cnt] = REX[cnt]
    ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
    ' It's apparently impossible to multiply an array element with other variables or constants, 
    ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
    ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
    Curr_REX = REX[cnt]
    Y_FFT[cnt] =  Curr_REX / 2
 
    'TextWindow.WriteLine("CNT "+ cnt + " - Y_FFT  " + Y_FFT[cnt] + " - REX[cnt] " + REX[cnt] + " - Curr_REX " + Curr_REX)
    'TextWindow.Read()
 
  EndFor
EndIf
 
 
' ====================
' Output IMAGINARY PART
' ===================
If OPTION = "IMAGINARY" Then
  GraphicsWindow.BrushColor = "Green"
 
  ' ===================================
  ' SCALING FOR VIDEO PLOTTING
  ' ===================================
  'Empirically, it was observed that - if compared to the Visual basic original - some output options result 
  'in a much smaller vertical scale. It's necessary, at the price of doubling the time of execution, to know in advance 
  'the biggest and the smallest values encountered, so that the Y range could be suitably scaled o fit on screen
  'Lowest_Value=0
  'Highest_Value=0
  'For cnt = 1 To 256
  '  Curr_Val = IMX[cnt] 
  '  If  Curr_Val > Highest_Value Then
  '    Highest_Value = Curr_Val
  '  EndIf  
  '  If  Curr_Val < Lowest_Value Then
  '    Lowest_Value = Curr_Val
  '  EndIf     
  'EndFor
  'Value_Range = Highest_Value-Lowest_Value
  'TextWindow.WriteLine("Highest_Value "+ Highest_Value + " - Lowest_Value  " + Lowest_Value + " - Value_Range " + Value_Range)
  'TextWindow.Read()     
  '
  ' Based on the code above, with around 500 pixel above and 500 pixel under the reference line, a scaling factor
  ' of 0.5 is applied to the IMAGINARY option
  ' -----------------------------------------------------------------------------------------------------------------------------------------------------------
 
  For cnt = 1 To 256
    outputarray[cnt] = IMX[cnt]
    ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
    ' It's apparently impossible to multiply an array element with other variables or constants, 
    ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
    ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
    Curr_IMX = IMX[cnt]
    ' value scaled to half, see above
    Y_FFT[cnt] =  Curr_IMX / 2
 
    'TextWindow.WriteLine("CNT "+ cnt + " - Y_FFT  " + Y_FFT[cnt] + " - IMX[cnt] " + IMX[cnt] + " - Curr_IMX " + Curr_IMX)
    'TextWindow.Read()
 
  EndFor
EndIf
 
 
' ====================
' Output MAGNITUDE
' apparently is equal to
' POWER VALUE
' ===================
If OPTION = "MAGNITUDE" Then
  GraphicsWindow.BrushColor = "Blue"
 
  ' ===================================
  ' SCALING FOR VIDEO PLOTTING
  ' ===================================
  'Empirically, it was observed that - if compared to the Visual basic original - some output options result 
  'in a much smaller vertical scale. It's necessary, at the price of doubling the time of execution, to know in advance 
  'the biggest and the smallest values encountered, so that the Y range could be suitably scaled o fit on screen
  'Lowest_Value=0
  'Highest_Value=0
  'For cnt = 1 To 256
  '  Curr_Val = Math.SquareRoot(((IMX[cnt]*IMX[cnt]) + (REX[cnt]*REX[cnt]))) 
  '  If  Curr_Val > Highest_Value Then
  '    Highest_Value = Curr_Val
  '  EndIf  
  '  If  Curr_Val < Lowest_Value Then
  '    Lowest_Value = Curr_Val
  '  EndIf     
  'EndFor
  'Value_Range = Highest_Value-Lowest_Value
  'TextWindow.WriteLine("Highest_Value "+ Highest_Value + " - Lowest_Value  " + Lowest_Value + " - Value_Range " + Value_Range)
  'TextWindow.Read()     
  '
  ' Based on the code above, with around 500 pixel above and 500 pixel under the reference line, a scaling factor
  ' of 0.33 is applied to the MAGNITUDE option
  ' -----------------------------------------------------------------------------------------------------------------------------------------------------------
 
  For cnt = 1 To 256
    outputarray[cnt] = Math.SquareRoot(((IMX[cnt]*IMX[cnt]) + (REX[cnt]*REX[cnt])))
    ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
    ' It's apparently impossible to multiply an array element with other variables or constants, 
    ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
    ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
    Curr_IMX = IMX[cnt]
    Curr_REX = REX[cnt]    
    ' value scaled to half, see above
    Y_FFT[cnt] = Math.SquareRoot(((Curr_IMX*Curr_IMX) + (Curr_REX*Curr_REX))) / 3
 
    'TextWindow.WriteLine("CNT "+ cnt + " - Y_FFT  " + Y_FFT[cnt] + " - REX[cnt] " + REX[cnt] + " - Curr_REX " + Curr_REX + " - IMX[cnt] " + IMX[cnt] + " - Curr_IMX " + Curr_IMX)
    'TextWindow.Read()
 
   EndFor
EndIf
 
 
 
' ====================
' Output POWER VALUE
' apparently is equal to
' MAGNITUDE
' ====================
If OPTION = "POWER" Then
  GraphicsWindow.BrushColor = "Magenta"
 
  ' ===================================
  ' SCALING FOR VIDEO PLOTTING
  ' ===================================  
  'Empirically, it was observed that - if compared to the Visual basic original - some output options result 
  'in a much smaller vertical scale. It's necessary, at the price of doubling the time of execution, to know in advance 
  'the biggest and the smallest values encountered, so that the Y range could be suitably scaled o fit on screen
  'Lowest_Value=0
  'Highest_Value=0
  'For cnt = 1 To 256
  '  Curr_Val = Math.SquareRoot(((IMX[cnt]*IMX[cnt]) + (REX[cnt]*REX[cnt]))) 
  '  If  Curr_Val > Highest_Value Then
  '    Highest_Value = Curr_Val
  '  EndIf  
  '  If  Curr_Val < Lowest_Value Then
  '    Lowest_Value = Curr_Val
  '  EndIf     
  'EndFor
  'Value_Range = Highest_Value-Lowest_Value
  'TextWindow.WriteLine("Highest_Value "+ Highest_Value + " - Lowest_Value  " + Lowest_Value + " - Value_Range " + Value_Range)
  'TextWindow.Read()     
  '
  ' Based on the code above, with around 500 pixel above and 500 pixel under the reference line, a scaling factor
  ' of 0.33 is applied to the POWER VALUE option
  ' -----------------------------------------------------------------------------------------------------------------------------------------------------------  
 
  For cnt = 1 To 256
    outputarray[cnt] = Math.SquareRoot(((IMX[cnt]*IMX[cnt]) + (REX[cnt]*REX[cnt])))
    ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
    ' It's apparently impossible to multiply an array element with other variables or constants, 
    ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
    ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
    Curr_IMX = IMX[cnt]
    Curr_REX = REX[cnt]     
    ' value scaled to half, see above   
    Y_FFT[cnt] = Math.SquareRoot(((Curr_IMX*Curr_IMX) + (Curr_REX*Curr_REX))) / 3
 
    'TextWindow.WriteLine("CNT "+ cnt + " - Y_FFT  " + Y_FFT[cnt] + " - REX[cnt] " + REX[cnt] + " - Curr_REX " + Curr_REX + " - IMX[cnt] " + IMX[cnt] + " - Curr_IMX " + Curr_IMX)
    'TextWindow.Read()
 
  EndFor
EndIf
 
 
'====================
' Output PHASE
' ====================
If OPTION = "PHASE" Then
  ' ===================================
  ' SCALING FOR VIDEO PLOTTING
  ' ===================================    
  'Empirically, it was observed that - if compared to the Visual basic original - some output options result 
  'in a much smaller vertical scale. It's necessary, at the price of doubling the time of execution, to know in advance 
  'the biggest and the smallest values encountered, so that the Y range could be suitably scaled o fit on screen
  GraphicsWindow.BrushColor = "DarkRed"
 
  'Empirically, it was observed that - if compared to the Visual basic original - some output options result 
  'in a much smaller vertical scale. It's necessary, at the price of doubling the time of execution, to know in advance 
  'the biggest and the smallest values encountered, so that the Y range could be suitably scaled o fit on screen
  'Lowest_Value=0
  'Highest_Value=0
  'For cnt = 1 To 256
  '  Curr_Val = Math.ArcTan(IMX[cnt]/REX[cnt]) 
  '  If  Curr_Val > Highest_Value Then
  '    Highest_Value = Curr_Val
  '  EndIf  
  '  If  Curr_Val < Lowest_Value Then
  '    Lowest_Value = Curr_Val
  '  EndIf     
  'EndFor
  'Value_Range = Highest_Value-Lowest_Value
  'TextWindow.WriteLine("Highest_Value "+ Highest_Value + " - Lowest_Value  " + Lowest_Value + " - Value_Range " + Value_Range)
  'TextWindow.Read()     
  '
  ' Based on the code above, with around 500 pixel above and 500 pixel under the reference line, a scaling factor
  ' of 100 is applied to the PHASE option
  ' -----------------------------------------------------------------------------------------------------------------------------------------------------------  
 
  For cnt = 1 To 256
    outputarray[cnt] = Math.ArcTan(IMX[cnt]/REX[cnt])
    Curr_IMX = IMX[cnt]
    Curr_REX = REX[cnt]    
    ' The existence of the Y_Pos array is explained under the Plot_Input_Wave Subroutine
    ' It's apparently impossible to multiply an array element with other variables or constants, 
    ' as the multiplication returns ZERO, so I could not use samp_buf[] to set the Y offset when
    ' plotting the curve on screen. I save the Y offset as soon as I calculate the value of samp_buf[cnt]
    Y_FFT[cnt] = Math.ArcTan(Curr_IMX/Curr_REX) * 100
 
    'TextWindow.WriteLine("CNT "+ cnt + " - Y_FFT  " + Y_FFT[cnt] + " - REX[cnt] " + REX[cnt] + " - Curr_REX " + Curr_REX + " - IMX[cnt] " + IMX[cnt] + " - Curr_IMX " + Curr_IMX)
    'TextWindow.Read()
 
  EndFor
EndIf
EndSub 
 
 
' =============================
' =============================
' PLOT_FFT - DISPLAY FFT RESULTS
' =============================
' =============================
Sub Plot_FFT
  GraphicsWindow.DrawBoundText(800,350,100,"Horizontal reference marks are spaced at approximately 10 Hz")
 
 
  GraphicsWindow.PenColor="Black"
  GraphicsWindow.DrawLine(0, 500, 1000,500)
 
  ' reference lines are placed at around 10 Hz
  For grid = 20 to 980 step 20
    GraphicsWindow.DrawLine(grid,490,grid,510)
  EndFor
 
  Baseline_FFT=500
 
  last_Y = Baseline_FFT
  last_X = 0
  For cnt = 1 To 256
    Curr_FFT_Value = Y_FFT[cnt]
    ' -------------------------------------------------------------
    ' Write a dot
    ' -------------------------------------------------------------
    'GraphicsWindow.SetPixel(cnt, Baseline_FFT-Curr_FFT_Value,"Green")
    ' -------------------------------------------------------------
    ' Write a line
    ' -------------------------------------------------------------
    GraphicsWindow.PenColor="Green"
    GraphicsWindow.DrawLine(last_X, last_Y, cnt*4, Baseline_FFT-Curr_FFT_Value)
    last_X = cnt * 4
    last_Y = Baseline_FFT-Curr_FFT_Value  
 
  EndFor
 
EndSub