From gos.ukc.ac.uk!harrier.ukc.ac.uk!ukc!mcvax!uunet!portal!cup.portal.com!dbell Tue Jan 10 18:41:42 GMT 1989
Article 1073 of rec.music.synth:
Path: gos.ukc.ac.uk!harrier.ukc.ac.uk!ukc!mcvax!uunet!portal!cup.portal.com!dbell
>From: dbell@cup.portal.com (David J Bell)
Newsgroups: rec.music.synth
Subject: Re: Re: Computer based FM editing
Message-ID: <13331@cup.portal.com>
Date: 8 Jan 89 21:23:38 GMT
References: <360@orawest.UUCP> <450005@hpcuhc.HP.COM> <12907@cup.portal.com> <727@islay.tcom.stc.co.uk>
Organization: The Portal System (TM)
Lines: 240

I've had enough requests for this that I'll go ahead and post it, 
rather than email'ing individual copies.  It's written in Borland
TurboBASIC, and would require extensive re-writing for MicroSoft
BASIC, but rather straight-forward for HP - I wrote it in HP BASIC
first... Anyway, the concepts are there, and I hope clear.

Oh yeah - it only concerns itself with sine-wave sources.

Dave


' This program will compute and graph the sidebands generated by
' frequency modulating a carrier oscillator, particularly by a modulating
' oscillator of near the same frequency as the carrier. The resulting
' waveform will also be displayed for reference.
'
' Two sideband plots are generated. The lower set is the mathematical
' result of the modulation, and shows sideband amplitudess in both positive
' and negative frequency domains. Amplitudes are shown as positive only,
' disregarding phase. The upper set is more practical, and shows only
' positive frequency components, with the components whose frequencies
' would have been negative "folded" across 0 (DC), and added in inverse
' phase to any positive-frequency components that happen to be at the
' same absolute frequency.
'
' The third plot is the result of directly mathematically modulating
' one sinusoidal signal with another.
'
' Some values to experiment with would be:
'
' Fc = 1 (all frequencies are relative - display is from -20 to +20)
' Fm = 1 (modulator at same frequency as carrier)
' Im = 1 (index of modulation)
'   This results in a distorted sawtooth wave, equivalent to summing
'   a considerable number of harmonics in additive synthesis.
'
'   Try reducing the Im to smaller values - the sawtooth wave will improve,
'   looking "best" somewhere around Im = .5
'
'
'
' Fc = 1
' Fm = 2 (modulator at twice the carrier frequency)
' Im = 1
'   This simple operator gives a square wave equivalent to 4 or 5 sine
'   generators added.
'
'   Again, reducing the Im removes some of the ripple, and results in
'   a smoother square wave.
'
'
'
' Fc = 1
' Fm = 1.4142 (modulator *not* in simple ratio to the carrier)
' Im = 2 (yields components extending farther out to the sides)
'
'   This shows some of the richness of components that can be generated
'   by FM synthesis, especially when the oscillators are not related simply.
'
'


dim amp(50), freq(50), ampf(50), freqf(50)

' Define graphics area
screen 8
window  (-20,-4)-(20,4)

fc=1
while fc<>0
  '  Get specs for the FM oscillators:
  locate 1,1
  input "Fc";fc  '  Center freq of (carrier) oscillator
  if fc>0 then
    input "Fm";fm  '  Freq of modulator
    input "Im";im  '  Modulation index (related to modulator output level)

    cls
    ' Get carrier oscillator frequency
    locate 15,50
    print "Fc = ";fc
    ' Get modulating oscillator frequency
    locate 16,50
    print "Fm = ";fm
    ' Get index of modulation
    locate 17,50
    print "Im = ";im

    amps=0
    s=1
    ' Number of significant sidebands depends on Im
    for sb=0 to im+3

      ' Set up arguments for Bessel function
      ord=sb
      arg=im
      gosub Bessel

      ' Calculate amplitudes and phase (0/180 deg)
      a1=j
      a2=j*s
      s=-s

      ' Set up arguments for component adder
      a=a1
      f=fc+sb*fm   ' Upper sideband
      gosub Add    ' Add into output bins
      '
      a=a2
      f=fc-sb*fm   ' Lower sideband
      gosub Add    ' Add into output bins
    next sb

    '  Set up axes
    gosub Axes

    '  Plot the sideband data
    gosub Bands

    '  Plot the waveform
    gosub Wave
  end if
wend

end

Bessel:'

'  Given parameters 'ord' and 'arg', returns the Bessel Function
'  J(arg), of order 'ord'

  m=18
  if arg=0 then x=1e-4 else x=arg  ' Can't handle x=0.0
  j1=0
  j0=1
  q=0
  for mm=m-1 to 0 step -1
    j=2*(mm+1)/x*j0-j1
    j1=j0
    j0=j
    if mm=ord then jj=j
    if mm/2=int(mm/2) then q=q+j0+j0   ' Add j0 on even passes only
  next mm
  q=q-j0
  j=jj/q   ' This (j) is the value we return
  return

Add:'

  '  Sums frequency/amplitude components into bins
  '  Inputs:  a, f: a frequency/amplitude pair

  done=0
  donef=0

  for i=0 to amps
    if freq(i)=f then amp(i)=amp(i)+a : done=1
  next i

  for i=0 to amps
    if freqf(i)=abs(f) then
      if f<0 then
        ampf(i)=ampf(i)-a
        donef=1
      else
        ampf(i)=ampf(i)+a
        donef=1
      end if
    end if
  next i

  if done=0 then
    amps=amps+1
    freq(amps)=f
    amp(amps)=a
  end if

  if donef=0 then
    if f<0 then
      freqf(amps)=abs(f)
      ampf(amps)=-a
    else
      freqf(amps)=f
      ampf(amps)=a
    end if
  end if

  return


Bands:'

  ' Plot the sideband data

  for i=0 to amps
    line(freq(i),0)-(freq(i),abs(amp(i))),15
  next i

  for i=0 to amps
    line(freqf(i),2)-(freqf(i),abs(ampf(i))+2),15
  next i

  return


Wave:'

  ' Plot the actual waveform from FM sine oscillators
  ' Inputs: Fc, Fm, and Im specify the carrier center frequency,
  '         the modulator frequency, and the modulator level

  line (-20,-2)-(-20,-2)
  for x=0 to 6.28318 step .0628318
    o1=im*sin(fm*x)
    o2=sin(fc*x+o1)
    line -(x/3.14159*10-20,2*o2-2),3
  next x

  return


Axes:'

    '  Set up axes
    line (-20,0)-(20,0),2
    for x=-20 to 20
      line(x,-.1)-(x,0),2
    next x

    line (-20,2)-(20,2),2
    for x=-20 to 20
      line(x,-.1+2)-(x,0+2),2
    next x

    line (0,-4)-(0,4),2
    for y=0 to 4
      line(-.5,y)-(.5,y),2
    next y

    return