MMBasic for Windows - betas

Author Message
Volhout
Guru

Joined: 05/03/2018
Location: Netherlands
Posts: 3340
 Posted: 02:28pm 01 Mar 2022

Hi Peter,

Please check the MATH FFT function. For CMM2 we had to fix it, and for MMB4W we seem to run into a similar issue (although not identical).

When I take the FFT magnitude from a sinesquare windowed sine wave, and display it in linear scale, it should look like this (CMM2 fft calculation from below program):

However, the MMB4W math function shows many spurious peaks. Below graph shows the sinesquare windowed sine signal (top), as well as the fft (lower graph)

I use this code to check it (the program can change the windowing function with "w", and change the number of sinewaves for fft with "a" and "z", and save the graph to disk with "s").

'check fft magnitude function

n=1023 'samples
cy=10  'cycles
sq=0   'sq=1 = square wave. sq=0 = sine wave

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!()

'convert to logaritm
'lin2log

'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)
if k\$="s" then save image "fft",0,0,800,600
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)
if sq=1 then
if a!(i)>0 then
a!(i)=1
else
a!(i)=-1
end if
end if
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) "triangle window"
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
else if w=3 then
for i=0 to n    'sine^2 (hann) window
a!(i)=sin(pi*i/n)*sin(pi*i/n)*a!(i)
next i
? @(60,5) "hann window"
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
? @(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!=1e6:ma!=-1e6
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
'apply some logaritmic scaling
'mi!=-120:ma!=60
'mi!=0:ma!=200
'mi!=max(mi!,-140):ma!=min(ma!,120)
'for i=0 to n
'  fm!(i)=min(ma!,max(mi!,fm!(i)))
'next i
'calculate scaling to fit the window
xgain=2*(mm.hres-50)/n
xoffs=50
ygain=-(mm.vres/2.3)/(ma!-mi!)
yoffs=mm.vres-mi!*ygain-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
? @(mm.hres/2,mm.vres/2+15) "magintude logarithm"
? @(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

'convert linear output to logarithmic
sub lin2log
for i=0 to n/2
fm!(i)= 20*log(fm!(i)+1e-300)/log(10)
next i
end sub

Edited 2022-03-02 17:27 by Volhout