1HARMINV(1) harminv HARMINV(1)
2
6 harminv - extract mode frequencies from time-series data
7
9 harminv [OPTION]... [freq-min-freq-max]...
10
12 harminv is a program designed to solve the problem of "harmonic inver‐
13 sion": given a time series consisting of a sum of sinusoids ("modes"),
14 extract their frequencies and amplitudes. It can also handle the case
15 of exponentially-decaying sinusoids, in which case it extracts their
16 decay rates as well.
17 harminv is often able to achieve much greater accuracy and robustness
19 than Fourier-transform methods, essentially because it assumes a spe‐
20 cific form for the input.
21 It uses a low-storage "filter-diagonalization method" (FDM), as
23 described in V. A. Mandelshtam and H. S. Taylor, "Harmonic inversion of
24 time signals," J. Chem. Phys. 107, 6756 (1997). See also erratum, ibid
25 109, 4128 (1998).
26
28 harminv reads in a sequence of whitespace-separated real or complex
29 numbers from standard input, as well as command-line arguments indicat‐
30 ing one or more frequency ranges to search, and outputs the modes that
31 it extracts from the data. (It preferentially finds modes in the fre‐
32 quency range you specify, but may sometimes find additional modes out‐
33 side of that range.) The data should correspond to equally-spaced time
34 intervals, but there is no constraint on the number of points.
35 Complex numbers in the input should be expressed in the format RE+IMi
37 (no whitespace). Otherwise, whitespace is ignored. Also, comments
38 beginning with "#" and extending to the end of the line are ignored.
39 A typical invocation is something like
41 harminv -t 0.02 1-5 < input.dat
43 which reads a sequence of samples, spaced at 0.02 time intervals (in
45 ms, say, corresponding to 50 kHz), and searches for modes in the fre‐
46 quency range 1-5 kHz. (See below on units.)
47
49 harminv writes six comma-delimited columns to standard output, one line
50 for each mode: frequency, decay constant, Q, amplitude, phase, and
51 error. Each mode corresponds to a function of the form:
52 amplitude * exp[-i (2 pi frequency t - phase) - decay t]
54 Here, i is sqrt(-1), t is the time (see below for units), and the other
56 parameters in the output columns are:
57 frequency
60 The frequency of the mode. If you don't recognize that from the
61 expression above, you should recall Euler's formula: exp(i x) =
62 cos(x) + i sin(x). Note that for complex data, there is a dis‐
63 tinction between positive and negative frequencies.
64 decay constant
66 The exponential decay constant, indicated by decay in the above
67 formula. The inverse of this is often called the "lifetime" of
68 the mode. The "half-life" is ln(2)/decay.
69 Q A conventional, dimensionless expression of the decay lifetime:
71 Q = pi |frequency| / decay. Q, which stands for "quality fac‐
72 tor", is the number of periods for the "energy" in the mode (the
73 squared amplitude) to decay by exp(-2 pi). Equivalently, if you
74 look at the power spectrum (|Fourier transform|^2), 1/Q is the
75 fractional width of the peak at half maximum.
76 amplitude
78 The (real, positive) amplitude of the sinusoids. The amplitude
79 (and phase) information generally seems to be less accurate than
80 the frequency and decay constant.
81 phase The phase shift (in radians) of the sinusoids, as given by the
83 formula above.
84 error A crude estimate of the relative error in the (complex) fre‐
86 quency. This is not really an error bar, however, so you should
87 treat it more as a figure of merit (smaller is better) for each
88 mode.
89
91 Typically, harminv will find a number of spurious solutions in addition
92 to the desired solutions, especially if your data are noisy. Such
93 solutions are characterized by large errors, small amplitudes, and/or
94 small Q (large decay rates / broad linewidths). You can omit these
95 from the output by the error/Q/amplitude screening options defined
96 below.
97 By default, modes with error > 0.1 and Q < 10 are automatically omit‐
99 ted, but it is likely that you will need to set stricter limits.
100
102 The frequency (and decay) values, both input and output, are specified
103 in units of 1/time, where the units of time are determined by the sam‐
104 pling interval dt (the time between consecutive inputs). dt is by
105 default 1, unless you specify it with the -t dt option.
106 In other words, pick some units (e.g. ms in the example above) and use
108 them to express the time step. Then, be consistent and use the inverse
109 of those units (e.g. kHz = 1/ms) for frequency.
110 Note that the frequency is the usual 1/period definition; it is not the
112 angular frequency.
113
115 -h Display help on the command-line options and usage.
116 -V Print the version number and copyright info for harminv.
118 -v Enable verbose output, printed to standard output as comment
120 lines (starting with a "#" character). Also, any "#" comments
121 in the input are echoed to the output.
122 -T Specify period-ranges instead of frequency-ranges on the command
124 line (in units of time corresponding to those specified by -t).
125 The output is still frequency and not period, however.
126 -w Specify angular frequencies instead of frequencies, and output
128 angular frequency instead of frequency. (Angular frequency is
129 frequency multiplied by 2 pi).
130 -n Flip the sign of the frequency (and phase) convention used in
132 harminv. (The sign of the frequency is only important if you
133 have complex-valued input data, in which case the positive and
134 negative frequency amplitudes can differ.)
135 -t dt Specify the sampling interval dt; this determines the units of
137 time used throughout the input and output. Defaults to 1.0.
138 -d d Specify the spectral "density" d to search for modes, where a
140 density of 1 indicates the usual Fourier resolution. That is,
141 the number of basis functions (which sets an upper bound on the
142 number of modes) is given by d times (freq-max - freq-min) times
143 dt times the number of samples in your dataset. A maximum of
144 300 is used, however, to prevent the matrices from getting too
145 big (you can force a larger number with -f, below).
146 Note that the frequency resolution of the outputs is not limited
148 by the spectral density, and can generally be much greater than
149 the Fourier resolution. The density determines how many modes,
150 at most, to search for, and in some sense is the density with
151 which the bandwidth is initially "searched" for modes.
152 The default density is 0.0, which means that the number of basis
154 functions is determined by -f (which defaults to 100). This
155 often corresponds to a much larger density than the usual
156 Fourier resolution, but the resulting singularities in the sys‐
157 tem matrices are automatically removed by harminv.
158 -f nf Specify a lower bound nf on the number of spectral basis func‐
160 tions (defaults to 100), setting a lower bound on the number of
161 modes to search for. This option is often a more convenient way
162 to specify the number of basis functions than the -d option,
163 above, which is why it is the default.
164 -f also allows you to employ more than 300 basis functions, but
166 careful: the computation time scales as O(N nf) + O(nf^3), where
167 N is the number of samples, and very large matrices can also
168 have degraded accuracy.
169 -s sort
171 Specify how the outputs are sorted, where sort is one of fre‐
172 quency/error/Q/decay/amplitude. (Only the first character of
173 sort matters.) All sorts are in ascending order. The default
174 is to sort by frequency.
175 -e err Omit any modes with error (see above) greater than err times the
177 largest error among the computed modes. Defaults to no limit.
178 -E err Omit any modes with error (see above) greater than err.
180 Defaults to 0.1.
181 -F Omit any modes with frequencies outside the specified range.
183 (Such modes are not necessarily spurious, however.)
184 -a amp Omit any modes with amplitude (see above) less than amp times
186 the largest amplitude among the computed modes. Defaults to no
187 limit.
188 -A amp Omit any modes with amplitude (see above) less than amp.
190 Defaults to no limit.
191 -Q q Omit any modes with |Q| (see above) less than q. Defaults to
193 10.
194
196 Send bug reports to S. G. Johnson, stevenj@alum.mit.edu.
197
199 Written by Steven G. Johnson. Copyright (c) 2005 by the Massachusetts
200 Institute of Technology.
201harminv June 4, 2005 HARMINV(1)