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. 