System and method for interpretation and visualization of acoustic spectra, particularly to discover the pitch and timbre of musical sounds

[0002] The invention relates to the analysis and understanding of acoustic spectra, particularly the spectra of musical sounds. More specifically, the invention relates to the extraction of pitch and timbre from the spectra of musical sounds.

[0003] Musicians and others (piano tuners, acousticians) are often concerned with evaluating musical sounds in real time, particularly with regard to pitch and tone quality (timbre). During training, rehearsal, and performance, both singers and instrumentalists make these evaluations continuously, and adjust their technique accordingly to improve the sound. Music teachers, orchestral conductors and choral directors make similar evaluations, and by gesture or verbal instruction indicate how performance should be improved.

[0004] In all of these endeavors, the human ear and brain are used to evaluate the sound. Although this mechanism is necessary during performance, and marvelous for judging the “higher” qualities of music such as “expressiveness”, it is hardly ideal for evaluating purely mechanical aspects of sound such as pitch and timbre, because human judgment is subjective. This problem is particularly acute for performing musicians, because the person evaluating the sound is busy producing it. Thus singers and instrumentalists often sing and play off-key while swearing they are in tune, or produce a poor tone quality (timbre) while imagining they are producing a good one.

[0005] The tendency to misjudge can be remedied by training, but in the absence of a teacher, in the hours practicing alone, there is typically no objective measure of pitch and timbre. Several techniques may be used but have limitations and drawbacks. For example, a keyboard instrument may be used intermittently to check pitch, but it dos not give continuous pitch feedback to the musician, and says nothing about timbre. Alternatively, recording and playing back may be used to separate the process of sound production from that of sound evaluation, but this is tedious because it is not real time.

[0006] To solve these problems, a mechanism is needed to provide real-time visual feedback of pitch and timbre to the musician, based on objective and consistent measurements. Visual feedback is ideal because does not interfere with the auditory feedback that the musician must ultimately use in performance. Rather, the visual feedback should help train the auditory system by showing the musician when pitch and tone quality are good. A personal-computer-based software tool would be ideal, since it is flexible, improves automatically as computer technology progresses, and avoids the cost of dedicated instrumentation.

[0007] To analyze sound, particularly musical sound, it is essential to begin, as the ear does, with a spectral analysis. All subsequent analysis, such as the extraction of pitch and timbre, depends on the spectral analysis. Yet, as shown below, it is at this fundamental level of spectral analysis that the prior art is deficient. The prior art's technique for doing spectral analysis is the Discrete Fourier Transform (DFT), and its efficient implementation known as the Fast Fourier Transform (FFT). See Numerical Recipes in C; The Art of Scientific Computing, William H. Press, Brian P. Flannery, et. al., Cambridge University Press, 1988, ISBN 0-521-35465-X, pp. 403-418, which is herein incorporated by reference in its entirety.

[0008] To demonstrate the deficiency of the DFT, it is helpful to summarize some of the mathematics involved. Using the DFT, a signal g(t) (e.g. sound pressure as a function of time) is windowed by a windowing function W(t), such as the Welch Window, which is defined to be non-zero only over the time interval [0, Δt].

[0009] The windowed function

[0010] ĝ(t)≡g(t)W(t)  (1)

[0011] is sampled at N discrete times in the interval [0, Δt], namely $\begin{array}{cc}{t}_n\equiv \frac{n}{S}=\frac{n}{N}\Delta t,n=0,\dots ,N-1,& \left(2\right)\end{array}$

[0012] where S is the sampling rate in Hz. Therefore the total time to measure the N-fold ensemble of samples is $\begin{array}{cc}\Delta t_{\mathrm{meas}}=\frac{N}{S}.\left(\mathrm{Sound}\text{-}\mathrm{measurement}\mathrm{time}\right)& \left(3\right)\end{array}$

[0013] Furthermore, using the DFT, the frequency content of ĝ(t) at frequency f, given by $\begin{array}{cc}\hat{G}\left(f\right)=\frac{1}{S}\sum _{n=0}^{N-1}\hat{g}\left(t_n\right)e^{2\pi if{t}_n},& \left(4\right)\end{array}$

[0014] is evaluated only at certain discrete values of the frequency f namely at the values $\begin{array}{cc}{f}_k\equiv \frac{S}{N}k,k=-\frac{N}{2},\dots ,\frac{N}{2}.& \left(5\right)\end{array}$

[0015] (5)

[0016] Therefore the frequency granularity of the DFT (the difference between two adjacent frequencies, f_k+1−f_k), is $\begin{array}{cc}\Delta f=\frac{S}{N}.\left(\mathrm{Frequency}\mathrm{granularity}\mathrm{for}\mathrm{DFT}\right)& \left(6\right)\end{array}$

[0017] The deficiency of the DFT is summarized by equations (3) and (6), which together imply that

[0018] Δf Δt_meas=1. (Property of Discrete Fourier Transform)  (7)

[0019] That is, with the DFT, it is impossible to achieve both a short sound-measurement time Δt_measand a fine frequency granularity Δf . For example, if it is desired to have a short sound-measurement time of 0.1 seconds, then (7) implies that Δf must assume the rather coarse value of 10 Hz. Conversely, if a fine frequency granularity of 1 Hz is desired, Δt_measmust assume the large value of 1 second. In light of equation (7), it may be concluded that the DFT is inadequate for applications requiring both real-time data acquisition and precise spectral analysis in real time, because such applications require both small Δt_measand small Δf.

[0020] For example, in applications where musical sound need to be measured and analyzed in real time, small Δt_measis necessary to achieve the “real-time” objective. In particular, since fast musical notes are on the order of 80 to 100 milliseconds in length, the application demands

[0021] Δt_meas≦0.1 seconds.  (8)

[0022] For the same type of application, small Δf (frequency granularity) is necessary to achieve accurate results in the computation of pitch. The frequency ratio between two musical notes a half-step apart on the equally tempered scale is $\begin{array}{cc}\frac{f_{+}}{f}=\sqrt[12]{2},& \left(9\right)\end{array}$

[0023] where f₊ is the upper of the two notes and f is the lower of the two notes. Thus the frequency difference between two notes a half-step apart is $\begin{array}{cc}\Delta f_{\mathrm{halfstep}}=f\left(\sqrt[12]{2}-1\right).& \left(10\right)\end{array}$

[0024] For example, at C131 (i.e. 131 Hz, a note in the middle of the range of a human baritone voice), Δf_halfstepis 7.8 Hz. Thus to achieve good pitch resolution of, say, an eighth of a half step, the application demands roughly

[0025] Δf≦1 Hz.  (11)

[0026] Thus, the requirements of such an application with regard to data-acquisition time and frequency granularity, typified by equations (8) and (11), are an order of magnitude more demanding than the capability (7) offered by the DFT. Therefore the DFT is inadequate for such applications, and any prior art that uses it is likewise inadequate. This inadequacy is not dependent on the speed of the computer used to implement the DFT; even if the computer were infinitely fast, the inadequacy would remain the same, because it is inherent in the DFT algorithm itself.

[0027] Because the prior art is thus deficient in its ability to perform real-time data acquisition and finely-resolved spectral analysis simultaneously, it is therefore also deficient in its ability to perform accurate, real-time “note analysis”, wherein the pitch and timbre of the sound are extracted, since note analysis uses the output of spectral analysis as its starting point.

[0028] PC Programs to acquire and analyze sounds using the DFT certainly exist, such as CoolEditby Syntrillium Software and Spectrum Analysisby Sonic Foundary. However, these programs are not typically aimed at real-time applications, and make no attempt to extract pitch and timbre information. As such, they fail to provide useful information to a musician or other user requiring instantaneous, continuous feedback on the pitch and quality of live sound.

[0029] One PC program aimed specifically at musicians is Soloist by Ibis Software. This program provides nothing related to timbre feedback. Moreover it provides only a limited form of pitch feedback; for example, it cannot distinguish notes an octave apart. Furthermore the pitch feedback is not truly “real-time”; only one sound sample is analyzed per metronome beat.

[0030] An object of this invention is a system and method for analyzing the frequency spectrum of a signal in real time, particuarly the spectrum of an acoustic signal having a musical nature, the method providing real-time means to identify the pitch and timbre of the musical note represented by the spectrum, and also providing means to visualize the pitch and timbre, thereby providing real-time visual feedback of musical sounds, particularly to singers and instrumental musicians.

[0031] The invention comprises a transducer, computer hardware, and software. The computer hardware may be a standard, IBM-compatible Personal Computer containing a waveform-input device, such as a Creative Labs' SoundBlaster™ or equivalent. The transducer (such as a microphone) converts a signal (such as sound waves) into a time-varying voltage. The waveform-input device periodically samples this voltage and digitizes each sample, thereby producing an array of N numbers in the memory of the computer that represent a small snippet of the signal measured over a time interval Δt_meas. Snippets are typically measured one after the other at a repetition rate that is inversely related to Δt_meas. The software, also stored in the memory of the computer, and executed using its central processing unit, includes a spectral-analysis process that analyzes the frequency content of each snippet and produces an associated spectrum. The software also includes a novel note-analysis process that analyzes the spectrum and extracts from it the pitch and timbre of the principal musical note contained therein. The process works for any spectrum, including cases where the fundamental frequency of the note is missing. The software further includes novel processes to visualize graphically the pitch and the timbre.

[0032] FIG. 1 is a block diagram of one preferred embodiment of the present invention.

[0033] FIGS. 2A and 2B combined, depict a flow chart of an overall process, including a discrete-Fourier-transform (DFT) process, an alternative logarithmic-frequency-decomposition (LFD) process, a note-analysis process, and various display processes, all executed by the present system.

[0034] FIG. 3 is an example of the output produced by a waveform-display process.

[0035] FIG. 4 is an example of the output produced by the LFD process.

[0036] FIG. 5 is an example of the output produced by a pitch-display process.

[0037] FIG. 6 is another example of the output produced by the pitch-display process.

[0038] FIG. 7 is an example of the output produced by a timbre-display process.

[0039] FIG. 8 is a graph comparing the DFT process to the LFD process with regard to two figures of merit, frequency granularity and process time.

[0040] FIG. 9 is an example of the output produced by the DFT process, showing the DFT's typical, coarse frequency granularity.

[0041] FIG. 10 is an example of the output produced by the LFD process, analogous to FIG. 9, showing the LFD's relatively finer frequency granularity.

[0042] FIGS. 11A and 11B combined, depict a flow chart of the note analysis process.

[0043] In a preferred embodiment, the system described in FIGS. 1 through 10 is used for real-time signal acquisition and spectral analysis. This system is further disclosed and claimed in U.S. Patent Application XXX, entitled “System and Method for Real-Time Signal Acquisition and Spectral Analysis, Particularly for Musical Sounds” to Hall, which is filed on the same day as this disclosure and herein incorporated by reference in its entirety.

[0044] FIG. 1 is a block diagram of one preferred embodiment of the present invention. Sound from a live sound source 105, such as a human voice, musical instrument or other vibrating object, is converted to a time-varying electrical voltage 115 using a microphone 110, such as a Shure BG 4.0, or other appropriate transducer, and this voltage 115 is connected to the computer 120, for example an IBM-Compatible Personal Computer containing a Waveform Input Device 125 such as a Creative Labs' SoundBlaster™ or equivalent. In an alternative embodiment, the voltage 115 is provided by the line output of a tape recorder playing a pre-recorded tape, or by the line output of a compact-disc player playing a pre-recorded disc. In yet another embodiment, the voltage 115 is provided by the line output of an electronic musical instrument, or by the output of an accelerometer attached to a vibrating object. Any voltage representing a vibration is contemplated.

[0045] In addition to the Waveform Input Device 125, the computer 120 comprises a central-processing unit (CPU) 140 such as an Intel Pentium processor, a random-access memory 145, a magnetic disk 150, and a video display 160. The random-access memory stores a software process 200 executed by the CPU; details of this process will be described below. The Waveform Input Device includes an Analog-to-Digital (A/D) Converter 130 that is capable of periodically sampling the time-varying voltage 115 at S samples per second, where S is typically either 8000, 11025, 22050, or 44100. The A/D converter 130 converts each sampled voltage to a signed integer having B bits of precision, and stores it either in the sample memory 135 or in the random-access memory 145. Typically, either B=8 (i.e. each signed integer is in the range −127 to +128), or B=16 (i.e. each signed integer is in the range −32767 to +32768). To obtain data that accurately reflect the sound source, it is important to strive for the highest possible signal-to-noise ratio at the input to the A/D Converter, by means of shielded cables and appropriate amplification if necessary. For example, if the A/D Converter digitizes to 16 bits (i.e. B=16), then a peak ambient-noise level less than ±3000 A/D units is preferred.

[0046] FIG. 2 is a flow chart of the overall software process 200, which is typically stored in the Random-Access Memory 145. The process consists of a number of sub-processes 210 through 270, including sound acquisition 210, analytical processes 230, 235, and 250 that interpret the sound, and display processes 220, 245, 260, and 270 that visualize the results. Following acquisition, analysis and display, the entire process may be repeated, depending on the user-selectable parameter represented by decision box 275. Thus the sound emanating from the sound source 105 may be monitored repeatedly and periodically, the period T being the time necessary to traverse the loop represented by sub-processes 210 through 270. In a typical application, the entire process will be repeated many times, and T will be a fraction of a second.

[0047] The loop begins with Sound Acquisition 210: on the request of a user, the CPU 140 triggers the A/D Converter 130 to begin the acquisition of N samples of the voltage 115, where N is a positive integer and the voltage 115 represents the sound produced by source 105. The values of the integer parameters N, S, and B are selected by commands sent from the software process to the Waveform Input Device via the CPU prior to or simultaneous with the command that triggers the sample acquisition. Software to achieve the sample acquisition is well known in the art. For example, if the computer 120 is an IBM-Compatible Personal Computer (PC), then the software used to achieve sample acquisition typically employs calls to functions in Microsoft's Win32 API. In particular, the API functions waveInOpen( ), waveInPrepareHeader( ), waveInAddBuffer( ), waveInStart( ), and waveInStop( ) may be used.

[0048] When sound acquisition is complete (a process requiring N/S seconds), the raw waveform data (i.e. the integer values logged by the A/D Converter 130 vs. time) may be plotted on the display 160 using the Waveform Display process 220. An example of such a display is shown as FIG. 3; the sound this waveform represents is the note A220 played on the “Oboe” stop of a commercial, electronic keyboard. In FIG. 2, the Waveform Display process 220 is optional, depending on the user-selectable parameter represented by decision box 215. Means for performing such a display process are well known in the art. For example, if the computer is an IBM-Compatible Personal Computer, various functions in a plotting package such as ProEssentials from GigaSoft Inc. or Chart FX from Software FX Inc. may be used.

[0049] Following waveform acquisition, one of two types of spectral-analysis processes is performed on the raw waveform data: either Discrete Fourier-Transform (DFT) analysis 230, or a novel Logarithmic-Frequency Decomposition (LFD) analysis 235. Either type of spectral analysis produces an array of spectral amplitudes Φ(f) that describe the frequency content of the raw data at various frequencies f. The choice of spectral-analysis method depends on the user-selectable parameter represented by the decision box 225. DFT analysis is well known in the art of signal processing. See Numerical Recipes in C; The Art of Scientific Computing, incorporated above. LFD analysis 235 is an alternative to standard DFT analysis that avoids the DFT's shortcomings described earlier.

[0050] When spectral analysis is complete, the spectral data (i.e. the spectral amplitudes Φvs. frequency f) may optionally be plotted on the display 160 using the Spectrum Display process 245. An example of such a display is shown as FIG. 4. This spectrum is the LFD output corresponding to the raw waveform shown in FIG. 3 (i.e. the note A220 played on the “Oboe” stop of a commercial, electronic keyboard).

[0051] In FIG. 2, this Spectrum Display process 245 is optional, depending on the user-selectable parameter represented by decision box 240. Means for performing this display process are well known in the art, as discussed above in connection with the Waveform Display process 220.

[0052] Following spectral analysis, in a preferred embodiment, a novel Note Analysis 250 is performed to extract from the spectrum the pitch and timbre of the musical note contained in the sound. This analysis assumes that the sound contains at most one musical note. If the sound contains more than one, only the most prominent one is extracted. Essentially, note analysis extracts from the sound's spectrum a group of spectral peaks whose peak frequencies are all low integral multiples (to within a tolerance) of a common, fundamental frequency. These frequencies and amplitudes of the spectral peaks are then used to compute the note's pitch and timbre. For example, suppose a spectrum contains spectral peaks at 101, 199, 301, 330, 380, 501, 650, and 702 Hz, and the peaks at 301 and 501 are the largest and second largest in amplitude, respectively. The note-analysis process will determine that the peaks at 101, 199, 301, 501 and 702 Hz are related by being low-integral multiples of a common fundamental at approximately 100 Hz (to within a tolerance). The note analysis therefore identifies these five peaks as belonging to the principle note in the sound, and recognizes them as representing overtones 1, 2, 3, 5, and 7 (“overtone 1” being synonymous with “the fundamental”). The other two peaks, 330 and 380, are rejected as extraneous; perhaps they are attributable to background noise.

[0053] The final part of note analysis is to determine the pitch and timbre of the extracted note. In one preferred embodiment on the invention, the pitch of the note is defined to be an amplitude-weighted average of the estimates that the various overtones make of the fundamental frequency. For the example postulated above, these “estimates of the fundamental frequency” are 101/1, 199/2, 301/3, 501/5, and 702/7, respectively, for overtones 1, 2, 3, 5, and 7. If the amplitudes (in some arbitrary units) of these five note peaks are 10.0, 8.0, 6.0, 4.0, and 2.0 respectively, then the pitch of the note is computed to be $\begin{array}{cc}p=\frac{10.0\left(101/1\right)+8.0\left(199/2\right)+6.0\left(301/3\right)+4.0\left(501/5\right)+2.0\left(702/7\right)}{10.0+8.0+6.0+4.0+2.0}=100.31\mathrm{Hz}.& \left(12\right)\end{array}$

[0054] A useful measure of “timbre” is obtained by defining it to be the vector Q whose i^thelement Q_iis the amplitude of the note's i^thovertone, all amplitudes being normalized by the largest one. Thus, for the above example, the timbre of the note is described by the vector(1.0, 0.8, 0.6, 0, 0.4, 0, 0.2, 0, 0, 0, 0, 0, . . . ), where the non-zero elements are the amplitudes of overtones 1, 2, 3, 5, 7 (normalized by the largest amplitude 10.0), and the zero elements represent the missing overtones in the note. For practical purposes, the length of the vector Q is truncated to some finite value such as 24.

[0055] When note analysis 250 is complete, the computed pitch p may be plotted (vs. the time t that the sound record was acquired) on the display 160 using the Pitch Display process 260. This display process is optional, depending on the user-selectable parameter represented by decision box 255. In reality of course, as described above, the “sound record” comprises N raw samples of the pressure waveform acquired over a time interval N/S, so t is taken to be the center of that interval. Thus, the Pitch Display process 260 involves, for each traversal of the loop 210-275, the plotting of just one point, (p, t), on a graph of pitch (Hz) versus time. In one preferred embodiment of the invention, as sounds are acquired in sequence during multiple traversals of the loop 210-275, the Pitch Display process may be arranged to show the accumulated series of pitches (P1, P2, P3, . . . ) versus times (t₁,t₂,t₃, . . . ). That is, the various points (p_i,t_i) may be plotted simultaneously on a single graph, with new points being added in real time, thereby providing a “live” measure of the pitch that also includes recent history, in the manner of a strip-chart recording. An example of such a display is shown in FIG. 5. This result was obtained on a system for which the loop period T (and therefore the horizontal distance between data points) is about 186 ms when 600 pressure samples comprise each loop's sound record. This figure shows the pitch record for the oboe stop of a commercial electronic keyboard playing the theme from J. S. Bach's “Jesu, Joy of Man's Desiring”. The vertical space on the plot is annotated by lines that document the frequencies of musical “notes”; solid lines represent the keyboard's “white notes”, and dotted lines represent “black notes”. The letter names of the white notes are also given on the right axis. All of the data points in FIG. 5 are exactly on the “note lines”, because the electronic keyboard produces precisely correct pitches. For non-electronic instruments, however, this is rarely the case. FIG. 6, for example, shows the same melody as FIG. 5 performed by a singer instead of a keyboard instrument. The singer's pitch is decidedly less accurate, and the pitch plot quantifies this instant by instant. Such a plot is extremely useful for musicians seeking to monitor pitch in real time.

[0056] Likewise, when note analysis is complete, the timbre Q of the extracted note may be displayed. This display process is optional, depending on the user-selectable parameter represented by decision box 240. In one preferred embodiment of the invention, the elements of the vector Q are displayed in the form of a bar chart. For example, FIG. 7 shows the timbre associated with the spectrum shown in FIG. 4. Thus, the shape of the bar chart illustrates the overtone content (i.e. timbre) of the sound. As sounds are acquired in sequence, the bar chart is updated in real time to reflect the timbre of the most recent sound. Such a plot is extremely useful for musicians seeking to monitor tone quality in real time.

[0057] Two of the analytical sub-processes discussed above, LFD Analysis 235 and Note Analysis 250, form the heart of the invention, and each requires elaboration.

[0058] To remedy the problem discussed earlier under “Problems With the Prior Art”, the set of evaluation frequencies f_knormally used with Fourier analysis, given by eq. (5) above, must be modified, since this choice of frequencies creates the problem expressed by eq. (7). The choice of frequencies f_kin eq. (5) is normally made for two reasons, but it is important to realize that neither applies to the current invention:

[0059] 1. The f_kin eq. (5) produce a set of N numbers Ĝ(f_k) in eq. (4) that can be “inverse transformed” to recover exactly the original N data points ĝ(t_n). Although this is essential in many areas of technology, it is irrelevant for the current invention because, as in the human ear, once conversion to the frequency domain is accomplished, there is no need to retrieve the original pressure signal g(t,). In other words, an invertible “transform” is not needed for this invention.

[0060] 2. The f_kin eq. (5) permit the Discrete Fourier Transform to be efficiently computed using the Fast Fourier Transform (FFT) algorithm, since the latter requires uniform spacing of the evaluation frequencies f_k. For the current invention, however, it is important to execute quickly the entire process loop in FIG. 2, which includes data acquisition as well as computation, so computational efficiency alone is not enough. From this holistic viewpoint, the advantage offered by the FFT's computational efficiency is more than outweighed by the disadvantage of eq. (5), which, as explained above, implies that to use the FFT with good (fine) frequency granularity, data acquisition must be very slow.

[0061] Thus, for purposes of this invention, there is no reason to use the set of evaluation frequencies f_kin eq. (5). In fact, the Fourier-style sums Ĝ may be computed at whatever set of evaluation frequencies makes sense—the f_kmay be distributed as desired, uniformly or not, and there may be any number of them, not necessarily N. For each chosen value of f_k, the square of the complex magnitude of Ĝ will still tell the spectral power density Φ of g(t) at that frequency, namely

[0062] Φ(f_k)≡|Ĝ(f_k)|².  (13)

[0063] Therefore, in accordance with the present invention, Logarithmic Frequency Decomposition (LFD) analysis 235 is provided as a superior alternative to conventional Discrete Fourier Transform (DFT) analysis 230. LFD chooses the following set of evaluation frequencies f_kinstead of eq. (5):

[0064] f_k≡f₀2^k/M, k=0, . . . , K−1.  (14)

[0065] In this equation, f₀is a fixed frequency to be specified, M is the number of evaluation frequencies per octave, and K is an integer to be specified. The idea behind the logarithmic nature of eq. (14) is to distribute the evaluation frequencies uniformly in a musical sense. As is well known in the art of music theory (see R. H. Cannon, Jr., Dynamics of physical Systems, McGraw Hill, 1967, ISBN07-009754-2, herein incorporated by reference), Western music uses a system of equal temperment having 12 “notes” per octave, and the notes have frequencies

[0066] (f_note)_j=f_ref2^j/12, j=0,±1,±2  (15)

[0067] where f_refis a pitch reference such as 440 Hz. The interval between two adjacent notes is called a half step. Thus, eq. (14) chooses Fourier evaluation frequencies f_kthat are uniformly spaced with respect to this musical system. Other equally tempered musical systems (with different numbers of notes per octave) are also accommodated by eq. (14), because the octave (frequency ratio 2:1) is common to all systems. To be specific, however, the following focuses on the 12-tone system given by eq. (15).

[0068] It is useful to choose the free parameters f₀and M in eq. (14) so that one of the evaluation frequencies f_kaligns with each of the musical notes f_note. For this purpose, let f₀be equal to the lowest f_noterequired for a particular application. For example, if sounds from a piano are to be analyzed, choose f₀=27.5 Hz (lowest A on the piano); if sounds from a violin are to be analyzed, choose f₀=196.0 Hz (G below middle C). Since M is the number of evaluation frequencies per octave, it should be chosen as a multiple of 12 to accommodate the 12 half steps per octave,

[0069] M=12 m.  (16)

[0070] Thus m is the number of evaluation frequencies per half step: if m=1, the frequency granularity is a half step; if m=2, the granularity is one half of a half step; if m=5, the granularity is one fifth of a half step, and so on.

[0071] To summarize, LFD analysis 235 has three advantages over traditional DFT analysis 230, for purposes of this invention:

[0072] 1. The LFD's frequency-evaluation points given by eq. (14) are based on uniform frequency ratios, and are therefore distributed evenly in the musical sense. In contrast, the DFT's evaluation points given by eq. (5) are based on uniform frequency differences, and are therefore distributed unevenly in the musical sense, such that the frequency granularity is too coarse at low frequency.

[0073] 2. The LFD's frequency-evaluation points may be located exactly on every musical note, as explained above; it is impossible to achieve this with the DFT.

[0074] 3. As shown in Table 1, the LFD admits five parameters S, N, K, f₀, and m, which may all be selected independently of each other. This leads to flexibility in practice; for example, the minimum frequency may be selected independently of the frequency granularity. Most importantly, the LDF's frequency granularity Δf is independent of the sound acquisition time Δt_meas, since they depend on different independent parameters. In contrast, the DFT admits only two parameters, S and N. This leads to inflexibility; for example, the minimum frequency must be the same as the frequency granularity. Most importantly, having only two parameters leads to the essential dilemma given by eq. (7), ΔfΔt_meas=1, a dilemma that the LFD was designed to solve.

[0075] To obtain the above advantages of the LFD, there is a countervailing disadvantage with regard to computational efficiency. The DFT is typically implemented using the Fast Fourier Transform (FFT) algorithm, whose speed on processors such as a 200 MHz Intel Pentium is so remarkable that the computation time is insignificant compared to the sound acquisition time Δt_meas. Unfortunately, the LFD cannot be implemented by an FFT-like algorithm, because the frequency spacing is logarithmic rather than linear. Nevertheless, reasonable computational efficiency may be obtained as follows.

[0076] The problem is to compute efficiently the Fourier sums $\begin{array}{cc}\hat{G}\left(f_k\right)=\frac{1}{S}\sum _{n=0}^{N-1}\hat{g}\left(t_n\right)e^{2\pi i{f}_k{t}_n},& \left(17\right)\end{array}$

[0077] for the LFD's array of evaluation frequencies

[0078] f_k≡f₀2^k/M, k=0, . . . , K−1.  (18)

[0079] Begin by defining the array of coefficients

[0080] α_k,n≡exp(2πif₀2^k/Mt_n), k=0, . . . , K−1; n=0, . . . , N−1.  (19)

[0081] such that $\begin{array}{cc}\hat{G}\left(f_k\right)=\frac{1}{S}\sum _{n=0}^{N-1}{a}_{k,n}\hat{g}\left(t_n\right),& \left(20\right)\end{array}$

[0082] The key observation, which is clear from eq. (19), is

[0083] α_k+M,n=α²_k,n;  (21)

[0084] This is true because, by adding M to k, a factor of 2 is introduced into the exponential. Physically, eq. (21) says that the coefficient α_k+M,ncorresponding to an evaluation frequency in a given octave is equal to the square of the corresponding coefficient α_k,none octave below. This implies that only the lowest octave's worth of the coefficients need to be computed using trigonometric functions; higher octaves may be computed by recursion. Assuming that C_k,nand S_k,nare the real and imaginary parts of α_k,nrespectively, i.e.

[0085] α_k,n=C_k,n+iS_k,n,  (22)

[0086] then the recursion is as follows:

[0087] α_k+M,n=(C²_k,nS_k,n)+i(2C_k,nS_k,n)  (23)

[0088] Thus, each complex number α_k,ncorresponding to an evaluation frequency in an octave other than the lowest is computed as the square of the corresponding complex number one octave below, thereby avoiding the computationally costly evaluation of trigonometric functions.

[0089] A second observation is that the α_k,ndepend on constant parameters only, not on the data g(t_n). Therefore the lowest octave's α_k,f, which must be computed to seed the recursion, may in fact be computed only once and stored, implying that no trigonometric functions need to be computed during real-time data acquisition and visualization. The amount of storage required for the lowest octave's α_k,nis modest by today's standards: if the number of samples in a sound record is N=1000 and the number of frequency-evaluation points per octave is M=60 (i.e. m=5 divisions per half step), then the storage for one octave of C_k,nand S_k,nis about 469 kilobytes.

[0090] A third observation, assuming that the recursion scheme (23) is used, is that efficient computation of the sums in eq. (20) requires care in the arrangement of loops, to avoid handling the stored numbers C_k,nand S_k,nmore than once. The pseudo-code below shows the proper arrangement of loops. (As given, the pseudo-code assumes that the frequency range covers is an integral number of octaves. If instead the highest octave is incomplete, the pseudo-code must be modified slightly).

[0091] The essential idea in the above pseudo-code is to recall each of the stored complex number a[k,n] (the lowest-octave coefficients that seed the recursion) only once, thereby avoiding memory and cache inefficiencies. Accordingly, the outermost loop on k processes each set of “octave-equivalent” frequency points together, since each such set depends on just one row of the matrix a[k,n]. Likewise, for each octave-equivalent set of frequency points, the middle loop on n processes all terms in the Fourier sums that depend on the same element a[k,n] of the matrix. Finally, the innermost loop on j applies the octave-squaring recursion formula (23). When the loop over the time index n is finished, the Fourier sums for the set k of octave-equivalent frequency points are complete, so the power spectral densities for this set of points may be calculated.

[0092] For example, suppose that the frequency evaluation covers exactly J=4 octaves, starting at middle C(262 Hz), with M=12 frequency-evaluation points per octave—one point located exactly on each true musical note. This implies that there are 48 frequency evaluation points in all. On the first iteration of the outermost loop, with k=0, the above code accesses the real and imaginary parts of each of the numbers a[0,n] only once to compute the power spectral densities for the octave-equivalent frequency evaluation points 0, 11, 23, and 35. These frequencies correspond the four C's (middle C plus 3 higher octaves). Within the k=0 loop, the first iteration of the loop on n handles all terms in the Fourier sums that depend on a[0,0]; the second iteration of the loop on n handles all terms that depend on a[0,1], and so on. On the second iteration of the outermost loop, with k=1, the above code computes the power-spectral densities for frequency evaluation points 1, 12, 24, and 36, which correspond to the four C#'s, and each of the numbers a[ 1,n] is accessed once in the loop on n. Subsequent iterations of the k loop proceed similarly. In this fashion, each complex number a[k,n] is retrieved from memory only once, in memory-contiguous order, to ensure the most efficient use of cache.

[0093] To compare this invention's LFD process to the prior art's DFT/FFT process, it is useful to consider the tradeoff between two figures of merit:

[0094] 1. Frequency granularity normalized in units of musical halfsteps, $\begin{array}{cc}\Delta h\equiv \frac{\Delta f}{\Delta f_{\mathrm{halfstep}}},& \left(24\right)\end{array}$

[0095] where Δf_halfstepis given by eq. (10), and

[0096] 2. Processing time

[0097] Δt_process≡Δt_meas+Δt_spec,  (25)

[0098] where Δt_specis the time required to perform spectral analysis, and Δt_measis given by eq. (3). (Other sub-processes in FIG. 2 are common to the two spectral-analysis methods, so are not included in Δt_process).

[0099] As discussed in the paragraph surrounding eqs. (8) through (11), it is desirable that both of these figures of merit be as small as possible.

[0100] Because the FFT is so efficient, Δt_specfor the DFT/FFT is very small compared to Δt_measand may be neglected, hence

[0101] Δt_process≈Δt_meas·(DFT/FFT only)  (26)

[0102] Recalling eq. (7), substituting eq. (25) and dividing by eq. (10) produces the following equation expressing the DFT/FFT's tradeoff between the two figures of merit: $\begin{array}{cc}\Delta h\approx \frac{1}{f\Delta t_{\mathrm{process}}\left(\sqrt[12]{2}-1\right)}.\left(\mathrm{DFT}/\mathrm{FFT}\mathrm{only}\right)& \left(27\right)\end{array}$

[0103] Thus, for the DFT/FFT, the two figures of merit tradeoff against each other in the form of rectangular hyperbolas, with frequency f as parameter. Four of these hyperbolas, Δh vs. Δt_processfor f=110, 220, 440, and 880 Hz, are plotted as the solid curves in FIG. 8. Each curve is annotated with the value of f and the corresponding musical note name “A”. The small dots on the curves in FIG. 8 represent locations that are typically possible in practice, inasmuch as Δt_meas=N/S, N is an integral power of 2 (FFT requirement), and waveform-input devices typically provide S=8000, 11025, 22050, or 44100 samples/sec.

[0104] In FIG. 8, it is desirable to be close to the origin. In fact, as explained earlier in connection with eqs. (8) through (11), values on the order of 0.1 or less for both Δh and Δt_processare desirable for the musical application addressed by this invention. But the hyperbolas are the best that the DFT/FFT can do to approach the origin. Clearly, this is not good enough for purposes of this invention, particularly at low frequencies such as A110.

[0105] For the LFD, it is impossible to draw generic curves analogous to the DFT's hyperbolas, because Δt_measand Δt_specare both significant for the LFD, and Δt_Specmust be measured experimentally inasmuch as it is dependent on the type of computer 120 and especially upon the type of CPU 135. Thus, the LFD data on FIG. 8 is shown as four sets of discrete data points. For each point, the abscissa is the sum given by eq. (25), where Δt_specwas measured on an IBM Intellisation computer containing a 200 MHz Pentium Pro processor.

[0106] The results show clearly that the LFD vastly outperforms the DFT/FFT in its ability to approach the origin of FIG. 8; that is, to achieve simultaneously fast processing time and fine frequency granularity. In spite of less-efficient spectral computation, the LFD wins because it is not hampered by eq. (7) (the unnormalized form of eq. (27)). This is particularly true at low frequency, such as A110, which is in the middle of the human male's bass voice range. However, it is even true at higher frequency, such as A880, which is at the top of the human female's soprano voice range.

[0107] Notice that each LFD data set approaches a vertical asymptote as the frequency granularity becomes large. This asymptote is Δt_meas; that is, as the frequency granularity becomes large, Δt_specfalls to zero because fewer frequency-evaluation points implies less computation, and so Δt_processapproaches Δt_meas=N/S. Thus, the amount that each LFD data set curls downward to the right is a reflection of spectral computation time Δt_spec. This may be improved in two ways:

[0108] 1. Use a faster computer. Δt_specwill decrease as CPU speed increases.

[0109] 2. Separate data acquisition into a separate computational thread. This allows the spectral analysis of a previously acquired sound to be handled by the CPU 135 while data for the next sound is simultaneously being acquired by the waveform input device 125. This works because data acquisition does not burden the CPU; in fact, in the single-thread implementation, the CPU is essentially idle during Δt_meas.

[0110] In contrast, the curves given for the DFT/FFT on FIG. 8 are not subject to improvements in computer technology. Even if the CPU 135 is infinitely fast, the curves remain as shown, because equation (7) is intrinsic to the DFT algorithm.

[0111] The superiority of the LFD over the DFT/FFT may be further demonstrated by observing the spectra they produce. FIGS. 9 and 10 show a typical comparison, where each frequency-evaluation point in the spectrum is plotted as a data point, and the points are connected by a line for ease of viewing. For this comparison, the sampling rate is S=11025 samples/sec in both cases. The number of samples N (2048 for the DFT/FFT, 1200 for the LFD) and the LFD's frequency granularity (m=10 divisions per halfstep; c.f. eq. (16)) and other parameters (f₀=55 Hz, K=780) have been chosen to produce roughly the same Δt_process(about 185 ms) in the two cases. The sound being analyzed in each case is G below middle C (196.0 Hz) played on the “flute” stop of a commercial keyboard. In both cases, the frequency scale has been expanded to show only the fundamental spectral peak at 196 Hz. The result is clear. On the one hand, in FIG. 9, the DFT/FFT's frequency granularity is very coarse, as shown by the spacing of data points near the horizontal axis. To be sure, this coarse granularity may be reduced by increasing N, but then Δt_meas=N/S, and hence Δt_process, will increase in accordance with eq. (27). The DFT/FFT's coarse frequency granularity in FIG. 9 leads to an artificially truncated spectral peak that fails to find the true peak frequency, because there is no frequency-evaluation point sufficiently near the true peak. On the other hand, in FIG. 10, the LFD's frequency granularity is very fine, allowing the true peak frequency to be accurately found.

[0112] In summary, for purposes of this invention, the LFD is a superior algorithm to the DFT/FFT for the three reasons given earlier in connection with Table 1. Notably, it allows the parameters related to spectral analysis (frequency granularity m and frequency-range parameters f₀and K) to be specified independently of the parameters related to data acquisition (N and S). This superiority has been shown above both in general (FIG. 8) and by specific example (FIGS. 9 and 10).

[0113] Note Analysis

[0114] Note Analysis 250, explained briefly in the foregoing (near eq. (12)), is described in more detail in FIG. 11. The input to Note Analysis is the power spectrum Φ(f_k), a set of real numbers that is easily derived, via eq. (13), from the complex numbers Ĝ(f_k) produced by spectral analysis. If the spectral-analysis method is DFT/FFT, then k=0, . . . , N−1; if the spectral-analysis method is LFD, then k=0, . . . , K−1.

[0115] Note Analysis 250 assumes that the sound emanating from sound source 105 contains at most one “musical note”, where a musical note is defined as a collection of one or more spectral peaks in the power spectrum Φ(f_k) whose peak frequencies are all integral multiples (to within a tolerance) of a common fundamental frequency. An example of such a power spectrum has been given previously as FIG. 4, where all of the peak frequencies, at approximately 220, 440, 660, 880, 1100, 1320, and 1540 Hz, are low integral multiples of the fundamental at 220 Hz. If the sound emanating from sound source 105 contains more than one musical note, only the most prominent one (i.e. the one owning the largest spectral peak) will be found.

[0116] Many vibrating objects—and some musical instruments—do not produce “musical notes” in the sense just defined; that is, they do not produce spectral peaks whose frequencies are integral multiples of a common fundamental. The simplest example is a tuning fork, which is a clamped-free beam satisfying the biharmonic equation. The spectral peak frequencies for such a vibrating object are not integral multiples of each other; for example, the second and third natural frequencies of the tuning fork are 6.267 and 17.55 times the fundamental. Such vibrations are well known in the art of Mechanical Engineering; see, for example, Cyril M. Harris and Charles E. Crede, Shock and Vibration Handbook, 2^ndEdition, pp. 7-11 to 7-15, which is herein incorporated by reference in its entirety. Nevertheless, many important acoustic musical instruments are described by the one-dimensional wave equation to a good approximation. For this equation, the natural frequencies (and therefore the peak frequencies of the musical instrument's power spectrum) are, in fact, integral multiples of a common fundamental. Consequently, Note Analysis 250 deals exclusively with this case. The integral-multiple frequencies are designated “overtones”, and the multiple is designated the “overtone number”. Thus the fundamental itself is designated “overtone 1”, the octave above the fundamental is “overtone 2”, etc.

[0117] It is important to recognize that the fundamental peak may actually be missing in the spectrum; the human ear/brain will nevertheless “hear” a note whose “pitch” is that fundamental frequency. For example, in FIG. 4, if the peak at 220 Hz were absent, the human ear/brain would still hear a “note” with pitch 220 Hz, because the array of overtones imply it. Note Analysis 250 must handle this “missing fundamental” case properly, because it commonly occurs, particularly in the lower ranges of acoustic musical instruments. In general, a “note” at frequency f will be heard if either (1) the fundamental at f is present and zero or more higher overtones are also present, or (2) the fundamental at f is absent and two or more higher overtones of f are present. Note Analysis 250 has been designed to handle both of these cases.

[0118] Refer to FIG. 11 for a detailed description of Note Analysis. In step 405, “peaks” are extracted from the series of numbers Φ(f_k); that is, values of k are sought for which

[0119] Φ(f_k)>Φ(f_k−1) and Φ(f_k)>Φ(f_k−1) and Φ(f_k)>Φ_min.  (28a)

[0120] In other words, a value of Φ is a “peak” if it is bigger than both of its neighbors and also bigger than some user-defined threshold Φ_min.

[0121] In one preferred embodiment of the invention, Φ_minis set to a user-selectable fraction β of the largest value in the array Φ(f_k), i.e.

[0122] Φ_min=βmax{Φ(f_k)}.  (28b)

[0123] The latter condition is useful to filter out noise and other small, insignificant local maxima. Suppose that J values of k, denoted k₀, k₁, k₂, . . . , k_J−1are found to satisfy the conditions in eq. (28). Thus f_kis the frequency of the j^thpeak, and Φ(f_k) is the amplitude of the j^thpeak, where j=0, 1, 2, . . . , J−1.

[0124] In step 410, a data structure P is set up for each of the J peaks found in step 405. Let P[j] denote the data structure for the j^thpeak, and let three members of each data structure be defined:

[0125] P[i].f≡f_k, the frequency of the j^thpeak  (29a

[0126] P[j].h≡Φ(f_k), the amplitude of the j^thpeak  (29b)

[0127] P[j].n, the overtone number of the j^thpeak, to be determined  (29c)

[0128] The primary objective of note analysis is to determine which of the J peaks belong to the principle note contained in the sound, and for those that do belong, to determine the overtone numbers P[j].n.

[0129] In step 415, the peaks are sorted by amplitude so that P[0] refers to the peak with largest amplitude, P[1] refers to the peak with second-largest amplitude, and so on. In FIG. 4 for example, the peaks would be sorted so that P[0] referred to the peak at 1100 Hz, P[1] to the peak at 1320 Hz, P[2] to the peak at 880 Hz, etc.

[0130] In step 420, the overtone number for all peaks except the largest one are initialized to 0 (i.e. P[j].n=0 for j=1, . . . ,J), where 0 indicates that the peak has not yet been found to belong to the note. The overtone number of the largest peak is initialized to 1 (i.e. P[0].n=1), indicating that the largest peak is assigned to the note (n≠0), and is temporarily regarded as the fundamental (n=1). In FIG. 4 for example, the peak at 880 Hz would be assigned to the note and temporarily regarded as the fundamental.

[0131] In step 425, an integer index b is initialized to 1. Throughout the remainder of note analysis, b will represent the index of the peak currently under consideration (i.e. the peak being compared to peak 0) in the loop formed by Steps 430, 435, 440, 445, 450, 455, and 460. The reasons for the various computations in these seven steps, described below, are perhaps difficult to fathom without an example. Therefore, an example is given following the exposition.

[0132] In step 430, a series of tests is made to determine whether the ratio between the frequency of peak 0 and the frequency of peak b is close to the ratio between some pair of low, positive integers i0 and ib. That is, let (i0, ib) range over a two-dimensional array of positive integer pairs, starting with low integers and proceeding to higher ones, each dimension ranging from 1 to n_max. For each pair of integers, determine the truth of the following: $\begin{array}{cc}\left|\frac{P\left[0\right]·f}{P\left[b\right]·f}-\frac{\mathrm{i0}}{\mathrm{ib}}\right|<\delta ,& \left(30\right)\end{array}$

[0133] where δ is a small tolerance. If such a pair of integers (i0, ib) is found to satisfy eq. (30), then the Note Analysis proceeds immediately to Step 435. Otherwise, if the entire two-dimensional array of integers in the range 1 to n_maxis exhausted without satisfying eq. (30), the analysis proceeds to Step 455.

[0134] In practice, n_maxshould be a user-selectable parameter, but by default it may be 24 in a preferred embodiment. Clearly, n_maxcannot be infinite, because then the ratio of any two floating-point numbers in a computer would qualify as the ratio of two integers, since all numbers in a computer are rational, and the test (30) would be meaningless. The idea is to provide only for overtone numbers that are typically seen in practice. Using n_max=24 seems to be a good compromise.

[0135] The tolerance δ in test (30) should account for the fact that the difference on the left-hand side of the inequality may be attributable solely to the finite frequency granularity of the spectral analysis technique (i.e. DFT/FFT or LFD) that produced the power spectrum Φ(f_k). With this in mind, one preferred embodiment of the invention uses the following formulas for δ (wherein P[0].f is abbreviated f₀and P[b].f is abbreviated f_b): $\begin{array}{cc}\delta =\frac{f_0+\frac{\Delta f}{2}}{f_b-\frac{\Delta f}{2}}-\frac{f_0}{f_b}\left(\mathrm{for}\mathrm{DFT}/\mathrm{FFT}\mathrm{analysis}\right)& \left(31a\right)\\ \delta =\frac{f_0{2}^{1/2M}}{f_b{2}^{-1/2M}}-\frac{f_0}{f_b}=\frac{f_0}{f_b}\left(2^{1/M}-1\right)\left(\mathrm{for}\mathrm{LFD}\mathrm{analysis}\right)& \left(31b\right)\end{array}$

[0136] As an example of the above (step 430), consider the spectrum shown in FIG. 4, which was obtained with LFD analysis using M=120 frequency divisions per octave. The actual peak frequencies for the two largest peaks are 1102.36 Hz and 1318.53 Hz respectively. These peaks, or course, represent overtones 5 and 6 of the note A220, but the computer doesn't know this a priori—it must discover the answer based on test (30). Using eq. (31b), test (30) for this case becomes $\begin{array}{cc}\left|\frac{1102.36}{1318.53}-\frac{\mathrm{i0}}{\mathrm{ib}}\right|<\frac{1102.36}{1318.53}\left(2^{1/120}-1\right),\mathrm{which}\mathrm{reduces}\mathrm{to}& \left(32a\right)\\ \left|0.83605-\frac{\mathrm{i0}}{\mathrm{ib}}\right|<0.0048432.& \left(32b\right)\end{array}$

[0137] As the computer applies this test for various values of the integers i0 and ib, it encounters the case i0=5, ib=6, and finds that $\begin{array}{cc}\left|0.83605-\frac{5}{6}\right|<0.0048432& \left(32c\right)\end{array}$

[0138] is in fact satisfied, so the Note Analysis proceeds to step 435 with i0=5, ib=6.

[0139] In step 435, the values of i0 and ib obtained in step 430 are used to compute a proposed new overtone number for peak 0 that accommodates the new peak b, and also to compute a proposed overtone number for peak b itself. These computations are only “proposed” because peak b is not yet accepted into the note. If any of the proposed overtone numbers computed in steps 435 and 440 is too large, peak b will be rejected at step 445, and the proposed new overtone numbers will be discarded. The first computation in step 435 is to compute the proposed new overtone number for peak 0 as the lowest common multiple of the old value P[0].n and the integer i0 found in step 430:

[0140] nNew[0]=LowestCommonMultiple(P[0].n, i0).  (33b)

[0141] Next, compute the proposed new overtone number for peak b:

[0142] nNew[b]=ib*nNew[0]/i0.  (33c)

[0143] In step 440, the value of nNew[0] found in step 435 is used to compute proposed overtone numbers for peaks other than 0 and b (i.e. 0<k<b) that have previously been accepted into the note. Recall that P[k].n=0 means that peak k has not yet been accepted into the note, so the “previously accepted” condition is imposed by looking for peaks having P[k].n>0. Therefore, step 440 performs the following computation:

[0144] for(0<k<b){if(P[k].n>0) nNew[k]=P[k].n*nNew[0]/P[0].n}  (34)

[0145] In step 445, the numbers nNew[ ] found in steps 435 and 440 are checked to see if any is larger than n_max, where n_maxis the user-selectable parameter discussed earlier in connection with eq. (30). That is, the new, smaller-amplitude peak b is rejected if its inclusion in the note would imply that the overtone number n for any of the previously accepted, larger-amplitude peaks would be unreasonably high. If peak b is rejected on this basis, then the Note Analysis proceeds from step 445 to step 455; otherwise, it proceeds to step 450.

[0146] In step, peak b is accepted into the note by copying the proposed overtone numbers nNew[ ] found in steps 435 and 440 into the data structure P[b]. Thus, the following computations are performed:

[0147] P[b].n=nNew[b]  (35a)

[0148] for(0≦k<b){if(P[k].n>0)P[k].n=nNew[k]}.  (35b)

[0149] In step 455, the peak index b is incremented by 1, thereby preparing to consider the next peak in the amplitude-sorted list on the next iteration of the loop containing steps 430 through 460.

[0150] In step 460, the newly incremented value of b is tested to see if the end of the J-element list of peaks P[j] has been reached. If so, computation proceeds to step 465; otherwise, steps 430 through 460 are repeated.

[0151] Before proceeding to a discussion of step 465, consider the following example that illustrates steps 430 through 460 described above. Suppose, on entry to step 420, that a spectrum has been found to contain the following six peaks:

[0152] The task of steps 430 through 460 is to fill in the last column of Table 2 via pair-wise comparison of the peaks. Since this example has been contrived to illustrate the various features of steps 430-460, it should be obvious that all peaks except 4 have been chosen to have frequencies close to multiples of 50 Hz. Therefore, the expected result for the last column of the table is as shown in Table 3, wherein P[4].n=0 implies that peak 4 does not belong to the note. It is expected that all other peaks will be accepted into the note and recognized as the overtones numbers given in the last column.

[0153] In pursuit of the result shown in the last column of Table 3, step 420 begins by artificially filing in the last column as shown in Table 4. In this Table and those that follow, the “Amplitude” column is omitted because, inasmuch as the array of peaks P[ ] is already sorted by amplitude, it is unneeded for the analysis of overtone numbers.

[0154] The initialization of overtone numbers shown in Table 4 implies that the largest peak (j=0) is automatically accepted into the note, and is tentatively considered the fundamental (overtone 1). The other peaks, pending the pair-wise comparisons below, are tentatively considered not to belong to the note, as indicated by the flag n =0.

[0155] On the first iteration of the loop (steps 430 through 460), with b=1, peaks 0 and 1 are compared. Step 430 asks, using eq. (30), if the frequency ratio 603/302 is close to a ratio of two small integers i0/ib within the tolerance given by equations (31). Suppose the tolerance is such that the answer is yes: i0=2, ib=1. Then step 435 computes, via eqs. (33),

[0156] nNew[0]=LowestCommonMultiple(1, 2)=2   (36a)

[0157] nNew[1]=1*2/2=1.  (36b)

[0158] Step 440 requires no computation in this case, because there are no “previously accepted peaks”, so the contents of the “for” loop on eq. (34) is never executed. Since neither value of nNew[ ] in eqs. (36) exceeds n_max=24, peak 1 is accepted into the note (step 445 succeeds). Thus, step 450 is executed, writing the new results into the data structure:

[0159] P[0].n=2; P[1].n=1.  (37)

[0160] In other words, the note is now regarded as containing two peaks, 0 and 1, with overtone numbers as shown in Table 5. This implies that the fundamental frequency of the note is now believed to be roughly 300 Hz.

[0161] On the second iteration of the loop, with b=2, peaks 0 and 2 are compared. Analogous to the first iteration, the ratio test in step 430, using eq. (30), finds 603/198 is approximately equal to 3/1, so i0=3, ib=1. Then step 435 computes, via eqs. (33),

[0162] nNew[0]=LowestCommonMultiple(2, 3)=6  (38a)

[0163] nNew[2]=1*6/3=2.  (38b)

[0164] Step 440 adjusts the overtone number of the previously accepted peak 1 in accordance with eq. (34):

[0165] nNew[1]=1*6/2=3.  (38c)

[0166] Since none of the values of nNew[ ] in eqs. (38) exceeds n_max=24, peak 2 is accepted lo into the note (step 445 succeeds). Thus, step 450 is executed, writing the new results into the data structure:

[0167] P[0].n=6; P[1].n=3; P[2].n=2;.  (39)

[0168] In other words, the note is now regarded as containing three peaks, 0, 1, and 2, with overtone numbers as shown in Table 6. The fundamental is now believed to be approximately 100 Hz, even though no peak has yet been encountered, in the pairwise comparisons, corresponding to that frequency.

[0169] On the third iteration of the loop, with b=3, peaks 0 and 3 are compared. Analogous to the first iteration, the ratio test in step 430, using eq. (30), finds that 603/450 is approximately equal to 4/3, so i0=4, ib=3. Then step 435 computes, via eqs. (33),

[0170] nNew[0]=LowestCommonMultiple(6, 4)=12  (40a)

[0171] nNew[3]=3*12/4=9.  (40b)

[0172] Step 440 adjusts the overtone numbers of the previously accepted peaks 1 and 2 in accordance with eq. (34):

[0173] nNew[1]=3*12/6=6  (40c)

[0174] nNew[2]=2*12/6=4.  (40d)

[0175] Since none of the four values of nNew[ 9 in eqs. (38) exceeds n_max=24, peak 3 is accepted into the note. Thus step 450 is executed , writing the new results into the data structure:

[0176] P[0].n=12; P[1].n=6; P[2].n=4; P[3].n=9 .

[0177] In other words, the note is now regarded as containing four peaks, 0, 1, 2, and 3, with overtone numbers as shown in Table 7. The fundamental is now believed to be approximately 50 Hz.

[0178] On the fourth iteration of the loop, with b=4, peaks 0 and 4 are compared. Analogous to the first iteration, the ratio test in step 430, using eq. (30), finds that 603/320 is approximately equal to 15/8, so i0=15, ib=8.

[0179] Then step 435 computes, via eqs. (33),

[0180] nNew[0]=LowestCommonMultiple(12, 15)=60  (41a)

[0181] nNew[4]=13*60/15=52.  (41b)

[0182] Step 440 adjust the overtone numbers of the previously accepted peaks 1, 2, and 3 in accordance with eq. (34):

[0183] nNew[1]=6*60/12=30  (41c)

[0184] nNew[2]=4*60/12=20  (41d)

[0185] nNew[3]=9*60/12=45.  (41e)

[0186] However, the test in step 445 fails, since nNew[0]=60 exceeds n_max=24. The values of nNew[1], nNew[3], and nNew[4] also exceed n_max; any one of these would cause the test in step 445 to fail. Therefore, peak 4 is rejected as not properly belonging to the note, and the overtone numbers computed in eqs. (41) are discarded. Thus, the table is unchanged except to emphasize that peak 4 has been discarded, as indicated by the bold 0 in Table 8. The fundamental is still believed to be approximately 50 Hz.

[0187] On the fifth and final iteration of the loop, with b=5, peaks 0 and 5 are compared. Analogous to previous iterations, the ratio test in step 430, using eq. (30) finds that 603/99 is approximately equal to 6/1, so i0=6, ib=1. Then step 435 computes

[0188] nNew[0]=LowestCommonMultiple(12, 6)=12  (42a)

[0189] nNew[5]=1*12/6=2,  (42b)

[0190] and step 440 computes

[0191] nNew[1]=6*12/12=6  (42c)

[0192] nNew[2]=4*12/12=4  (42d)

[0193] nNew[3]=9*12/12=9.  (42e)

[0194] Since none of the five values of nNew[ ] in eqs. (42) exceeds n_max=20, peak 5 is accepted into the note. Thus, step 450 is executed, writing the new results into the data structure:

[0195] P[0].n=12; P[1].n=6; P[2].n=4; P[3].n=9; P[5].n=2.

[0196] Actually, nothing changes for the previously accepted peaks; the only new result is P[5].n. The final result is Table 3, as expected. Thus the note finally contains 5 peaks, and the fundamental is believed to be approximately 50 Hz, even though no spectral peak exists at that frequency. Thus, this example illustrates explicitly how the Note Analysis can deal with the case of a missing fundamental, as described earlier just prior to eq. (28).

[0197] In step 465, the pitch and timbre of the note are determined from the frequencies, overtone numbers, and amplitudes of the peaks it comprises, as determined in earlier steps of the analysis.

[0198] Let f_jbe the vector of frequencies P[j].f of the L spectral peaks comprising the note, let n_jbe the corresponding vector of overtone numbers P[j].n found in steps 430 through 460, and let Φ_jbe the corresponding vector of spectral amplitudes P[j].h, where j=0, . . . , L−1.

[0199] For example, using the case given above in Table 3, L=5 (not 6, because peak 4 does not belong to the note), and the vectors are as follows:

[0200] f_j=(603, 302, 198, 450, 99)  (43a)

[0201] n_j=(12, 6, 4, 9, 2).  (43b)

[0202] Φ_j=(64, 56, 41, 29, 16)  (43c)

[0203] The “pitch” of the sound, in units of Hz, may be computed from these three arrays. In one preferred embodiment of the invention, pitch is defined as follows: $\begin{array}{cc}\mathrm{Pitch}\equiv \frac{\sum _{j=0}^{L-1}{\Phi }_j\left(\frac{f_j}{n_j}\right)}{\sum _{j=0}^{L-1}{\Phi }_j}& \left(44\right)\end{array}$

[0204] This definition states that the pitch of a note is the amplitude-weighted average of the estimates, f_j/n_j, that the various spectral peaks make of the fundamental frequency.

[0205] For example, using the case given by eqs. (43), these various estimates of the fundamental frequency are 603/12, 302/6, 198/4, 450/9 and 99/2, each of which is near 50 Hz. Applying (44) to compute a weighted average of these estimates: $\begin{array}{cc}\mathrm{Pitch}=\frac{64\left(\frac{603}{12}\right)+56\left(\frac{302}{6}\right)+41\left(\frac{198}{4}\right)+29\left(\frac{450}{9}\right)+16\left(\frac{99}{2}\right)}{64+56+41+29+16}=50.03\mathrm{Hz}.& \left(45\right)\end{array}$

[0206] As expected, the computed pitch is near 50 Hz.

[0207] In step 465, the timbre of the note is also computed from the arrays n_jand Φ_j. In reality, “timbre” is a complex, difficult-to-define term that includes transient characteristics of a sound such as attack and decay. However, in one preferred embodiment of this invention, timbre for a steady-state note is defined simply as a vector Q≡(Q₁, Q₂, . . . , Qn_Max) of normalized overtone amplitudes, where n_maxis the user-selectable parameter discussed earlier in connection with eq. (30). That is, $\begin{array}{cc}{Q}_{n_j}=\frac{{\Phi }_j}{{\Phi }_{\max }},j=0,\dots ,L-1& \left(46\right)\end{array}$

[0208] where elements of Q not represented in (46) (because the vector n_jonly contains certain integers) are assigned the value 0.

[0209] For example, using the case given by eqs. (43), the non-zero elements of Q in eq. (46) are

[0210] Q₁₂=64/64=1.00,  (47a)

[0211] Q₆=56/64=0.875,  (47b)

[0212] Q₄=41/64=0.641,  (47c)

[0213] Q₉=29/64=0.453,  (47d)

[0214] Q₂=16/64=0.250,  (47e)

[0215] and all other elements are zero. Thus the “timbre” of this note is represented by the vector

[0216] Q=(0.0, 0.250, 0.0, 0.641, 0.0, 0.875, 0.0, 0.0, 0.453, 0.0, 0.0, 1.00, 0.0, 0.0, . . . , 0.0).