'GFXterm - magnitude FFT

'init gfxterm parameters
Const GFX=Chr$(16)
Const ENQ=Chr$(5)
Const ACK=Chr$(6)
Print GFX "?"
Input Gw, Gh
Print GFX "clear" 0, 0, Gw, Gh
GoSub drw

'defines
n=1023   'samples for fft (multiple of 2 - 1)
cy = 10  'cycles of sinewave in n samples

Dim a!(n),fm!(n),fp!(n)
GoSub sinfill
'GoSub sqrfill
window=3      '1=trapezoid, 2=triangle, 3=sin^2
GoSub window  'if required
'GoSub sinlist
'End

GoSub sinshow
Math fft magnitude a!(),fm!()
'GoSub magnitudelog
'Math fft phase a!(),fp!()

GoSub fftshow
'GoSub phshow

End

'window the input signal to remove spurious
window:
 If window=1 Then
   For i=0 To n/cy
     a!(i)=a!(i)*i/(n/cy)
     a!(n-i)=a!(n-i)*i/(n/cy)
   Next i
 ElseIf window =2 Then
   For i=0 To n/2
     a!(i)=a!(i)*i/(n/2)
     a!(n-i)=a!(n-i)*i/(n/2)
   Next i
 ElseIf window=3 Then
   For i=1 To n
     a!(i)=a!(i)*Sin(Pi*i/n)*Sin(Pi*i/n)
   Next i
 EndIf
Return

'show magnitude fft on graphics terminal
fftshow:
 Print GFX "ink" 255, 255, 0, 2
 'scan magnitude array
 mi!=0:ma!=0
 'note: the fft array is symetrical, you only need look at first half
 For i=0 To n/2
   mi!=Min(mi!,fm!(i))
   ma!=Max(ma!,fm!(i))
 Next i
 'Print mi!, Ma!
 ygain=-Gh/(2.2*(ma!-mi!))
 yoffset=Gh-ygain*mi!
 xgain=2*Gw/(n+1)
 xoffset=0
 For i=0 To n/2-1
   x1=i*xgain+xoffset
   x2=(i+1)*xgain+xoffset
   y1=fm!(i)*ygain+yoffset
   y2=fm!(i+1)*ygain+yoffset
   Print GFX "line" x1,y1,x2,y2
 Next i
Return

magnitudelog:
 Print "magintude log"
 For i=0 To n
   fm!(i)=20*Log(fm!(i))/Log(10)
 Next i
Return

'show the phase of the signal
phshow:
 Print GFX "ink" 0, 0, 255, 1
 'scan phase array
 mi!=0:ma!=0
 'note: the fft array is symetrical, you only need look at first half
 For i=0 To n/2
   mi!=Min(mi!,fp!(i))
   ma!=Max(ma!,fp!(i))
 Next i
 'Print mi!,ma!
 ygain=-Gh/(2.2*(ma!-mi!))
 yoffset=3*Gh/4
 xgain=2*Gw/(n+1)
 xoffset=0
 For i=0 To n/2-1
   x1=i*xgain+xoffset
   x2=(i+1)*xgain+xoffset
   y1=fp!(i)*ygain+yoffset
   y2=fp!(i+1)*ygain+yoffset
   Print GFX "line" x1,y1,x2,y2
 Next i
Return

'fill input array a() with sine wave samples
sinfill:
For i=0 To n
 a!(i)=Sin(i*2*Pi*cy/n)
'  a!(i)=a!(i)/2+0.5 'offset
'  Print a!(i),
Next i
Return

'fill input array a() with square wave samples
sqrfill:
period=Int((n+1)/cy)
For i=0 To cy-1
 For j=0 To period/2 - 1
   a!(i*period+j)=1
 Next j
 For j=period/2 To period
   a!(i*period+j)=-1
 Next j
'  Print a!(i),
Next i
Return

'show the graph of the input signal
sinshow:
 Print GFX "ink" 0, 255, 0, 2
 ygain=-Gh/5
 yoffset=Gh/4
 xgain=Gw/(n+1)
 xoffset=0
 For i=0 To n-1
   x1=i*xgain+xoffset
   x2=(i+1)*xgain+xoffset
   y1=a!(i)*ygain+yoffset
   y2=a!(i+1)*ygain+yoffset
   Print GFX "line" x1,y1,x2,y2
 Next i
Return

sinlist:
 For i=1 To n
   Print a!(i)
 Next i
Return

'draw frames
drw:
 Print GFX "ink" 255, 255, 255, 2
 Print GFX "line" 0,Gh,0,Gh/2+1
 Print GFX "line" 0,Gh,Gw,Gh
 Print GFX "line" 0,Gh/2+1,Gw,Gh/2+1
 Print GFX "line" Gw,Gh,Gw,Gh/2+1
 Print GFX "line" 0,0,0,Gh/2-1
 Print GFX "line" 0,0,Gw,0
 Print GFX "line" 0,Gh/2-1,Gw,Gh/2-1
 Print GFX "line" Gw,0,Gw,Gh/2-1
Return

'check fft magnitude function

n=1023 'samples
cy=10  'cycles

dim a!(n),fm!(n)
w=0  'window type, 0=none

do
 cls

 'draw 2 frames using mm.hres mm.vres
 box 50,0,(mm.hres-51),mm.vres/2-2
 box 50,mm.vres/2,mm.hres-51,mm.vres-1

 'prepare input signal
 sinfill
 if k$="w" then w=(w+1) mod 4
 if w then window
 
 'plot input signal
 plot_in

 'do fft
 math fft magnitude a!(),fm!()

 'plot fft  
 plot_fft

 'keyboard commands
 do
   k$=inkey$
 loop until k$<>""
 if k$="a" then cy=cy+1
 if k$="z" then cy=cy-1:cy=max(cy,2)
loop until k$="q"

end

'subroutines

'fill the input array a!() with cy cycles of sinewave
sub sinfill
 for i=0 to n
   a!(i)=sin(2*pi*cy*i/n)
 next i
end sub

'apply a window over the input signal (triangle for now)
sub window
 if w=1 then
   for i=0 to n/2  'linear
     a!(i)=(i/n)*a!(i)
     a!(n-i)=(i/n)*a!(n-i)
   next i
   ? @(60,5) "window linear"
 else if w=2 then
   for i=0 to n    'quadratic (do linear twice)
     a!(i)=(i/n)*a!(i)
     a!(n-i)=(i/n)*a!(n-i)
   next i
   ? @(60,5) "window quadratic"
 else if w=3 then
   for i=0 to n    'cosine window
     a!(i)=(1-cos(2*pi*i/n))*a!(i)
   next i
   ? @(60,5) "window cosine"
 end if
end sub

'plot the input signal in the upper window
sub plot_in
 'find minimum and maximum
 mi!=0:ma!=0
 for i=0 to n
   mi!=min(mi!,a!(i))
   ma!=max(ma!,a!(i))
 next i
 '? mi!,ma!
 'calculate gain and offset to fit it into the window
 xgain=(mm.hres-50)/n
 xoffs=50
 ygain=-(mm.vres/2.3)/(ma!-mi!)
 yoffs=mm.vres/4-ygain*(ma!+mi!)/2
 'plot the actual samples  using linear interpolation
 for i=0 to n-1
   line i*xgain+xoffs,a!(i)*ygain+yoffs,(i+1)*xgain+xoffs,a!(i+1)*ygain+yoffs
 next i
 'add legend
 ? @(mm.hres/2,5) "input linear"
 ? @(5,5) str$(ma!,2,2)
 ? @(5,mm.vres/2-20) str$(mi!,2,2)
end sub

'plot the fft output in the lower window
sub plot_fft
 'find min and max and where max is
 mi!=0:ma!=0
 for i=0 to n/2  'fft output is symetrical
   mi!=min(mi!,fm!(i))
   ma!=max(ma!,fm!(i))
   if ma!=fm!(i) then cf=i':? cf
 next i
 '? mi!,ma!
 'calculate scaling to fit the window
 xgain=2*(mm.hres-50)/n
 xoffs=50
 ygain=-(mm.vres/2.3)/(ma!-mi!)
 yoffs=mm.vres-20
 'plot the fft graph
 for i=0 to n/2-2
   line i*xgain+xoffs,fm!(i)*ygain+yoffs,(i+1)*xgain+xoffs,fm!(i+1)*ygain+yoffs
 next i
 'add legend
 ? @(mm.hres/2,mm.vres/2+15) "magintude linear"
 ? @(5,mm.vres/2+5) str$(ma!,2,2)
 ? @(5,mm.vres-30) str$(mi!,2,2);
 'put marker ticks in multiples of centre frequency
 for i=0 to n/2-2 step cf
   line i*xgain+xoffs,mm.vres/2,i*xgain+xoffs,mm.vres/2+5
   if ((i mod(5*cy)) = 0) then
     line i*xgain+xoffs,mm.vres/2,i*xgain+xoffs,mm.vres/2+15
   end if
 next i
end sub

> option list
OPTION SDCARD GP5
OPTION SYSTEM SPI GP2,GP3,GP4
> files
A:/
     1351  bench.bas
    26240  ccpico.bas
      201  circles.bas
     3164  gameoflife.bas
      618  gettime.bas
      636  grain.bas
    11008  gui.bas
       22  hello.bas
       35  hello0.bas
     1444  mcp23017.bas
      817  menu.bas
      364  pelican.bas
      309  settime.bas
      367  speed.bas
     4515  ta.bas
     4700  tain.txt
     3017  taout.txt
0 directories, 17 files