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
- Introduction
- Installation
- User's guide
- Adding instruments or functions
- Scripts and interfaces
- The core and program layout
- References
Appendices
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 -
- 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 .
- ln -s $HOME/harmbal/tracpar . (discontinued)
- export PATH=$HOME/bin:$PATH
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]]
| -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:
where M = Kωr2/ωe2, R = Kgrωr/ωe2, and K = μeHωe2/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:
where![Zr(w)=i tan[w/4+(1-i)a(w)]](eqs/tubeimp1.png)
- Params: resfreq (resonance or reference frequency), nu (dimensionless loss parameter), and disper.
- Impedance function:
- 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
- Function:
- 103: Model 101 + vocal airways by Johnston et al.
- 104: Impedance of Guillemain et al. 2003
- 105: Approximate model 101
- 104: Impedance of Guillemain et al. 2003
- 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
- 115: Model impedance given by experimental values given by a file
- 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
Quantities applicable for many of the clarinet and saxophone models are summarised in the table below:
3.5 Examples of use:
- runs Harmbal on default parameter file, params.pmt, adjusting the initial frequency to 100 Hz, and return the solution, if found, in pout.pmt.
- 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.
- 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.
harmbal 100
harmbal -f clar1.pmt -o clar2.pmt -c gamma 0.5 -t 64 100 0 0.1 0 0
harmbal -f pout.pmt -c gamma 0.4
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 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.
- Remove the line double *Z;
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); - Assure that the line running initfunc() becomes the last of the three commands (just after the getparams() line) and change it to
- 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.
- Remove the line double *u;
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
- 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]
- 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)
- 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)
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.
- 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
- A simple saxophone model has been added.
- 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.