Almost there. Added a straightforward Fourrier transform and its inverse, and some output stuff.
I will add a FFT function too, although it's speed is not needed at these quantities and fast devices. But it should be there.
The dsp_util file contains all functions and is included from the same directory where the calling test program resides (you can of course modify the include clause).
Code: Select all
' test program for dsp functions in dsp_util
'
graphics ! graphics clear .8,.8,.8
dim x(512),reX(257),imX(257)
N=256 ! M=N/2 ! pi=4*atan(1)
signal(N,x,2,50,0,0,"0") ' generate input signal, 3 components
signal(N,x,4,10,-pi/2.2,0,"+")
noise(N,x,5,"u","+")
sc=graph("input signal, time domain",N,x,24,0)
dft(N,x,reX,imX) ' perform fourrier transform
graph_magn(N/2,reX,imX,400) ' output in polar form
end
{dsp_util}
' def signal (N,v(),ampl,freq,shift,vdc,mode$)
' def noise (N,v(),ampl,distr$,mode$)
' def sigstat(N,v(),st()) (min, max, mean, st_dev)
' def convolute (N,M,x(),h(),y())
' def convol_out(N,M,x(),h(),y())
' def dft(N,x(),reX(),imX())
' def inv_dft(N,reX(),imX(),x())
' def graph(txt$,N,v(),ytop,sc)
' def graph_system(M,h(),ytop)
' def graph_freqs(M,rx(),ix(),ytop)
' def graph_magn(M,rx(),ix(),ytop)
' def box(ytop)
' def box2(ytop)
' def systembox(ytop)
Code: Select all
' sine generator
' N = # of samples
' v() contains (cumulative) generated signal
' ampl = amplitude
' freq = frequency (# of complete cycles)
' shift = phase shift in radians
' vdc = base value for signal
' mode = "+" (add to v()), "0" (init v()), or "-" (subtract from (v))
'
def signal(N,v(),ampl,freq,shift,vdc,mode$)
if mode$<>"+" and mode$<>"0" and mode$<>"-" then return 0
dt=2*freq*.pi/N
for i=0 to N-1
sig=vdc+ampl*sin(i*dt+shift)
if mode$="0" then ! v(i)=sig
else ! if mode$="+" then v(i)+=sig else v(i)-=sig
end if
next i
end def
' random noise generator
' distr$ = Gaussian ("n") or uniform ("u")
' mode$ = add, init, or subtract ("+", "0", or "-")
'
def noise(N,v(),ampl,distr$,mode$)
randomize
if distr$<>"n" and distr$<>"u" then return
if mode$<>"+" and mode$<>"0" and mode$<>"-" then return
for i=0 to N-1
if distr$="u" then ! p=rnd(1)
else ! p=0 ! for j=1 to 12 ! p+=rnd(1) ! next j ! p/=12
end if
p=2*ampl*p-ampl
if mode$="0" then ! v(i)=p
else ! if mode$="+" then v(i)+=p else v(i)-=p
end if
next i
end def
' base statistics of signal
' returns st(4), resp. minimum, maximun, mean, standard deviation
'
def sigstat(N,v(),st())
mini=999999 ! maxi=-999999 ! mean=0 ! st_dev=0
for i=0 to N-1
num=v(i) ! mini=min(mini,num) ! maxi=max(maxi,num) ! mean+=num
next i ! mean/=N
for i=0 to N-1 ! variance+=(v(i)-mean)^2 ! next i
st(0)=mini ! st(1)=maxi ! st(2)=mean ! st(3)=sqrt(variance/(N-1))
return
end def
' Convolution, input side algorithm
' N = # samples in input x()
' M = # samples in impulse response h()
' y() = output
'
def convolute(N,M,x(),h(),y())
for i=0 to N+M-1 ! y(i)=0 ! next i
for i=0 to N-1
for j=0 to M-1 ! y(i+j)+=x(i)*h(j) ! next j
next i
return
end def
' Convolution, output side algorithm
' N = # samples in input x()
' M = # samples in impulse response h()
' y() = output
'
def convol_out(N,M,x(),h(),y())
for i=0 to N+M-1
y(i)=0
for j=0 to M-1
if i>=j and i-j<=N-1 then y(i)+=h(j)*x(i-j)
next j
next i
end def
' simple Fourier transform (gives same results as FFT)
' N = # samples in input x()
' x() input samples (N samples)
' reX() output: N/2+1 amplitudes of cosine components
' imX() output: N/2+1 amplitudes of sine components
'
def dft(N,x(),reX(),imX())
pi=4*atan(1) ! fac=2*pi/N ! M=N/2
for k=0 to M
f=fac*k ! reX(k)=0 ! imX(k)=0
for i=0 to N-1
angl=f*i
reX(k)+=x(i)*cos(angl) ! imX(k)-=x(i)*sin(angl)
next i
next k
end def
' inverse discrete fourrier transform
' N = # of samples in resulting time domain
' reX(),imX() = input data frequency domain
' resulting time domain
'
def inv_dft(N,reX(),imX(),x())
pi=4*atan(1) ! fac=2*pi/N ! M=N/2
for i=0 to M ! reX(i)/=M ! imX(i)/=-M ! next i
reX(0)/=2 ! reX(M)/=2
for i=0 to N-1 ! x(i)=0 ! next i
for k=0 to M
f=fac*k
for i=0 to N-1
angl=f*i
x(i)+=reX(k)*cos(angl) ! x(i)+=imX(k)*sin(angl)
next i
next k
end def
' graphs a series of N samples in v()
' ytop = top of graph-box
' sc = scale, if scale=0, maximum scale will be calculated
' function return the scale used.
'
def graph(txt$,N,v(),ytop,sc)
dim stat(4)
sigstat(N,v,stat) ! maxi=stat(1) ! mini=stat(0)
box(ytop) ! yc=150+ytop ! draw size 1 ! draw color 1,0,0
if sc=0 then
if maxi>0 then sc1=150/maxi else sc1=9999
if mini<0 then sc2=-150/mini else sc2=9999
sc=min(sc1,sc2)
end if
dx=720/(N-1)
for i=0 to N-1
if i=0 then draw to 24,yc-sc*v(0) else draw line to 24+i*dx,yc-sc*v(i)
next i
draw font size 12 ! draw color 0,0,1
draw text txt$ at 24,ytop-16
draw font size 20
return sc
end def
' graphs an impuls responce of M samples in h()
' ytop = top of graph-box
'
def graph_system(M,h(),ytop)
systembox(ytop)
sc=50 ! yc=ytop+100 ! dx=200/(M-1)
fill color 1,0,0
for i=0 to M-1
yy=yc-sc*h(i)
if abs(h(i))<=2 then fill circle 284+i*dx,yy size 3
next i
draw font size 12
draw text "system's impulse response" at 284,ytop-16
draw font size 20
end def
' graphs the transform result in cartesian form
' M = # of cosine and sine components = length of frequency axe.
' rx() contains the cosine components, ix() the sine components
' ytop = top of graph-box
'
def graph_freqs(M,rx(),ix(),ytop)
dim rxabs(M+1),riabs(M+1),stat(4)
for i=0 to M ! rx(i)=abs(rx(i)) ! ix(i)=abs(ix(i)) ! next i
sigstat(M,rx,stat) ! max1=stat(1)
sigstat(M,ix,stat) ! max2=stat(1)
maxi=max(max1,max2) ! sc=200/maxi
box2(ytop) ! ybot=200+ytop ! draw size 1 ! draw color 1,0,0
dx=360/(M-1) ! x1=16 ! x2=393
for i=0 to M-1
if rx(i) then draw line x1+i*dx,ybot to x1+i*dx,ybot-sc*(rx(i))
if ix(i) then draw line x2+i*dx,ybot to x2+i*dx,ybot-sc*(ix(i))
next i
draw font size 12 ! draw color 0,0,1 ! draw alpha 1
draw text "frequency domain, cosine functions" at x1,ytop-15
draw text "frequency domain, sine functions" at x2,ytop-15
draw font size 20 ! draw color 0,0,0
end def
' graphs the transform result in polar form
' M = # of cosine components in rx() and sine components in ix()
' output shows magnitudes and phase shifts
' ytop = top of graph-box
'
def graph_magn(M,rx(),ix(),ytop)
dim mag(M+1),phase(M+1),stat(4)
pi=4*atan(1) ! eps=0.00001
for i=0 to M
re=rx(i) ! im=ix(i)
mag(i)=sqrt(re*re+im*im)
if re=0 then phase(i)=pi/2 else phase(i)=atan(im/re)
if re<0 then ! phase(i)+=pi
else ! if im<0 then phase(i)+=2*pi
end if
next i
sigstat(M,mag,stat) ! sc1=200/stat(1)
sigstat(M,phase,stat) ! sc2=200/stat(1)
box2(ytop) ! ybot=200+ytop ! draw size 1 ! draw color 1,0,0
dx=360/(M-1) ! x1=16 ! x2=393
for i=0 to M-1
if mag(i) then draw line x1+i*dx,ybot to x1+i*dx,ybot-sc1*(mag(i))
if abs(ix(i))>eps then draw line x2+i*dx,ybot to x2+i*dx,ybot-sc2*(phase(i))
next i
draw font size 12 ! draw color 0,0,1 ! draw alpha 1
draw text "frequency domain, magnitudes" at x1,ytop-15
draw text "frequency domain, phase shifts" at x2,ytop-15
draw font size 20 ! draw color 0,0,0
end def
def box(ytop)
refresh off
get screen size sw,sh ! xc=384 ! yc=150+ytop ! ybot=150+yc
draw size 3 ! draw color 0,0,0
draw rect 24,ytop to 744,ybot
draw size 1 ! draw alpha .3
for x=24 to 744 step 20 ! draw line x,ytop to x,ybot ! next x
for y=ytop to ybot step 15 ! draw line 24,y to 744,y ! next y
draw size 2 ! draw color 0,0,1 ! draw alpha 1
draw line 24,yc to 744,yc
refresh on
end def
def box2(ytop)
refresh off
x1=16 ! x2=393 ! ybot=ytop+200
draw size 3 ! draw color 0,0,0
draw rect x1,ytop to x1+360,ybot
draw rect x2,ytop to x2+360,ybot
draw size 1 ! draw alpha .3
for x=x1 to x1+360 step 18
draw line x,ytop to x,ybot
draw line x+376,ytop to x+376,ybot
next x
for y=ytop to ybot step 20
draw line x1,y to x1+360,y
draw line x2,y to x2+360,y
next y
refresh on
end def
def systembox(ytop)
refresh off
xc=384 ! yc=100+ytop ! ybot=100+yc
draw size 3 ! draw color 0,0,0
draw rect 284,ytop to 484,ybot
draw size 1 ! draw alpha .3
for x=284 to 484 step 20 ! draw line x,ytop to x,ybot ! next x
for y=ytop to ybot step 10 ! draw line 284,y to 484,y ! next y
draw font size 12
draw text 2 at xc,ytop ! draw text 1 at xc,yc-55
draw text -1 at xc,yc+45 ! draw text -2 at xc,yc+88
draw font size 20
draw size 2 ! draw color 0,0,1 ! draw alpha 1
draw line 284,yc to 484,yc
refresh on
end def
DSP lib in progress or already done?
