Support Forum / Форум поддержки

Code: Select all

' test program for dsp functions in dsp_util
'
graphics ! graphics clear .8,.8,.8
dim x(4096),reX(2049),imX(2049)
N=256 ! M=N/2 ! pi=4*atan(1) ! time reset
signal(N,x,2,50,0,0,"0")        ' generate input signal, 4 components
signal(N,x,4,10,-pi/2.2,0,"+")
signal(N,x,3,13,0,0,"+")
noise(N,x,5,"n","+")
sc=graph("input signal, time domain",N,x,24,0)
r_fft(N,x,reX,imX)                ' perform fourrier transform
graph_magn(N/2,reX,imX,400)     ' output in polar form
button "1" text int(1000*time())&" msec." at 20,650 size 120,40
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 r_fft(N,x(),reX(),imX())
' def c_fft(M,x())
' 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) real 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

' real FFT (using complex FFT)
' passed parameters: see DFT
'
def r_fft(N,x(),reX(),imX())
c_fft(log2(N),x)
for i=0 to N/2 ! reX(i)=real(x(i)) ! imX(i)=imag(x(i)) ! next i
end def

' complex FFT (converted Fortran function)
' M = power of 2 -> N=2^M
' x() is array of type complex
' internally, the function works with option base 1
' input AND results via complex array x()
'
def c_fft(M,x())
pi=4*atan(1) ! N=2^M ! ob=option_base() ! option base 1
for k=1 to M
  ke=2^(M-k+1) ! ke2=ke/2 ! u=1+0i ! angl=pi/ke2
  s=cos(angl)-sin(angl)*1i
  for j=1 to ke2
    for i=j to N step ke
      ip=i+ke2 ! t=x(i)+x(ip)
      x(ip)=u*(x(i)-x(ip)) ! x(i)=t
      next i
    u*=s
    next j
  next k
nd2=N/2 ! nm1=N-1 ! j=1
for i=1 to nm1
  if i<j then ! t=x(j) ! x(j)=x(i) ! x(i)=t ! end if
  k=nd2
  while k<j ! j-=k ! k/=2 ! end while
  j+=k
  next i
option base ob
end def

' inverse discrete fourrier transform
' N = # of samples in resulting time domain x()
' reX(),imX() = input data frequency domain (# = N/2)
'
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 returns 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