# Seismological Research Letters

- © 2013 by the Seismological Society of America

*Online Material: *Interactive numerical model of the Mader–Ott Harmonic Analyzer.

## Introduction

For contemporary seismological observatory practice and seismic exploration work, the fast Fourier transform (FFT) has become such a common tool that we rarely think about single transforms. Often, we switch from time to frequency to time domain without even examining the spectral data. Restituting seismometer response, calculating source spectra, checking frequency‐dependent attenuation, and doing *f*−*k* analysis are only a few applications in seismology, in which the FFT processing steps are an essential tool. Since the advent of modern digital computers, much effort has been made to make the Fourier transform efficient and fast.

The basics for the modern transform algorithms were first described by the French mathematician and physicist (Baron) Jean‐Baptiste‐Joseph Fourier (born 21 March 1768 in Auxerre; died 16 May 1830 in Paris). In 1822, he published his seminal work on heat transport in solid‐state bodies, *The Analytic Theory of Heat*, while at the Ecole Polytechnic Institute in Paris. Fourier not only presented the derivation of the heat transport equations (later called Fourier’s law), but he also proposed a method of resolution, including what we today call a Fourier series. This paper was awarded the prize in mathematics in 1811 by the institute (Hérivel, 1975; Bracewell, 1986).

Spectral analysis has a long history in the geophysical sciences, with work on the earth’s free oscillations as early as the mid‐nineteenth century (e.g., Lamé, 1853). The first recorded seismograph appeared a few decades after Lamé’s work and has been credited to Cecchi in 1887 (Dewey and Byerly, 1969). Time series analysis, the analysis of a sequence of signals characterized by a row vector with (usually) real components (Kanasewich, 1981), did not gain wide theoretical treatment until published works by Wiener (1930, 1949) and Kolmogorov (1939) made calculations with transient signals feasible for geophysical research. Various methods for increased efficiency in the computing of Fourier transforms in the first half of the twentieth century (Cooley *et al.*, 1967) culminated in the first widely disseminated FFT algorithm from Cooley and Tukey (1965). As described by Kanasewich (1981), time series analysis became the basis for significant breakthroughs in mineral exploration with early work, for example, by Wadsworth *et al.* (1953).

Gubbins (2006) described the development of instrumental geophysical investigation as going through four major stages since the late nineteenth century: mechanical, optical, analog magnetic, and digital. Prior to digital computational efforts, scientists working with periodic functions relied heavily on mechanical methods of computation, such as the harmonic analyzer. A nicely comprehensive summary of the evolution of such devices is given by Otnes (2008).

Perhaps the earliest example of a harmonic analyzer is the Kelvin’s tide machine, based on a conception accredited to Lord Kelvin (William Thomson) and constructed in 1872 by the Légé Engineering Company (Parker, 2011). This machine, operated by a system of pulleys and wheels, could be used to compute the first 10 tidal harmonics. As noted by Parker (2011), this device was followed by other designs on both sides of the Atlantic. Parker (2011) gave an excellent and riveting account about the use of these machines, especially during wartime.

Harmonic analyzers were also used to analyze periodic motions of a variety of periodic functions. Two examples include graphical analysis of solar apex and velocity (Nassau and Morse, 1927) and the controversial ether investigations of Miller (1933) as chronicled by Fickinger (2005). Both studies were carried out on the machine designed by Olaus Henrici (and first built by G. Coradi in 1889), which utilized glass spheres rotating brass dials (Henrici, 1894; Miller, 1916). The original Henrici harmonic analyzer is still on display at the Science Museum in London. (The periodical *Philosophical Magazine and Journal of Science*, in which Henrici’s description of his machine appeared, was coproduced by Lord Kelvin.)

Nearly concurrently, as described by Otnes (2008), a much simplified and economical design, based on the implementation of different‐sized cogwheels attached to a planimeter, was published by Mader (1909), based on suggestions from Yule (1894) and implemented with the construction from Ott (1931). It is the operation of this device that will be explored in the present paper.

## Harmonic Analyzer

The Mader–Ott Harmonic Analyzer (MOHA) is a high‐precision mechanical instrument for determining the coefficients of a Fourier series (Fig. 1). Measurements with this instrument are based on a mechanical planimeter device to mechanically determine the surface area of a plane of arbitrary shape. (This is also a useful tool to determine the radius of an equivalent circle for isoseismals in a macroseismic study.) An important aspect of the analog era was the uniqueness of plotted data. Any reproduction or change of scales had to be made by redrawing the plot by hand.

Mader (1909) published the idea of “a simple harmonic analyzer with arbitrary base.” The instrument was designed with several practical advantages: (1) easy handling without the need of reduction or auxiliary calculations; (2) applicable to time series with arbitrary physical length of the baseline (within practical limits), avoiding the need for a transformation of original measurements (always a potential source of error); and (3) economy. Mader (1909) mentioned in a footnote that the instrument was manufactured by Gebrüder Staerzl, Mechanische Werkstätte in Munich and was available at a price of 120 mark. (Following an estimate of Deutsche Bundesbank, this is equivalent to up to 1700 euro.)

The following description of the theory of MOHA is based on the instructions for the use of the instrument given in the German language by Ott (1931). The harmonic analysis of an oscillation requires the decomposition into a fundamental mode and higher modes for a given periodic function *y*=*f*(*x*) with the period *p*. Following Fourier’s 1822 treatise, this decomposition can be achieved by an expansion of trigonometric functions: (1)and the absolute term *a*_{0}/2 corresponds to the average ordinate of *f*(*x*): (2)Fourier’s coefficients *a*_{n} and *b*_{n} are given as (3)

Instead of integrating over the function from 0 to *p*, any other starting point can be chosen as long as the integration is made over a full period. Important for the applicability of MOHA is the simplification of equations (1), (2), and (3) for the case *p*=2*π*: (4)(5)(6)

Superposition of the fundamental harmonic and the higher harmonics with the static component, the latter in seismology commonly called offset, results in function *f*(*x*). If the independent variable is time, as is the case for any seismogram, *x*=*t*, then the period *p* equals the duration of the oscillation *T*, and we get (7)in which *f* is the frequency and *ω* is the angular frequency.

MOHA accomplishes the harmonic analysis following equation (1) by delivering the numerical values of *a*_{n} and *b*_{n} for *n*=1,2,…,*N*, with *N*=25, by following a trace of *y*=*f*(*x*), with a pin as many times as coefficients are to be determined. The instrument is composed of an analyzer guide and a common planimeter (polar or linear). The guide is constructed so that when the pin follows the curve *f*(*x*), the attached planimeter delivers the surface area proportional to the integrals in equation (3), in which *f*(*x*) is multiplied with and . The absolute term can be determined by tracing *f*(*x*) directly with the planimeter, without the use of the analyzer guide. The analyzer guide is based on two principal ideas. (1) A 90°‐angle‐lever *ACB* (Fig. 2), which can rotate around *C* within a limited‐angle *β* due to the two arrestors, projects an arbitrary period *p*=*OO*′ to a fixed period length *l*=*BB*′. The pivot *C* can move along the perpendicular bisector *CD* of the period axis *OO*′. The moving tracer pin *P* can be adjusted in its position along the instrument arm *AC* in such a way that when the arm *CB* reaches its end positions at the arrestors, *P* is at the beginning and end of the period at *O* and *O*′, respectively. Figure 2b shows that movement of the pin *P* along the *x* axis provokes the proportional movement of *B* along the *y* axis by the amount *z*=−(*l*/*p*)*x*.

(2) A cogwheel rolls along a cograil of length *l* exact *n* times. For displacement of the cograil of distance *z*, the cogwheel will rotate with the proportional‐angle *ω*=−*n*(2*π*/*l*)*z*. For better clarity in Figure 3, the cograil is fixed, and the cogwheel is free to move. The angle lever, cograil, and cogwheel are arranged on a carriage whose movement is limited to the *y* direction. During operation of the instrument, the pivot of the angle lever *C* must move along the perpendicular bisector *CD* of the period axis *OO*′. The end point *B* of the angle lever moves the cograil by the amount *l* when the moving pin *P* moves along the curve of *y*=*f*(*x*) from start point *O* to end point *O*′, independent of the length of *p*=*OO*′. This is technically achieved by a parallel sliding block. For movement of *P* by the amount *x* in the *x* direction, the cograil moves *z*=−(*l*/*p*)*x*, which is equal to a cogwheel rotation of (8)

The cogwheel of radius *R*(*l*=2*πnR*) has two small boreholes positioned at a 90° angle that are marked with the letters *c* and *s* indicating the positions for the attachment of the planimeter pin to determine the cosine and sine coefficients, respectively.

To understand the movements of the planimeter attached to the boreholes *c* and *s* of the cogwheel, a new coordinate system with the axes *ξ* and *η* (Fig. 3b) is introduced. At the beginning of a measurement, the cogwheel is positioned in such a way that the *s* boring lays on the negative *ξ* half axis at *ξ*=−*R* and the *c* boring lays on the positive *η* axis at *η*=*R*. For a fixed carriage and guiding the moving pin *P* from *O* to *O*′, *P* will move along a circular arc with the ordinates (Fig. 4). If *P* is steered along the function *y*=*f*(*x*), the movement of the carriage in the *y* direction is . Because of this displacement and because of the rotation of the cogwheel of *ω*=*n*(2*π*/*p*)*x* due to the displacement *x* in *x* direction, the coordinates of the *s* boring are changed (Fig. 3): (9)

If we move the pin *P* along *f*(*x*) from *O* to *O*′ and afterward back along the straight axis from *O*′ to *O*, we have *f*(*x*)=0. The *s* boring with respect to the stationary *ξη* coordinate system describes a closed curve with surface area: (10)in which (11)

Separating forward and backward paths, we can write (12)

The entire second integral cancels out with the corresponding parts of the first integral from the movement of the tracer pin back from *O*′ to *O*, which is independent of the forward path, and the remaining part of equation (12) becomes (13)

The term *πRn* is the device constant *K*. Consequently, the planimeter with its moving pin in the *s* boring of the cogwheel shows the surface area, (14)which is just a constant factor times the Fourier coefficient *b*_{n}. When the planimeter pin is placed in the *c* boring, the corresponding coefficient *a*_{n} will result.

The MOHA used in this study has a *K*=100 mm; a planimeter value of 0.1 cm^{2} corresponds to a Fourier coefficient of 0.01 cm. The ingenuity of the MOHA construction is the projection of all period lengths *p* to a fixed lateral movement of the cograil of *l*, which makes the resulting Fourier coefficients *a*_{n} and *b*_{n} independent of *p*. Our model of MOHA allows period lengths between 2.5 and 36.0 cm and is equipped with 22 cogwheels for the first 25 Fourier coefficients.

For coefficients higher than 6, a conversion cogwheel is used in addition to the specific cogwheel of the particular coefficient. A unique feature of the instrument is the option of installing two parallel coefficient cogwheels. In other words, with two planimeters attached to the two wheels (Fig. 1), two coefficients can be determined with one run of the pin along the curve, cutting the total measuring time in half.

Theoretically, by stepwise changing the period length *p* of a record, hidden periodicity could be uncovered. Furthermore, by partitioning a record in *u* parts of equal length, analyzing each part separately, summing up *a*_{n} and *b*_{n}, and dividing by *u*, the coefficients for the *u*·*n*th harmonic of *f*(*x*) can be determined.

## Practice

To gain experience with MOHA, we determined coefficients of an analytic signal. The function shown in Figure 5 is the superposition of seven harmonic functions, including the fundamental harmonic and those of the order 2, 3, 8, and 9, with amplitudes between 1.0 and 3.0. Students and staff of the seismological station Bensberg were challenged by one of the authors to determine the first 10 cosine and sine coefficients of the analytic function using the MOHA device. Figure 6 shows the result of the seven test runs. With the exception of the coefficient 9 s, the measurements not considered as outliers are all better than 8% error. This is compatible with earlier tests made by Korschunow (1955). Mader’s (1909) own test with a sawtooth function resulted in errors between 0% and 9% for the first seven coefficients.

A test with a more complicated function was made using a seismogram from a local mine tremor with *M*_{L} 4.1 recorded at 154‐km distance. Because of the shallow source (<2 km), the mine tremor caused distinct surface waves. The restituted displacement of the vertical component is shown in Figure 7. The original digital record was sampled at 125 Hz. For comparison, a time window containing 4096 samples (Fig. 7) was transformed with a standard FFT. The normalized spectral amplitudes of the first 25 harmonics and the offset at the zero frequency are shown in the bottom of Figure 7. The time window was plotted for the analysis with a horizontal scale of 36 cm, the maximum the MOHA can handle. The first 25 Fourier spectral amplitudes are also shown in Figure 7. The value for the zero frequency was determined with the planimeter without the use of MOHA by measuring separately the area under the positive and negative amplitudes. Coefficients 20, 22, and 24 were determined in sections by dividing the seismogram into two parts of equal length and using cogwheels 10, 11, and 12, respectively. The averages of the two cosine and sine coefficients of the two sections are calculated, which is a necessity as no cogwheels for these coefficients exist. However, as stated earlier, this method also enabled derivation of higher coefficients than the first 25.

## mathMOHA

As a demonstration tool, but also to pay tribute to the engineering ingenuity of Mader and Ott, we generated mathMOHA (Fig. 8), a fully interactive numerical model of MOHA in an open CDF format (http://www.wolfram.com/cdf/, last accessed January 2013). It is attached to this paper as an electronic supplement and can be run with the freely available CDF player (http://www.wolfram.com/cdf-player/, last accessed January 2013). mathMOHA illustrates the determination of the Fourier series coefficients *a*_{1} to *a*_{10} and *b*_{1} to *b*_{10} for a set of selected time series, which includes the ones discussed in this paper. mathMOHA can operate in two modes of operation.

In the Freehand mode, the users control the trace point *P* (see Fig. 8) with their mouse. Once the end of the trace is reached, the trace point can be automatically returned to the origin by first moving vertically to the zero line and then back to the origin by hitting the Automatic Finish button (see Fig. 8). This way, the tracking error is due purely to the imprecision of the user in following the actual trace and not from closing the loop to the origin.

In the Autotracking mode, the actual tracking error is further reduced to the numerical and sampling error because the user moves the trace point with a slider along the time series and back to the origin. This mode is particularly instructive for the illustration of the working principle of MOHA. Simultaneous with the movement of the trace point, the simulated planimeter (the framed yellow rectangle attached to the planimeter pivot point in Fig. 8), which is a simple implementation of equation (10), displays the momentary area within the trace of the selected cogwheel pin. The cogwheel pins and the corresponding pin traces are colored red and blue for the *b* and *a* coefficients, respectively. Figure 8 displays the virtual model corresponding to the actual setup as shown in Figure 1a.

## Conclusion

With the reactivation of the mechanical harmonic analyzer, we have been able to reproduce earlier results (Mader, 1909; Korschunow, 1955) that showed that at least the lowest 25 Fourier coefficients can be determined with the analog MOHA, usually with errors smaller than 8%. Considering the fact that the analysis of the first 10 coefficients of our test signal took on average about 2 hr, the title of this contribution seems justified. We suggest that students in the course of learning about signal processing should, as an exercise, perform a slow Fourier transform as a way to deeply understand that Fourier integrals are defined by the area under the curve to be analyzed. As original MOHAs have become rare instruments, the electronic version serves as a viable alternative. FFT routines can easily be added to analysis codes and are indispensable to modern signal processing. However, in many ways similar to a finely crafted instrument in the hands of a well‐trained musician, there is something quite special and intrinsically satisfying about the feeling of working with a masterpiece of craftsmanship such as MOHA.

Harmonic analysis of signals is of interest in many disciplines in addition to seismology, and as alluded to in the Introduction, it is of special interest in gravity and tidal studies. Parker (2011) reported that harmonic analysis and synthesis of tides were key features for the success of the D‐day invasion during World War II. Spies behind the German lines had carefully studied the local tides along the coast of Normandy. These empirical time series were analyzed, most probably with an instrument such as the one presented in this contribution. The resulting Fourier coefficients were the key for the precise prediction of the height and timing of the tide at beaches Omaha, Gold, and Juno. For the synthesis of those data, mechanical computing devices that are the size of a closet and much more impressive than MOHA were used.

## Acknowledgments

We thank S. Falter, C. Fleischer, S. Gerz, H. Kehmeier, G. Schweppe, and J. Tzislakis for taking their time to make a slow Fourier transform. Comments from the Editor, J. Ebel, helped to significantly improve the manuscript. Most of all, we thank Ali Akasya for making the instrument available that was used in our study.