' ==================================================== ' FFT_LOGIX4_SMALL_BASIC ' ' Version 5.0 - Date 4 January 2024 ' ' In version 4.0 - Harmonics survive the windowing functions ' 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 - Added Time-Frequency Spectrogram - ' Works with nsamples = 512 ' Works with FFT size = 127 and better if windowing ' function is set to NONE. ' https://smallbasic.com/program/?CFDR662.000 ' 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.logix4u.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 +"========================================="+CRLF Message = Message +"==> WORKS WELL ONLY WITH NSAMPLES = 512 <=="+CRLF Message = Message +"========================================="+CRLF Message = Message +"==> WORKS WELL WITH SPECTROGRAM FFT SIZE = 128 <=="+CRLF Message = Message +"========================================="+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 = " NSAMPLES = 512 - SPECTROGRAM FFT SIZE = 128"+CRLF Message = 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."+CRLF Message = Message +" LAST WARNING!!!!"+CRLF Message = Message +"Some loops may be controlled by hard-coded values. If you"+CRLF Message = Message +"decide to edit any constant or variable, don't forget to"+CRLF Message = Message +"check whether loops exist where a number was hard-coded in"+CRLF Message = Message +"place of the constant / variable name!!" 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 #20 H1_Frequency = 80 ' optional harmonic 1 #80 H2_Frequency = 160 ' optional harmonic 2 #160 sample_rate = 1000 ' sample rate #1000 nsamples = 512 ' samples that will be generated #512 Message = "Non-normalized parameters for fundamental sine wave " Message = Message +"generation are as follows: Sample Rate "+sample_rate+", " Message = Message +" Number of Samples "+nsamples+", Frequency "+Frequency+" Hz, Amplitude "+ampl+"." Message = Message +"The two optional harmonic frequencies are "+H1_Frequency+" and "+H2_Frequency+"." GraphicsWindow.DrawBoundText(525,5,450,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))) * 0.66 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))) * 0.66 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 ' ================================= ' PLOT A SPECTROGRAM FROM THE INITIAL ' WAVEFORM ' Even though it was developed as a last ' complement, and even though it is based on the ' FFT() function, still unexplained up here in the ' source code, I decided to present the ' spectrogram here, before any intervention ' on the initial waveform. ' ================================= Plot_SPECTROGRAM() ' ================================= ' 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 nsamples-1 'For cnt = 1 To nsamples ' ------------------------------------------------------- ' 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 = nsamples - 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 = nsamples - 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) GraphicsWindow.FillRectangle(565,55,265,435) 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 / nsamples)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / nsamples))) 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 / nsamples)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / nsamples))) * 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 / nsamples)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / nsamples))) ' 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 / nsamples)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / nsamples))) / 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 / nsamples)) + (0.08 * Math.Cos(4 * Math.Pi * cnt / nsamples))) / 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 / nsamples))) 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 / nsamples))) * 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 / nsamples))) ' 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 / nsamples))) / 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 / nsamples))) / 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 / nsamples))) 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 / nsamples))) * 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 / nsamples))) ' 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 / nsamples))) / 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 / nsamples))) / 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 nsamples-1 'For smpl = 1 To nsamples ' ------------------------------------------------------------------------------------------------------------------------------------ ' 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 nsamples/2 ' 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 nsamples/2 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 nsamples/2 ' 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 nsamples/2 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 nsamples/2 ' 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 nsamples/2 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 nsamples/2 ' 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 nsamples/2 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 nsamples/2 ' 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 nsamples/2 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(10,515,400,"Vertical 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 nsamples/2 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 ' ============================= ' ============================= ' PLOT_SPECTROGRAM ' ============================= ' ============================= Sub Plot_SPECTROGRAM '======================================= ' WORKS BETTER WHEN WINDOWING FUNCTION ' IS SET TO "NONE" '======================================= ' The samp_buf array will be divided in small discrete chunks ' Each chunk contains Samples_per_column entries ' FFT will be performed on each subset ' The resulting FFT will be used to draw a vertical column in the spectrogram ' [only half the subscripts will be used for reasons clarified elsewhere, ' check FFT() and Prepare_FFT_For_Plotting()] Samples_per_column = nsamples/4 'FFT size - Works Well with 128, 64 halves the column height Sample_From = 0 'Works well with 0 Sample_To = (nsamples/4)-1 'Works well with 127 Iterations = (nsamples / Samples_per_column)-1 '------------------------------------------------------------------- ' Loop for each spectrogram column '------------------------------------------------------------------- For Spectrogram_Columns = 0 To Iterations-1 ' Copy into REX[] the next Samples_per_column samples from samp_buf index_REX = 0 'Always fill the first Samples_per_column columns in REX For ind = Sample_From To Sample_To REX[index_REX] = samp_buf[ind] IMX[index_REX] = 0 index_REX = index_REX + 1 EndFor 'Now the first Samples_per_column entries of REX are filled 'Execute the FFT N = Samples_per_column fft() '------------------------------------------------------------------------------------- ' Prepare FFT for POWER plot as in Sub Prepare_FFT_for_Plotting() ' Do NOT forget that the FFT returns an array that is symmetrically ' divided in two halves, and only one half is used to generate the ' spectrogram and the other plots. ' [only half the subscripts will be used for reasons clarified elsewhere, ' check FFT() and Prepare_FFT_For_Plotting()] '------------------------------------------------------------------------------------- 'TextWindow.WriteLine("Values in the spectrogram array") MIN_SPECT_POWER = 99999 MAX_SPECT_POWER = -99999 For Spectrogram_Rows = 0 To Samples_per_column-1 Spectrogram_Array[Spectrogram_Columns][Spectrogram_Rows] = Math.SquareRoot(((IMX[Spectrogram_Rows]*IMX[Spectrogram_Rows]) + (REX[Spectrogram_Rows]*REX[Spectrogram_Rows]))) If Spectrogram_Array[Spectrogram_Columns][Spectrogram_Rows] < MIN_SPECT_POWER Then MIN_SPECT_POWER = MATH.Round(Spectrogram_Array[Spectrogram_Columns][Spectrogram_Rows]) EndIf If Spectrogram_Array[Spectrogram_Columns][Spectrogram_Rows] > MAX_SPECT_POWER Then MAX_SPECT_POWER = MATH.Round(Spectrogram_Array[Spectrogram_Columns][Spectrogram_Rows]) EndIf ' TextWindow.WriteLine(Spectrogram_Array[Spectrogram_Columns][Spectrogram_Row]) EndFor 'TextWindow.WriteLine("Values in the Spectrogram Array") 'For list = 0 To Samples_per_column-1 ' TextWindow.WriteLine(REX[list]) 'EndFor 'TextWindow.Pause() '------------------------------------------------------------------- ' Next subset (next column) '------------------------------------------------------------------- Sample_From = Sample_From + Samples_per_column Sample_To = Sample_To + Samples_per_column EndFor 'TextWindow.WriteLine("MIN_SPECT_POWER "+ MIN_SPECT_POWER) 'TextWindow.WriteLine("MAX_SPECT_POWER "+ MAX_SPECT_POWER) 'TextWindow.Pause() '------------------------------------------------------------------- ' Brutal color mapping ' Assign a greytone for each possible value ' of power in the spectrogram array ' (lightest value = higher power) '------------------------------------------------------------------- HEXDIGIT[0]="0" HEXDIGIT[1]="1" HEXDIGIT[2]="2" HEXDIGIT[3]="3" HEXDIGIT[4]="4" HEXDIGIT[5]="5" HEXDIGIT[6]="6" HEXDIGIT[7]="7" HEXDIGIT[8]="8" HEXDIGIT[9]="9" HEXDIGIT[10]="A" HEXDIGIT[11]="B" HEXDIGIT[12]="C" HEXDIGIT[13]="D" HEXDIGIT[14]="E" HEXDIGIT[15]="F" COUNT_SLOW = 15 COUNT_FAST= 15 For i = MAX_SPECT_POWER To MIN_SPECT_POWER Step -1 SPECT_COLOR[i] ="#"+HEXDIGIT[COUNT_SLOW]+HEXDIGIT[COUNT_FAST]+HEXDIGIT[COUNT_SLOW]+HEXDIGIT[COUNT_FAST]+HEXDIGIT[COUNT_SLOW]+HEXDIGIT[COUNT_FAST] COUNT_FAST = COUNT_FAST - 1 If COUNT_FAST = -1 Then COUNT_FAST = 15 COUNT_SLOW = COUNT_SLOW - 1 If COUNT_SLOW = -1 Then COUNT_SLOW = 15 EndIf EndIf 'TextWindow.WriteLine(SPECT_COLOR[i]) EndFor 'TextWindow.Pause() GraphicsWindow.BrushColor = "Red" GraphicsWindow.PenColor = "Red" GraphicsWindow.DrawBoundText (670,150,140,"SCALE FOR FFT SIZE = 127 - COLUMN SIZE 64 - CAUTION! THE SPECTROGRAM WORKS BETTER WITH WINDOWING FUNCTION SET TO NONE - Spectrogram Horizontal reference marks are spaced at approximately 20 Hz") GraphicsWindow.DrawRectangle(570,60,250,405) skip=1 Hor_line_step = 16 Ref_line_X = 575 For Ref_line_Y = 455 to 58 Step -Hor_line_step GraphicsWindow.DrawLine(Ref_line_X, Ref_line_Y, Ref_line_X + 20, Ref_line_Y) EndFor For Spectrogram_Columns = 0 To Iterations-1 ' Considering that the FFT algorithm For Spectrogram_Rows = 0 to math.Round((Samples_per_column-1)/2) Spectrogram_Pixel_X = 600 + (Spectrogram_Columns * 20) Spectrogram_Pixel_Y = 450 - (Spectrogram_Rows * 6) ASSIGNED_COLOR = MATH.Round(Spectrogram_Array[Spectrogram_Columns][Spectrogram_Rows]) 'TextWindow.WriteLine("Assegnato colore "+ASSIGNED_COLOR) GraphicsWindow.BrushColor = SPECT_COLOR[ASSIGNED_COLOR] GraphicsWindow.FillRectangle(Spectrogram_Pixel_X,Spectrogram_Pixel_Y,18,5) EndFor EndFor 'TextWindow.Pause() EndSub