Time Domain Model of the Flute CURWAVL.BAS - John W. Coltman

     The BASIC program CURWAVL.BAS models the currents in the flute
as a function of time.  Time (U%) proceeds in steps called ticks (
% with a variable indicates an integer).  The variable T% is a
local time running from 0 to 63, and a circular array 64 ticks long
contains past values of jet current, J1(T%), pipe current P(T%)
consisting of a downgoing part P2(T%) and an upgoing part (P1(T%),
and mouth current M(T%).  The arrays are updated at each tick with
a new value, and the value 64 ticks ago is dropped.  At any
instant, the mouth current is the sum of the pipe current and the
jet current.  The flute is assumed to be so long as to require R%
ticks for a round trip of the wave (I set this to 32 in line 40,
there is no need to change it).  At the far end is a reflection
function that is convoluted with 11 past values of the downgoing
pipe current to get the returning pipe current. The resulting mouth
current is used as a driver to deflect the jet, and the jet current
in turn is determined by a delay J% and a saturation function. The
jet current entering the flute splits into two parts, one going
down the pipe and the other (algebraically negative, since the
particle velocity is in the opposite direction) up, where it is
reflected by the end correction (takes three ticks round trip) and
follows the other part down the pipe.  Both of these are added to a
downgoing wave reflected from the up going wave.

     The program presents graphically the pipe current wave (above)
and the jet current wave (below) as they develop.  Time is tallied
in the left upper corner, and  points of upward zero crossings at
the right.  These are useful in deciding how many periods to
analyze later. It may take several thousand points for the wave to
settle into a steady form - longer for some critical cases. At any
time, pressing the space bar will interrupt the process, and a menu
is presented for analysis. One may display the numerical values for
the last 42 points, Fourier analyze the pipe wave, print out
values, or continue the simulation from the point last tallied. In
the Fourier analysis, the program selects points between successive
upward zero crossings, fits a curve to these, determines and
displays the period in ticks that the wave has settled down to
(need not be an integral number of ticks), and interpolates values
for 32 points.  This curve is displayed for inspection - sometimes
the curve fitting doesn't go properly, and one can check if
everything is OK before proceeding. If not, continue the wave for a
bit and try again. If there are two zero crossings per cycle, as
may occur with high harmonic content, one should select 2 periods
when asked this question before the analysis. A bar graph showing
the harmonic content compared to the largest component is
displayed, followed by a display of numerical values for the
harmonic components.

     Variables at the disposal of the experimenter are the jet
delay J%, which must have integral values form 0 to 63, the gain G
on the jet, which determine how far it swings in response to the
mouth current, the jet offset S (value zero for centered, up to 1
or 2 for large offsets).  The saturation jet current C1 may be left
at 1, it simply sets the scale of the output.

     Here I will discuss the important program lines by number.
Lines 10-30 are housekeeping. Line 40 is the approximate period of
the fundamental for the pipe in ticks, it should probably be left
at 32. Line 50 just sets the scale of things, leave it at 1. Line
60 , the jet gain, is an important variable, and should be changed
to explore various cases from weak jet deflection (perhaps .3) to
strongly saturated cases (2). Line 70, delay on the jet, is a
second important variable, representing in part the blowing
velocity. Only integer value are permitted here. As the delay
decreases, the period of oscillation shortens slightly, at some
point it will break over into the second mode, cutting the period
in half.  Line 80 is the offset of the jet, which brings in
interesting behavior with respect to harmonic generation.  Lines
60, 70 and 80 are the major ones you should change to explore the
behavior of the model.

     Lines 90 and 100 store the values of the reflection function.
The jagged nature of the points given in line 100 is due to the
fact that they have been weighted by iterated 1,4,1 weights
according to Simpson's Rule for finding areas to a 2cd degree
approximation.  Without these weights the  values would have been
0,-.85,0,2.5,1.65,.75,.4,.25,.15,.1,.1 a curve that dips negative,
rises quickly to a peak and falls off slowly. This curve was
arrived at by trial and error plus some analytical derivations of
of what reflection functions look like for inductive terminations.
It can be justified by changing line 190 in the program to
J1(T%) = SIN(2*PI*U%/PE).  This will disable the feedback and
substitute a sinusoidal forcing function for the jet current with
period PE. By choosing various values of PE one can demonstrate the
resonance responses of the pipe. For the reflection function chosen
the response is very like a flute, with the first 4 modes having
Q's in the neighborhood of 50, and slightly stretched from an exact
harmonic relationship. The entire curve of pipe current as a
function of frequency resembles very closely that which one gets
from an analytical derivation using a simple transmission line
analogy.

     Line 110 fills the downgoing wave array P2 with a sinusoidal
function to get the loop off to a fast start.  Without this it may
a take a long time to build up, or may just stall. For cases where
oscillation in more than one mode is possible, one can change the
number in the denominator from 32 to say 16 or 10, to encourage the
second or third mode to sound.

     Line 130 starts the loop for the simulation, allowing up to
9999 ticks, which is enough for almost any case to settle down.
Line 140 sets the local time, which is the actual time mod 64. Line
150 represents the jet displacement V at the mouth hole, being
proportional to the mouth current at an earlier (by the jet delay
J%) time, the jet gain G, and offset by an amount S1.

     The next four lines, 160 to 190, produce the jet current from
the jet deflection. For deflections >1 the current saturates at 1,
for deflections less than -1 the current has a value 0.  In between
the function follows a cubic that is linear near the origin, falls
off in slope and joins the saturation values of 1 and 0 at
deflections 1 and -1 with no discontinuity in slope, and has a
value of .5 at zero deflection. Line 190 allows one to change the
scale of the jet current.

     Line 200 starts the convolution by setting the initial area S
to zero. Line 210 performs the multiplication of the reflection
function C(F1%) with points on the downgoing wave P2, incrementing
the area S at each point.  The index of P2 depends on the present
time U%, diminished by the round trip time R%, increased by 8 to
center the convolution curve properly, and indexed by the local
time variable F%. In Line 220 the reflected (upgoing) wave is given
by the sum S normalized by dividing by 17.03, an important part of
the reflection function that sets the resistive (frequency
independent) losses and therefore in part the Q.

     Line 230 gives the downgoing wave. It is made up of the
upgoing wave reflected from the upper end of the flute after a
delay of three ticks, plus half of the jet current of the moment,
plus the other half reflected from three ticks earlier. This half
of the jet current is given a negative sign because it is directed
initially in the negative x direction. Line 240 is the total pipe
current, the sum of the upgoing and down going pipe currents.  Line
250 is (temporarily) the mouth current, the pipe current less the
inwardly directed jet current. Line 252 and two following is a DC
filter applied to remove any steady flow from the mouth current so
as not to deflect the jet, leaving the offset at the control of the
experimenter rather than permitting the jet to be influenced by DC
terms. It is made up of 99% of the mouth current from the previous
loop, plus 99.5% of the difference in the new driving function and
the previous driving function. You can try this out in a separate
program to see how it acts as a high pass filter, eliminating DC
but hardly affecting any frequencies encountered in the simulation.
This filtered mouth current M(T%) is the value used in line 150 in
the next loop.

     Lines 260-330 are used to display the time, plot the two
waves, save recent zero crossings, and permit interruption of the
loop. Line 340 loops back to calculate the next set of values. All
of the rest of the program is for the subsequent analysis and
display of results, and need not be discussed here.