Snorre Ballière Farner

Harmbal User Guide

Harmbal is a free, open-source program in C under CeCILL, a French license compatible with the GNU General Public License. The last version can be downloaded from http://www.lma.cnrs-mrs.fr/logiciels/harmbal.

This user guide is updated for Harmbal version 2.1. See the Release notes if you upgrade.

The program has been realised with support from the MOSART IHP network of the EU, and was first presented at the MOSART Midterm Meeting in Esbjerg, Denmark, in August 2002. See References for theory about the method on which Harmbal is based.


Contents

  1. Introduction
  2. Installation
  3. User's guide
  4. Adding instruments or functions
  5. Scripts and interfaces
  6. The core and program layout
  7. References

Appendices
  1. Pitfalls
  2. Reporting bugs and improvements
  3. Release notes
  4. To be done

1  Introduction

Harmbal is a program in C for calculation of steady-state solutions to nonlinear physical models of self-sustained musical instruments by the harmonic balance method.

Harmbal consists of two parts, the core and a user part. The core part contains the harmonic-balance loop, the interface, numerics, and so on. The physical equations are programmed in the user part, a facility that enables the user to formulate and add new equations for new purposes without needing to worry about the details of the core of the program while obtaining fast calculation. A minimal knowledge of C is preferable, however.

Harmbal is written with financial support from the Europeen Union through the MOSART IHP network project at Laboratoire de Mécanique et d'Acoustique at CNRS, Marseilles, France by Postdoc Snorre Farner. A report of Harmbal was written for the mid-term meeting in Esbjerg in 2002 describing its principle [1]. However, the used is encouraged to refer to the more complete and recent paper by Farner, Vergez, Kergomard, and Lizée [3].

1.1  License

The software is governed by the CeCILL license under French law and abiding by the rules of distribution of free software. You can use, modify, and/or redistribute the software under the terms of the CeCILL license as circulated by CEA, CNRS, and INRIA at http://www.cecill.info.

As a counterpart to the access to the source code and rights to copy, modify and redistribute granted by the license, users are provided only with a limited warranty and the software's author, the holder of the economic rights, and the successive licensors have only limited liability.

In this respect, the user's attention is drawn to the risks associated with loading, using, modifying and/or developing or reproducing the software by the user in light of its specific status of free software, that may mean that it is complicated to manipulate, and that also therefore means that it is reserved for developers and experienced professionals having in-depth computer knowledge. Users are therefore encouraged to load and test the software's suitability as regards their requirements in conditions enabling the security of their systems and/or data to be ensured and, more generally, to use and operate it in the same conditions as regards security.

The fact that you are presently reading this means that you have had knowledge of the CeCILL license and that you accept its terms.

2  Installation

The program is only tested on Linux, but it does not use system-specific commands and should thus be compilable on any platform as long as there exists a C compiler.

Warning! If you upgrade from an older version and you have made changes to Harmbal, make a backup of these! Nonreported changes will not survive an update, and even if you have reported them, please take a backup to be sure.

2.1  Unix/Linux systems

Unzip and untar the package in a suitable directory, e.g. your home directory $HOME, by running either
tar xzvf harmbal-vX.XX-yymmdd.tgz
or
gunzip < harmbal-vX.XX-yymmdd.tgz | tar xvf -
This makes a new directory harmbal with all the files. Go to this directory and compile the program by running make. To be able to use Harmbal from other directories, make a soft link from your bin directory ($HOME/bin), i.e. go to this directory and run
ln -s $HOME/harmbal/harmbal .
ln -s $HOME/harmbal/tracpar . (discontinued)
ln -s $HOME/harmbal/hbmap .
ln -s $HOME/harmbal/cone .
ln -s $HOME/harmbal/normp .
(The periods at the end are necessary.) If $HOME/bin is not in your path (verify by echo $PATH), include it, e.g. in the bash environment:
export PATH=$HOME/bin:$PATH
Copy the parameter file params.pmt to the directory where you do your calculations.

2.2  Microsoft Windows

I have not tried this. Get yourself a C compiler (e.g. GCC) and follow its instructions. The file Makefile is probably usable with this compiler. If not, all .c files should be compiled to object files except tracpar.c and harmbal.c. These latter two are compiled when all the object files are available.

If you succeed in doing this, please tell me how you did it, and I will put it on this page. Thank you.

2.3  Macintosh

For Mac OS X, see the information for Unix/Linux. This has been tested. For Mac OS 9 or earlier, see the information for Microsoft Windows.

3  User's guide

Calculations with Harmbal are performed from the command prompt, thus a command-line window must be opened and the command harmbal followed by options and arguments must be executed. The parameters as well as a guess values of the spectrum and playing frequency are given in a parameter file (and possibly changed by the command-line arguments). The program returns the results in another parameter file of the same format so that a parameter can be changed and the program be re-executed based on the new parameter file.

This low-level interaction between user and program was chosen to give a high degree of flexibility; Harmbal can be used directly in connection with other programs, such as Matlab, or by shell scripts to do batch jobs such as finding solutions for a range of a parameter or iteratively searching a solution for a complicated case by starting from a simple one. And even if there should be a bug in Harmbal that causes it to crash, nothing is lost and the user may modify whatever she thinks caused the error and re-run Harmbal.

Finally, experienced C programmers may compile the functions of Harmbal into other C programs performing related tasks (see however licence agreement).

3.1  The usage of Harmbal

Harmbal is run by the command
harmbal  [options]  [frequency  [pressure harmonics]]
where brackets mean that the content is optional, options start with a minus sign (see examples), frequency is the guessed playing frequency, and partials are given as Re(P0), Re(P1), ..., Re(PNp), Im(P0), Im(P1), ..., Im(PNp). Possible options are:
-f  fn input name fn of parameter file (default: params.pmt)
-o  fn output name fn of parameter file (default: pout.pmt)
-p  Np change number of partials (excl. the constant component)
-t  Nt change number of time samples (should be a power of 2 and > 2Np
-c  p v change a parameter p to the given value v
-h help information
-v verbose iteration information
-i  fn use experimental data provided by a file fn for the resonator impedance, and set impmodel to 115
-b  fn use experimental impedance data provided by a file fn for an additional resonator impedance in parallel with the main resonator (given by the clarinet_tubeimp1 model (101), and set impmodel to 115 [this will be made more general in a later version]
-e  fn1 fn2 use experimental data provided by the files fn1 for the resonator impedance and fn2 for the additional impedance (see option -b), and set impmodel to 117 [this will be included in the -i and -b models in a later version]
-w writes Pc=ZU (poutPc.pmt) and Pm=Pc-P (poutPm.pmt) to file.
-M calculate until minimum of |G| is reached (default is only until |G| < maxerr)
-D direct Newton-Raphson, i.e. no backtracking
-F filter u by lowpass before FFT
-I show final impedances at end of calculations
-J write final Jacobian matrix to file (jacob.dat)

3.2  Some rules

  • If input file already is a solution, the result is merely copied to the output file. If Harmbal cannot find a solution, the output file is not made (nor touched if already existing).

  • The command-line arguments take precedence to the values in the file, and if Np is given, it is more important than the number of partials given in the file or even at the command line. Thus, there is no point in specifying at the same time Np and the partials.

  • The frequency must be given if the partials are given, and the partials must include the constant component, so for Np partials, you should have 2(Np+1) values plus the frequency.

3.3  The parameter and other files

3.3.1  The parameter file

The parameter file is simply a text file that lists all necessary parameters and their values, one per line in the format parameter value. The last item in the file must be P which is given by P alone on a line immediately followed the 2(Np+1) partials, one on each line. Blank lines and lines starting with a hash (#) are ignore and not rewritten in the output parameter file.

There are some parameters that that must be supplied: exmodel (or K, R, and M, for the default exmodel=101; see Equation (10) in the report), nlmodel and impmodel (nonlinear and impedance model numbers), pmax (=pM, a value of 1 for dimensionless), denom (which pressure compontent harmonic (default is 1, i.e. Re(P1)) to use in the denominator in Equation (13)), maxitno (maximum number of iterations), maxerr (maximum |G| and |f i - f i-1| /f i acceptable for convergence), freq (initial playing frequency f 0), resfreq (the lowest resonance frequency of the resonator ωt/2π), and of course Np and Nt. Both model numbers are the sum of the instrument number multiplied by 100 and the model number. For the nonlinear model (model 2) of the clarinet (instrument 1), for example, one line in the parameter file should read nlmodel 102.

Other parameters depend on the nonlinear and impedance models used. The use of the equations in Section 2.3 in the report requires nu (damping factor) and disper (dispersion flag; 0=none, 1 dispersion; with intermediate values allowed if needed to obtain convergence) for the tube, and zeta (ζ) and gamma (γ) for the nonlinear equation. An example.pmt: K 1.0000000000 R 0.000000000000e+00 M 0.000000000000e+00 zeta 2.000000000000e-01 gamma 3.100000000000e-01 nu 2.000000000000e-05 resfreq 100.0000000000 lenfact 2.0000000000 disper 0.000000000000e+00 pmax 1.0000000000 maxitno 200.0000000000 denom 1.0000000000 maxerr 1.000000000000e-06 Nt 512 Np 2 exmodel 101 nlmodel 102 impmodel 201 freq 99.9999999980 P 5.434994719237852203e-24 1.975789892011730609e-01 -1.163266761622221424e-01 0.000000000000000000e+00 0.000000000000000000e+00 8.680410426329590981e-11

3.3.2  Other files that Harmbal writes

Harmbal writes a few other helpful files:
  • P.dat: Contains the harmonics detached from the .pmt file.
  • uswept.dat: Contains 3 columns: p, x, and u for the period of Nt samples.
  • err_P.dat and err_uswept.dat: same as P.dat and uswept.dat, but they are returned in case of lack of convergence or other error. They are very useful if convergence is stopped due to a too low maxitno! (Ctrl-C doesn't make these files to be written.)

3.4  Instrument models

Currently, Harmbal incorporates instrument models for the clarinet (100) and for the saxophone (200). They are listed below, exciter-impedance models first, then resonator-impedance models, and finally nonlinear-coupling models. The variables are described below.

3.4.1  Exciter impedances exmodel

Clarinet impedence models
101: Simple spring impedance
Impedance function:
Ze(w)=K-M(w/2pi)^2+R(w/2pi)
where M = r2/ωe2, R = Kgrωr/ωe2, and K = μee2/pM
Params: resfreq = ωe/2π, M, R, K

3.4.2  Resonator impedances impmodel

[Details to be added]

Clarinet impedence models
101: General cylindrical tube impedance
clarinet_tubeimp1()
Impedance function:
Zr(w)=i tan[w/4+(1-i)a(w)]
where a(w)=psi*eta*sqrt(w/2pi)
Params: resfreq (resonance or reference frequency), nu (dimensionless loss parameter), and disper.
102: Model 101 + vocal airways by Claudia Fritz
clarinet.c: clarinet_totalimp()
Function:
Params: impmodel, resfreq, nu, disper, f1, f2, A1, A2, Q1, Q2, sndspeed
103: Model 101 + vocal airways by Johnston et al.
104: Impedance of Guillemain et al. 2003
105: Approximate model 101
Function: Z(2n) = eta*psi*sqrt(2n), Z(2n+1)=1/(eta*psi*sqrt(n))
107: Model of Raman
115: Model impedance given by experimental values given by a file
116: Model + vocal-airways impedance given by experimental values
117: Model 115 + 116

Saxophone impedance models
[To be added]

3.4.3  Nonlinear coupling functions nlmodel

[Details to be added]

Clarinet nonlinear models
[To be added]
Saxophone nonlinear models
[To be added]

3.4.4  Variable and constant descriptions

Originally, Harmbal was written to use dimensionless variables and equations. From version 1.30, the exciter equation can be changed freely, opening for using Harmbal with dimensional variables. To avoid misunderstandings, we use the following conventions for naming the variables. For historical reasons, variables are in originally dimensionless. Physical (dimensional) variables are then marked with a hat (e.g., â as in Farner et al. [3]) or a star (e.g., a* or astar), as a hat is not so simple to apply in HTML. Constants may be dimensional or dimensionless without need for a special mark. See also Farner et al. [3] for a complete description of the method and for equations related to the clarinet.

The three equations for the exciter, the coupling, and the resonator are coupled by three principal dimensionless variables xe(t), xc(t), and xr(t) with related Fourier variables Xe(ω), Xc(ω) , and Xr(ω). For the clarinet, the time domain variables are related to the following physical quantities:

  • Exciter variable: xe = x = y*/H + γ/K where y* is the reed opening
  • Coupling variable: xc = p = p*/pM where p* is the dynamical pressure in mouth piece
  • Resonator variable: xr = u = u*ρc/SpM where u* is the volume flow through mouth piece
The dimensionless time is defined as t = t*&omegar, which is important to keep in mind when taking the Fourier transform of an equation.

Quantities applicable for many of the clarinet and saxophone models are summarised in the table below:

Table of quantities

3.5  Examples of use:

harmbal 100
runs Harmbal on default parameter file, params.pmt, adjusting the initial frequency to 100 Hz, and return the solution, if found, in pout.pmt.
harmbal -f clar1.pmt -o clar2.pmt -c gamma 0.5 -t 64 100 0 0.1 0 0
runs Harmbal on an earlier parameter file clar1.pmt, changes gamma to 0.5, Nt to 64, the frequency to 100 Hz and the initial P to have one harmonic with amplitude 0.1, makes the appropriate iterations to find a solution, and finally writes the result to clar2.pmt.
harmbal -f pout.pmt -c gamma 0.4
reads pout.pmt, changes gamma, runs Harmbal, and writes to the same file. If no solution was found, pout.pmt is not changed. This is particularly practical when a gradual change of a parameter is necessary, or to collect solutions for a range of a parameter.

3.6  Error messages:

If harmbal is not capable of converging towards a solution, it gives an error message:
  • 1: General error
  • 2: No convergence within max number of iterations
  • 3: Frequency negative
  • 4: Minimum lambda reached too many times
  • 5: Solution explodes
If an error should interrupt the calculations, the last solution will be available in the three files err_pout.pmt, err_uswept.pmt, and err_P.dat. This is useful in at least two cases:
  • If you have chosen maxitno too small and you want to continue the iteration, simply rerun with harmbal -f err_pout.pmt
  • If you would like to know u(t) or Pi of a theoretical solution that you put into theory.pmt, then run

    harmbal -f theory.pmt -c maxitno 0

4  Adding instruments or functions

To add a nonlinear equation u or an impedance Z, it is necessary to go into the source code of the program and do some simple programming. However, this part of the code is separated from the core of the program to facilitate the change for unexperienced users of C. This is described in the following sections. Section 3.3 lists some pitfalls.

After whatever little change of the program, it must be recompiled. This is done in Unix/Linux systems by writing make in the directory of the program. If there are errors, these are often of the nature covered in the Section 3.3.

Originally, Harmbal was built for dimensionless equations and thus based on dimensionless parameters. In this case you should ensure that all equations employ the same dimensionless versions of the three variables xe (exciter variable x = y'/H + γ/K for the clarinet), xc (coupling variable p = p'/pM), and xr (resonator variable u = u'ρc/SpM), with respective Fourier transforms (Xe, Xc, and Xr). This applies to the choice of characteristic impedances Zrc and Zec. The prime (') or hat (^) denotes the corresponding dimensional quantity. See the references for the theory about the equations.

From version 1.30, it became possible to change the exciter equations as well as the resonator and coupling equations. This implies that dimensional equations may be used as well. Be carefull that all your equations use the same set of variables.

4.1  Adding a new instrument

Assume that you want to add saxophone functions to the program. Then you should make a new file saxophone.c with a header file saxophone.h (see below). But first, tell about its existence by opening the file named instr.h and adding a line #include "saxophone.h" to the list close to the end of the file: #ifdef INSTR #include "clarinet.h" #include "saxophone.h" #endif

Then decide a number that is not already taken, say 2, for your new instrument package and add a line #define SAXOPHONE 2 at the end of the file. Within the program, you should always refer to this instrument by the constant SAXOPHONE, not the number, in case it is necessary to change the number.

In the file instr.c, find the function general_resonator() and add a new case in the switch clause, just before the default case. This is where the instrument is chosen depending on the parameter impmodel (or nlmodel for the nonlinear coupling model and exmodel for the exciter model): case SAXOPHONE: reson = saxophone_resonator(impmodel % INSTRFAC,paramlist); break;

Repeat this in the functions general_exciter() and general_nonlin().

Now, it is time to make the file saxophone.c by copying the standard file stdinstr.c. Substitute saxophone for clarinet and update saxophone_exciter(), saxophone_resonator(), or saxophone_nonlin() every time a new function is added, see next sections. Make also a file saxophone.h (from stdinstr.h) containing the declaration of all the new functions you add.

In makefile you should add saxophone.o to the variable INSTOBJS, which you find in the preamble of the file. This tells of the existence of the new file to the compiler.

4.2  Adding a new exciter

To add a new model for the clarinet exciter, a new function must be created and put into clarinet.c. The function name should start with the name of the instrument file to avoid double-use of a name. The framework of an exciter-impedance function is illustrated by the function for the simple spring with mass and damping: /* Mx'' + Rx' + Kx = p; simple spring */ int clarinet_simplespring(double *H, int Nc, double freq, double *params) { static int k; static complex tmp; static double resfreq,M,R,K; // initialize each param variable // give a value to all parameters resfreq = params[0]; M = params[1]; R = params[2]; K = params[3]; freq=freq/resfreq; for(k=0;k<Nc;k++){ tmp = Cinv(Complex(K - M*pow(k*freq,2),R*k*freq)); H[k] = tmp.re; /* H = Y/P = 1/(-Mw^2 + iRw + K) = imp */ H[k+Nc] = tmp.im; } return 0; }

For Harmbal to recognize the new fuction, you must add its function declaration to clarinet.h, i.e. simply add the first line of the function followed by a semicolon. Leave also a little explication above as has been done for the other declarations: /* Mx'' + Rx' + Kx = p; simple spring */ int clarinet_simplespring(double *H, int Nc, double freq, double *params);

Close to the top of the file clarinet.c, there should be a function clarinet_exciter(). The new resonator must be added to the switch clause, just before default. Here the function to be used is chosen by the parameter exmodel. We look at the already existing case 1 in clarinet.c: case 1: // Ze(w)=K-M(w/2pi)^2+R(w/2pi) np = 4; params = getparams(paramlist,stringlist(np,"resfreq", "M", "R", "K"),YES); exciter = initfcnstruct(clarinet_simplespring, params, np); break; Text in red should be verified and changed if necessary. When you have made a new function, all the parameters that it needs should be listed with quotes as arguments to the getparams() function, and np in the line above should be changed to the number of parameters. It is important that resfreq is the first as this is used by other functions in the program. In the line starting with exciter, substitute the name of your new function. To write the impedance function as a comment after case makes it easier to find the right case later. The case number is the model number and should thus be changed accordingly.

Finally, all new variables should be added to the parameter file used for the problem. The exmodel parameter is set to the number after case. The Missing parameters will produce an error message, but they may also be defined at the command line in recent versions. For backward compatibility, a missing exmodel makes Harmbal assume that the original function 101 is chosen.

4.3  Adding a new resonator

To add a new clarinet resonator, a new function must be made and put into clarinet.c. The function name should start with the name of the instrument file to avoid double-use of a name. The framework of a resonator-impedance function is as follows: int clarinet_tubeimp(double *Z, int N, double freq, double *params) // Z=tanh[j*2*pi*f/(4*resfreq) + att*sqrt(f/resfreq)] // write the function clearly here { // declare all local variables int i; // harmonic counter, 0=constant component double att,resfreq; // parameter names complex Zk,arg; // other temporary complex... double fratio; // ...and real variables // put understandable names to the input parameters resfreq = params[0]; // resonance freq. of full tube length att = params[1]; // attenuation /* calculate Z */ fratio = freq/resfreq; // angular frequency of first harmonic for(i=0;i<N;i++){ // calculations, dimensionless and in complex form arg.re = att*sqrt(k*fratio); arg.im = 0.5*pi*k*fratio; Zk = Ctanh(arg); Z[k] = Zk.re; // convert to real array Z[k+N] = Zk.im; } return 0; }

For Harmbal to recognize the new fuction, you must add its function declaration to clarinet.h, i.e. simply add the first line of the function followed by a semicolon: int clarinet_tubeimp(double *Z, int N, double freq, double *params);

Close to the top of the file clarinet.c there should be a function clarinet_resonator(). The new resonator must be added to the switch clause, just before default. Here the function to be used is chosen by the parameter impmodel. We look at the already existing case 1 in clarinet.c: case 10: // Z=tanh(jkl + al) np = 2; params = getparams(paramlist,stringlist(np,"resfreq","att"),YES); reson = initfcnstruct(clarinet_tubeimp, params, np); break; Text in red should be verified and changed if necessary. When you have made a new function, all the parameters that it needs should be listed with quotes as arguments to the getparams() function, and np in the line above should be changed to the number of parameters. It is important that resfreq is the first as this is used by other functions in the program. In the line starting with reson, substitute the name of your new function. To write the impedance function as a comment after case makes it easier to find the right case later. The case number is the model number and should thus be changed accordingly.

Finally, all new variables should be added to the parameter file used for the problem. The impmodel parameter is set to the number after case. The Missing parameters will produce an error message, but they may also be defined at the command line in recent versions.

Note for upgrading from version 1.28 to 1.29:

Assuming that you have added a function to clarinet.c for version 1.28 of ealier, the following steps must be followed for the changes to be accepted by version 1.29 or higher:

In clarinet.h:
Change the declaration lines for impedance functions: Exchange double* by int and add double *Z as first argument, e.g.: int clarinet_tubeimp(double *Z, int N, double freq, double *params);
In the clarinet_resonator() function in clarinet.c:
Use the function name initfcnstruct instead of the old initresonator, e.g.: reson = initfcnstruct(clarinet_tubeimp, params, np);
In the impedance functions in clarinet.c:
Change the first line as in the *.h above, e.g.: int clarinet_tubeimp(double *Z, int N, double freq, double *params)
Remove the line double *Z;
Remove the line Z = allocvec(2*N);
Change the return argument from Z to 0 in the last line.

4.4  Adding a new nonlinear function

Adding a new nonlinear function is done in the same way as adding a new resonator, i.e., a nonlinear function that corresponds to a clarinet, should be put into the file clarinet.c. The function name should start with the name of the instrument file to avoid double-use of a name. The framework of a nonlinear function is as follows: int clarinet_nonlinmodel1(double *u, double *x, double *p, int Nt, double sampfreq, double *params, int np) // u=Ax+Bp+dx/dt { int n; double A,B; // parameters double deriv; // temporary variables A = params[0];// get parameter from params B = params[1]; // calculate dimensionless u nold = Nt-1; for(n=0; n<Nt; n++){ deriv = (x[n]-x[nold])*sampfreq; // the time derivative of x at n u[n] = A*x[n] + B*p[n] + deriv; // arbitrary example function! nold = n; } return 0; }

Similarly to the impedance, the new function must be made known to Harmbal by repeating the first line of the function in clarinet.h, with a trailing semicolon: int clarinet_nonlinmodel1(double *u, double *x, double *p, int Nt, double sampfreq, double *params, int np);

The function must also be added to the clarinet_nonlin() function in clarinet.c as a case in the switch clause, for instance just before the default statement: case 10: // u=Ax+Bp+dx/dt np = 2; // i.e., "A" and "B" params = getparams(paramlist,stringlist(np,"A","B"),YES); nonlin = initfcnstruct(clarinet_nonlinmodel, params, np); break; Text in red should be verified and changed if necessary. When you have made a new function, all the parameters that it needs should be listed with quotes in the getparams() function, and np in the line above should be updated to the number of parameters, A and B in this case. The number after case must be changed so that it is unique. This is the model number referred to by nlmodel when running Harmbal. In the line starting with nonlin, substitute the name of your new function. To write the nonlinear function as a comment after case makes it easier to find the right case later.

Finally, all new variables should be added to the parameter file used for the problem, and the nlmodel be changed to the number after case. Missing parameters will produce an error message, but they may also be defined at the command line in recent versions.

To update from version 1.27 to 1.28:

Assuming that you have added a nonlinear function to clarinet.c for version 1.27 of ealier, the following steps must be followed for the changes to be accepted by version 1.28 or higher:
In clarinet.h:
Change the declaration lines for impedance functions: Exchange double* by int and add double *u as first argument, e.g.: int clarinet_nonlinmodel(double *u, double *x, double *p, int Nt, double sampfreq, double *params, int np);
In the clarinet_nonlin() function in clarinet.c:
In the line invocing the getparams() function, change the variable name valuelist to params, e.g.:
params = getparams(paramlist,stringlist(np,"A","B"),YES);
Assure that the line running initfunc() becomes the last of the three commands (just after the getparams() line) and change it to
nonlin = initfcnstruct(clarinet_nonlinmodel1, params, np);
In the nonlinear functions in clarinet.c:
Change the first line as in the *.h above, e.g.: int clarinet_nonlinmodel(double *u, double *x, double *p, int Nt, double sampfreq, double *params, int np)
Remove the line double *u;
Remove the line u = allocvec(Nt);
Change the return argument from u to 0 in the last line.

5  Scripts and interfaces

5.1  Hbmap for continuation

Hbmap is a Perl script that runs Harmbal automatically for a range of a given parameter. The usage is: hbmap [options] outfile.dat param end [steps [start]] which changes param from its present value (found in the file pout.pmt or from the optional argument start, to end. Unless start is given, the starting point is excluded. The optional argument steps is the number of steps between start and end points, and if steps 0 (zero) or not given, Harmbal uses a variable step, i.e. it chooses 10 steps and decreases the step if a solution is not found. This is useful to find the end of a regime or to pass an area with bad convergence, even though this gives no garantee. All necessary parameters are taken from pout.pmt.

Hbmap writes one line of results for each step to outfile.dat. Default output is:

param freq P1 Re(P2) Im(P2) Re(P3) Im(P3)...

where param is the value of the parameter that changes, freq is the resulting playing frequency, and P1 etc are the (complex) pressure Fourier components. A few comments (lines starting with "#") are made in the file in case manual manipulation should be necessary. The data file is directly plottable using Gnuplot, for instance.

Options:

-a  print absolute values instead of giving both Re and Im:
 param freq |P1| |P2| |P3|...

5.2  Normp for normalizing P

normp is a Perl script for changing the parameter pmax and correspondingly modify the components of P, see The parameter file.

Usage: normp [pmax] < infile.pmt > outfile.pmt where the new file outfile.pmt is the same list of parameters as the original file infile.pmt, except that pmax and the components of P are modified. They are all divided by the value given either by the argument pmax or by the parameter pmax extracted from infile.pmt if the argument is not provided by the user. The latter case is the normalization case in which P is divided by pmax. The "<" and ">" are the Unix redirection symbols, which should also work for MS Windows.

5.3  Cone for theoretical stepped cones

cone is a Perl script to calculate the Np first harmonics of the Helmholtz (square-wave) oscillation for a clarinet-like nonlinearly coupled reed and resonator, where the resonator is a lossless stepped cone (see Kergomard et al., Acustica-Acta Acustica 2000 (86), pp.685-703) with N steps, N = 1 representing a cylinder.

Usage: cone Np {g gamma|p pmax} [> file.pmt] where p is the smallest amplitude in time domain while gamma is the wanted γ at which to calculate p. The output is in the Harmbal parameter-file format written to standard output (or optionally redirected to the file file.pmt).

5.4  Tracpar for tracing the minimizing function

First, Tracpar is not a very robust program, and no longer updated with the evolution of harmbal. Some modification in the code may therefore be necessary. It was made for debugging purposes, where the function to be minimized, G1(P), was traced with respect to P1 for Np= 1.

Usage: tracpar [-f infile.pmt] [-o outfile.dat] i param from to [steps] Tracpar is a test program for tracing Gi(P) or |G| versus a given parameter (param), for example γ (gamma) or P1 (P1) in the range from-to in steps steps (default 100). The rest of the parameters are taken from infile.pmt (default pout.pmt), and the results are written to outfile.dat (default tracpar.dat). The code must be slightly modified in order to trace |G| instead of Gi(P).

The output is written in three columns: the parameter value, Gi(P) or |G|, and finally its derivative with respect to the parameter.

6  The core and program layout

6.1  Files in the package

You should have received all the following files:
README                   # general and installation information
makefile                 # the compilation program
LICENSE-FR               # license agreement in French
LICENSE-EN               # English translation of license agreement

main.c                   # main program
tracpar.c                # assisting program for tracking a parameter
harmbal.c, harmbal.h     # harmonic balance loop with subroutines

calc.c, calc.h           # calculation subroutines
interface.c, interface.h # file and screen interface subroutines
mem.c, mem.h
instr.c, instr.h         # main instrument

clarinet.c, clarinet.h   # instrument-specific functions
saxophone.c, saxophone.h #
stdinstr.c               # template for instrument files

params.pmt               # example parameter file for saxophone

hbmap                    # Perl script to cover a range of a parameter
normp                    # Perl script to normalize P in a parameter file
cone                     # Perl script to calculate theoretic P components

7  References

  1. S. Farner. Harmbal: Open-source computer program in C for calculating steady-state solutions to non-linear physical models of wind and string instruments. MOSART Midterm Meeting, Esbjerg, August 2002 [PDF]
  2. C. Fritz, S. Farner, and J. Kergomard. Some aspects of the harmonic balance method applied to the clarinet. Applied Acoustics, 65 (12), pp. 1155-1180 (2004)
  3. S. Farner, C. Vergez, J. Kergomard, and A. Lizée. Contributions to harmonic balance calculations of periodic oscillation for self-sustained musical instruments with focus on single-reed instruments. J. Acoust. Soc. Am., 119 (3), pp. 1794-1804 (2006)

Appendices

A  Pitfalls

  • Since we use dimensionless quantities, the impedance should be divided by the characteristic impedance, for instance Zc=ρc/Ac, where Ac is the characteristic cross section of the tube. The definition of Zc is related to ζ. This is easy to forget when using more than one impedance function at the same time.

  • Arrays in C starts with the zeroth element, so for Nt elements, the elements of u are u[0], u[1],...,u[Nt-1].

  • A forgotten semi-colon after a statement or mismatched braces ({, }) may cause errors somewhere else.

  • All variable must be declared (floating point variables as double, integers as int, and complex variables as complex.

  • Calculations with complex variables must be performed with the functions given in calc.h, even adding or changing sign (since complex is a struct), e.g. aeiz+b may be written RCmul(a,Cexp(Cadd(ICmul(1,z), Complex(b,0)))), where a and b are double and z is complex. To access to the real and imaginary parts, add .re and .im, respectively.

B  Reporting bugs and improvements

First, the program is non-interactive, so should it crash, an error message is normally given, and the program stops with no loss of data. Correct the error and rerun.

However, please do not hesitate to send me an e-mail if you cannot solve a problem after having consulted this guide. Then you have probably found either a bug in the program or a point to be improved in this guide.

I also encourage everyone to send me bug reports and improvements to the code, and new model implementations. I will release it in a new version. In this way everyone can take profit of your improvements, and Harmbal starts to live its own life.

Please include:

  • version number
  • occured behaviour and/or error messages
  • command-line arguments and parameter file that provoked the error
  • program files that have been changed (if applicable)
If you have implemented a new model, I appreciate some sort of documentation along, to include in this guide.

Thank you.

C  Release notes

Version 2.1 (2008-02-06):
  • Added possibility to write finale Jacobian to file (option -J)
Version 2.0 (2007-08-13):
  • Bug-fix: Arbitrary values were reported for u in uswept.dat
  • Harmbal is now under CVS to facilitate parallel development.
Version 1.31c (2007-05-28): Bug fix for compiling with gcc 4
  • warning upon compilation: incompatible implicit declaration of built-in function
  • wrong numbers stated after "Lin.diff" at the end
  • segmentation fault when using -i after "Resonator" at the end. Now no number is given after the resonance frequency.
Other changes:
  • numbers listed at the end easier to read
Version 1.30 (2006-06-19):
  • The exciter is now a modular unit, where an arbitrary impendance function can be given, similarly to the resonator. Add exmodel to the .pmt file (or by -c). The original model Mx''+Rx'+Kx=p is available as exmodel=101. Harmbal is backward compatible and assume this model if exmodel is not given.
Version 1.29 (2006-03-05):
  • Avoid allocating memory for Z each time calcimp() is called. This implies that if you upgrade from 1.28 or earlier and your impedance functions are not in the official distribution, you need to make a few changes as explained in Sec. 4.2.
Version 1.28 (2005-05-11):
  • Avoid allocating memory for u each time calc_u() is called. This implies that if you upgrade from 1.27 or earlier and your nonlinear functions are not in the official distribution, you need to make a few changes as explained in Sec. 4.3.
  • New struct fcnstruct that contains all necessary information for each of the three equations. A later update will makes it possible to change the exciter equation also.
Version 1.27c (2005-02-08) with contributions from M. Demoucron/C. Fritz:
  • Added option -i for using experimental impedance data provided by a file (sets impmodel 115).
  • Added option -b for using an additional impedance provided by a file with clarinet_tubeimp1 (sets impmodel 116).
  • Added option -e for using both -i and -b
  • Added option -w that saves Pm and Pc in pmt-files
  • Added option -I to show the final impedances when the program finishes.
  • Version number appears now upon harmbal -h
  • Removed comment line in output data files for easier access in Matlab.
Note: In a later version, -e should be a result of -i + -b, and -b should be possible to combine with whatever resonator impedance.

Version 1.26 (2005-01-28):
  • License had to be changed from to CeCILL, so this is done.
  • It is now possible to define new parameters on the command line without defining them in infile.pmt first.
  • Added clarinet_raman() to clarinet.c
  • Renamed isnumber() to isnumb() in interface.c since isnumber() is defined by ctype.h, which also is used in interface.c.
Version 1.25 (2004-02-27):
Mostly minor improvements during testing and use.
A simple saxophone model has been added.
An optional low-pass filter has been put before the FFT (option -F).
Added the Perl scripts cone and normp
Added clarinet_imp_guillemain() to clarinet.c
Version 1.20 (2002-09-26):
First version available. No known bugs.

D  To be done

  • Support for dimensional equations (should be possible already, especially starting with version 1.30, by exchanging equations and quantities by dimensional counterparts)
  • Update the experimental impedance models so that -e is a result of -i + -b (thus making -e obsolete), and that -b can be combined with whatever resonator impedance.
  • Put all vectors in a single struct and let it be available to all functions needing it instead of sending a number of vectors in all function calls.