WO1994024662A1

WO1994024662A1 - Method of calculating filter weights for compression wave cancellation systems

Info

Publication number: WO1994024662A1
Application number: PCT/US1994/004368
Authority: WO
Inventors: J. Paul Martinez; Donald W. Cooley; Kenneth G. Oetzel; Carla P. Woodworth; V. Bradford Mason; Julia A. Oekin-Meza; Anne R. Hessing; D. Kent Peterson; Richard H. Steinberger
Original assignee: Sri International
Priority date: 1993-04-21
Filing date: 1994-04-20
Publication date: 1994-10-27

Abstract

A method of calculating filter weights for a multichannel compression wave cancellation system (10) is characterized by the following steps: system configuration (48), calibration data collection (50), calibration data preprocessing (52), filter weight calculation (54), system performance simulation (56), and filter weight downloading (58). The filter weight calculation is preferably both for forward filters (F) and neutralization filters (N). The system configuration step specifies the numbers and locations of the input transducers (18) and output transducers (22). The data collection step determines the responses of all input transducers and quiet zone transducers (24) to a noise source (16) and selected output transducers. The preprocessing step windows, averages, deconvolves, and filters the calibration data in preparation for filter weight calculation. The forward and neutralization filter weights are derived from the preprocessed data and may be constrained to limit power or energy outputs. The forward filter weights for small multichannel systems may be calculated by a direct method, but for large multichannel systems an iterative method is preferred. Finally, system performance is simulated with the calculated forward and neutralization filter weights to determine if the filter weights need further refinement before they are downloaded into a hardware signal processing unit.

Description

METHOD OF CALCULATING FILTER WEIGHTS FOR COMPRESSION WAVE CANCELLATION SYSTEMS

Description

Technical Field

This invention relates generally to sound dampening techniques and more particularly to methods and apparatus for active noise cancellation.

Background Art It is often desirable to reduce the ambient noise level in a particular environment. This is particularly true when the noise is loud and unpleasant, such as the noise produced by machinery. In fact, loud noise can be more than annoying: at certain sustained levels it can cause pain and permanent hearing loss.

One approach to noise reduction is active noise control, which refers to any electro-acoustical method in which an undesired sound wave is canceled by a second sound wave that has the same amplitude but is 180° out of phase. In practice, active noise control is most often used to attenuate low frequency noise and vibration.

There are a great number of patent disclosures describing active noise control systems. Examples of patents describing active methodologies for quieting heating, ventilation, and air conditioning (HVAC) installations include: U.S. patents 4,122,303; 4,171,465; 4,473,906; 4,480,333; 4,596,033;

4,665,549; 4,669,122; 4,677,676; 4,677,677; 4,783,817; 4,815,139; and 4,837,834. Some of these patents, such as patents 4,473,906 and 4,665,549, disclose the use of multiple input microphones to detect the noise to be canceled. Others of these patents, such as patents 4,171,465 and 4,669,122, disclose the use of multiple loudspeakers to cancel noise in a duct. U.S. patent

4,815,139 discloses both multiple input microphones to sense noise and multiple loudspeakers to cancel noise in a duct. Other examples include U.S. patent 4,637,048, which teaches the cancellation of noise from an automobile tail pipe; and U.S. patents 4,562,589, 4,689,821 and 4,715,559, which teach the cancellation of noise in the fuselage or cockpit of aircraft. All of the above-identified patents operate on noise within enclosed spaces such as ducts or airplane fuselages, and they all disclose single channel signal processors. As such, it is relatively straightforward to determine the correct filter weights for these prior art noise cancellation systems. Typically, microphones sample the noise at one or more locations within the enclosed space, and a simple set of simultaneous linear equations is solved to calculate the proper filter weights. These filter weight calculations can often be done by hand or with a simple calculator.

U.S. patent application 07/697,154 filed on May 8, 1991, on behalf of

Martinez et al, assigned to the assignee of this present application, and entitled "Method and Apparatus for the Active Reduction of Compression Waves" describes a multichannel active noise control system which includes a plurality of input microphones and at least one (but preferably more) output loudspeakers. The number of channels defined by the system is equal to the product of the number of input microphones and the number of output loudspeakers. This multichannel noise control system is very advantageous in that it provides three-dimensional noise cancellation in unbounded environments.

If the number of channels in a Martinez et al. system is small (e.g., less than 12), the calculation of filter weights is relatively straightforward, as was the case with prior art noise cancellation systems. For example, the Wiener least-squares minimization methods as taught in Optimum Signal Processing. An Introduction, by S. J. Orfanidis, Macmillan Publishing Company, 1988, or the Widrow-Hoff algorithm as taught in Adaptive Signal Processing, by B. Widrow and S. Stearns, Prentice-Hall, Inc., 1985, are suitable for calculating filter weights if the number of channels is small. However, as the number of channels increases, traditional methods for calculating filter weights become cumbersome because thousands of linear equations must be solved simultaneously. Using standard techniques, these computations can take days or even weeks. It would therefore be desirable to have a faster method for calculating filter weights for large multichannel noise control systems of the kind described in the Martinez et al. patent application.

Because it was relatively simple to calculate filter weights in prior art noise control systems, it typically was not necessary to worry about real world constraints on the system. For example, if certain filter weights would create a peak output power that might destroy the output loudspeakers, new filter weights could be easily calculated by changing some of the input parameters. This approach is not feasible for large multichannel systems. It is therefore also desirable to have a method for including system constraints into the filter weight calculations to eliminate impractical solutions.

It is fairly easy to optimize the performance of a single channel system of the prior art. Optimal locations of the input microphone(s) and output loudspeaker(s) can be determined experimentally by successive system calibrations without changing the signal processor. In a large multichannel system, however, the signal processor is a major variable that affects not only system performance but also hardware complexity and cost. In general, it is neither technically nor economically feasible to modify the design of a multichannel signal processor for experimental purposes. Therefore, it is desirable to have a method for simulating alternative processor designs and evaluating performance before a system is built.

Disclosure of the Invention

The present invention comprises a method of calculating filter weights for compression wave cancellation systems such as multichannel noise control systems. The method includes the steps of configuring a system, collecting calibration data, preprocessing the data, calculating the forward filter weights W and the neutralization filter weights Wc, simulating the system to predict performance, and loading a set of filter weights (W and Wc) into the signal processor of a multichannel noise control system.

In the configuration step, the number and placement of a plurality of input transducers (microphones) and output transducers (loudspeakers) is determined. The placement of input transducers depends on the directional characteristics of the noise source; the placement of the output transducers governs the size and shape of the quiet zone.

The calibration data collection step determines the transfer functions of four parameters which define the multichannel system. Input microphone calibration data M is the response of the input microphones to typically wideband signals emitted from a noise source. Quiet zone calibration data Q is the response of the quiet zone microphones to the noise source. Loudspeaker calibration data L is the response of the quiet zone microphones to signals emitted from the loudspeakers. Coupling calibration data C is the response of the input microphones to signals emitted from the loudspeakers.

The preprocessing step windows, averages, deconvolves, and filters the calibration data in preparation for filter weight calculations. The filter weight calculation step can optionally use power or energy constraints on the signals feeding the output transducers; it also offers a choice of a direct or an iterative method. In the present implementation of the invention, the direct method for calculating filter weights is used for systems with less than twelve channels, and the iterative method is used for systems with twelve or more channels.

The final step is a simulation of system performance based on the calculated forward filter weights W and the neutralization filter weights Wc. If system performance is determined to be inadequate, the configuration of the system can be modified, or the parameters used in the calculation of the forward filter weights W can be tweaked to improve system performance. The preprocessing and compute filter weights steps can then be repeated to produce a new set of filter weights. Once the filter weights in the simulated system appear to meet system requirements, they are loaded into the signal processor of the multichannel noise control system.

The principal advantage of the present invention is that it provides a means for calibrating large, multichannel active noise control systems that are computationally more demanding but capable of greater performance than single channel systems of the prior art. Furthermore, for large systems, the iterative method of the present invention for calculating forward filter weights is much faster than the direct method, greatly reducing computer processing time. System constraints can be taken into account when calculating filter weights, so that optimum performance can be realized with imperfect transducers. Finally, the invention facilitates system design in that different signal processors can be simulated and their performance compared before a design is committed to hardware.

Brief Description of the Drawings Figure 1 is a pictorial view of a compression wave cancellation system;

Figure 2 is a block diagram representing the compression wave cancellation system of Fig. 1; Figure 3 is a block diagram of the signal processing unit (SPU) 20 of Figs. 1 and 2;

Figure 4 is a flow diagram of a method for calculating filter weights in accordance with the present invention;

Figure 5 is an example of a noise control system produced by the Configure System step 48 of Fig. 4;

Figures 6a and 6b are graphs used to illustrate the Calibration Data Collection step 50 of Fig.4;

Figure 7 is a flow diagram of a process for calculating the forward filter weights W and neutralization filter weights Wc of the Compute W and Wc step

54 of Fig. 4;

Figure 8 is a flow diagram of a direct method for calculating forward filter weights W within the process illustrated in Fig.7;

Figure 9 is a flow diagram of an iterative method for calculating forward filter weights W within the process illustrated in Fig. 7; and

Figure 10 is a block diagram of a system simulation used to simulate performance in the Simulate Performance step 56 of Fig.4.

Best Modes for Carrying out the Invention

In Fig. 1, a multichannel apparatus 10 is used to reduce undesired compression waves 12 in a medium 14 produced by a source 16. A full description of apparatus 10 can be found in U.S. patent application serial number 07/697,154 filed on May 8, 1991, on behalf of Martinez et al., assigned to the assignee of this application and entitled "Method and Apparatus for the Active Reduction of Compression Waves", said application being incorporated herein by reference.

Briefly, apparatus 10 includes a number of input transducers such as microphones 18a and 18b (collectively input transducers 18), a signal processing unit (SPU) 20, and a number of output transducers such as loudspeakers 22a, 22b, 22c, and 22d (collectively output transducers 22). As used herein, the term "loudspeaker" means any electro-acoustical transducer capable of transforming electrical signals into compression waves in a medium. A quiet zone transducer such as microphone 24 can be used to detect the effectiveness of apparatus 10 in reducing undesired compression waves in a quiet zone 26 of medium 14. Microphone 24 can be moved to a number of positions 24' to evaluate the effectiveness of apparatus 10 at various angular positions relative to source 16. Alternatively, a number of quiet zone microphones can be used to simultaneously measure the residual compression waves in the quiet zone.

Fig. 2 illustrates the system of Fig. 1 in block diagram form. Medium 14 in this example is air, and the propagation paths through medium 14 are represented by the transmission path boxes 28, 30, 32, and 34. . Sounds traveling along transmission paths 28 (from source 16) and transmission paths

34 (from output transducers 22) are combined in medium 14 as indicated by summer 36 and picked up by input transducers 18. Signals 38 from input transducers 18 are presented to SPU 20, and signals 40 from SPU 20 are applied to output transducers 22. Sounds traveling along transmission paths 30 (from source 16) and transmission paths 32 (from output transducers 22) are combined in medium 14 as indicated by summer 41 and picked up by quiet zone transducer 24.

Input transducers 18 convert the combined compression waves from summer 36 to electrical signals 38. SPU 20 processes electrical signals 38 to produce complementary electrical signals 40. Output transducers 22 produce complementary compression waves in medium 14 along transmission paths 32 to at least partially cancel the compression waves created by source 16 along transmission paths 30.

Referring now to Fig. 3, SPU 20 includes a pair of input summers 42a and 42b, eight forward filters F, four output summers 44a, 44b, 44c, and 44d, and eight neutralization filters N. The two digit subscripts of the forward filters F identify the inputs and outputs they couple together. For example, forward filter F_{j j} couples input 1 to output 1, and forward filter F₂₃ couples input 2 to output 3. In a similar fashion, the eight neutralization filters have two digit subscripts where the first digit indicates which input the filter is coupled to and the second digit which output it is coupled to. SPU 20 further includes a pair of input buffers 43a and 43b coupled to the two input microphones 18a and 18b, respectively, and four output buffers 45a, 45b, 45c, and 45d coupled to the loudspeakers 22a, 22b, 22c, and 22d, respectively.

In the forward path of SPU 20, inputs 1 and 2 are presented to summers 42a and 42b, respectively, and the output of each summer is applied to four forward filters F. The signal from summer 42a on a line 46a is applied to forward filters F, _j, F₁₂, F₁₃, and F₁₄. Similarly, the signal from summer 42b on a line 46b is applied to forward filters F₂ι , F₂₂, F₂₃, and F₂₄. The outputs of the forward filters F are applied to summers 44a-d in the following fashion: the signals from filters F_{j j} and F₂₁ are applied to summer 44a, the signals from filters F,₂ and F₂₂ are applied to summer 44b, the signals from filters F_j and

F₂₃ are applied to summer 44c, and the signals from filters F₁₄ and F₂₄ are applied to summer 44d. The outputs of summers 44a-d are coupled to buffers 45a-d, respectively.

In the reverse path, output signals 1-4 are fed back through neutralization filters N to summers 42a and 42b. More specifically, neutralization filters N_{j j}, N₁ , N₁₃, and N₁₄ feed back the signals from outputs 1-4 to summer 42a, and neutralization filters N₂₁, N₂₂, N₂₃, and N₂₄ feed back the signals from outputs 1-4 to summer 42b.

The filters F and N can be made from discrete components such as inductors, capacitors, and resistors. Preferably, however, the filters F and N are digital filters and part of a digital SPU 20. In a digital SPU 20, input buffers 43a and 43b include analog-to-digital (A D) converters which convert the analog signals from microphones 18a and 18b into digital inputs 1 and 2, respectively. As is well known to those skilled in the art, the input buffers can also include pre-amplifiers, anti-aliasing (low-pass) filters, etc. Lines 46a and 46b carry the digital sum calculated by digital summers 42a and 42b to the digital forward filters F. The outputs of the digital forward filters are applied to digital summers 44a-d to produce digital outputs 1-4. Output buffers 45a-d include digital-to-analog (D/A) converters to convert the digital outputs 1-4 to analog signals capable of driving loudspeakers 22a - 22d. As is also well known to those skilled in the art, the output buffers can include reconstruction filters, power amplifiers, etc. The digital signals on outputs 1-4 are fed back through digital neutralization filters N to digital summers 42a and 42b.

For small multichannel systems (e.g., less than 12 channels), the weights of the forward filters F and the neutralization filters N can be calculated by the conventional techniques described previously. However, as the number of channels of SPU 20 increases, it becomes increasingly difficult to calculate the filter weights in a reasonable period of time. The present invention comprises methods for quickly and efficiently generating filter weights for large multichannel noise control systems.

As is well known to those skilled in the art, it is convenient to characterize a linear filter in the time domain as a multitap delay line, where the taps are coupled to a summer by a plurality of adjustable gain amplifiers. Any linear filter can be created by properly adjusting the length of the delay line and the gains of the amplifiers. This approach is particularly attractive for digital filters, because digital signal processors are well suited for implementing multitap delay lines. A textbook on digital filters is Discrete-Time Signal Processing written by Alan Oppenheim and Ronald Schafer, published by Prentice Hall (1989). In the following discussions, the characteristics of a filter will be described in terms of the length and tap gains of the respective delay line.

Figure 4 is a flow diagram illustrating a process for generating filter weights in accordance with the present invention. In step 48, the system is configured by positioning a number of input microphones and loudspeakers with respect to a noise source and a designated quiet zone (see Fig. 5). In step 50, calibration data are collected. The calibration data are then preprocessed in step 52, where the signals are windowed, averaged, deconvolved, and filtered in preparation for filter weight calculations. Forward and neutralization filter weights are computed in step 54, and system performance is simulated in step 56. If the simulated performance is substandard, steps 52, 54, and 56 are repeated, as indicated by broken line 57. In step 58, the values of W and Wc are downloaded into the SPU. When the calibrated system is activated, the antinoise substantially cancels the undesired noise in the quiet zone.

Figure 5 is a diagram showing the physical configuration of an active noise control system. The numbers and locations of the input microphones, loudspeakers, and quiet zone microphones are determined by the specific noise problem to be solved. To illustrate how system components are typically arranged in space, a system 10' is shown as consisting of a signal processing unit (SPU) 20', two input microphones Mi and M2, four loudspeakers Li, L2, L3, L4, and two quiet zone microphones Qi and Q2. A complex, wideband _, noise source S produces compression waves which are to be canceled by complementary compression waves produced by system 10'. Quiet zone 66 is the spatial region where noise reduction is to be achieved; the boundary of the quiet zone is indicated by a broken line 68. In SPU 20', every input microphone (Mj and M2) is preferably connected to every loudspeaker (Lι_t L2, L₃, L4) through a forward filter. The forward filter weights W will be chosen to minimize the residual noise in the quiet zone.

Figures 6a and 6b illustrate a preferred method for collecting calibration data in accordance with the present invention. To produce a set of filter weights, input microphone calibration data M, quiet zone calibration data Q, loudspeaker calibration data L, and coupling calibration data C are collected. SPU 20' is bypassed during the data collection process. To collect the M and Q data, a "Condition A" is used as illustrated in Fig. 6a, where the noise source S is active over its entire operating range and the loudspeakers are turned off. To collect the L and C data, a "Condition B" is used as illustrated in Fig. 6b, where the noise source S is turned off and the loudspeakers are energized, one at a time, with a probe signal p.

With reference to Fig. 6a, while noise source S is active, data are simultaneously collected at all input and quiet zone microphones; these data constitute the time series M]S, M2S, QjS, and Q2S. Each of these time series simply represents the response of a specific microphone to noise source S. In theory, the time series should be long enough to encompass the entire repertoire of source characteristics, including maximum amplitude, frequency range, and directionality. But because storage space on computers is limited, it is often necessary to select relatively brief "time slots" that contain only the most important source characteristics. In Fig. 6a, three time slots Si, S2, and S₃ are indicated by broken lines. Each of the time series is "windowed" at the selected time slots, and the resulting data are saved as a set of windowed time series. Thus, the M data now comprise the time series M1S1. M1S2, MiSβ, M2S1,

M2S2, and M2S3. Similarly, the Q data now comprise the time series Q1S1_. Q1S2, Q1S3, Q2S1, Q2S2, and Q2S3.

In Fig. 6b, noise source S is turned off and the loudspeakers are energized, one at a time, with a probe signal p. (In most situations, it is possible to shut down or disable the noise source while measurements are being taken.) First, L is activated and data are simultaneously collected at the input (Mι,M₂) and quiet zone (Qι,Q₂) microphones. Next, L₂ is activated and data are collected at the input and quiet zone microphones, and so on, until all four loudspeakers (Lι,L₂,L₃,L4) have been used. Probe signal p can be any well defined wideband signal of brief duration such as an impulse or, preferably, a burst of pseudo-random noise. The duration of probe signal p is usually less than 1 sec, so that it is feasible to repeat the measurement several times.

Figure 6b shows the subset of data resulting from the activation of loudspeakers Lj and L₂. Similar data are obtained from the activation of loudspeakers L₃ and L₄. The data collected in the quiet zone are the L data, which comprise the time series Qι _\p, Q₁L2/?, Q_{1 3}P, Q₁L₄/7, Q_{2 1} ?, Q₂ 2 , Q₂L₃P, and Q2L₄ J. The data collected with the input microphones are the C data, which comprise the time series MJLJ ?, M_\h₂p, MιL₃/?,MjL₄p, M₂L₁/?, M₂L2/?, M₂L₃/?, and M₂L₄P. The C data are used to compute the neutralization filter weights W_c. In SPU 20', the neutralization filters serve to reduce the effects of coupling (acoustic feedback) from the loudspeakers to the input microphones.

The collected data are preprocessed in several ways prior to the computation of W and W_c. As explained above, the original M and Q data obtained under Condition A (Fig. 6a) are windowed at the desired time slots SI, S2, and S3 to provide more compact descriptions of source characteristics. The L and C data obtained under Conditions B (Fig. 6b) are averaged to improve the signal-to-noise ratio of these measurements. Because the forward filter weights calculation code requires the impulse response of the loudspeakers

(rather than their response to probe signal p), the averaged L data are deconvolved with probe signal p. In addition, the Q data and either the L or the M data are filtered to select the optimal forward filter frequency band.

After windowing, averaging, deconvolving, and filtering have been completed in step 52 (Fig 4), the data records available for the forward and neutralization filter weight calculations are as depicted below: Q1L Q1L Q1 L Q1 L

Lij = 1 2 3 4

Q2L Q2L Q2L Q2L 1 2 3 4

Q1S Q1 S Q1S

Qit = 1 2 3

Q2S Q2S Q2S

1 2 3

M1L1p M1L2p M1L3p M1L4p

Ckj = M2L1p M2L2p M2L3p M2L4p

70&P4

where the subscripts have the following meanings: i = quiet zone location; j = loudspeaker number; k = input microphone number; t = time slot number.

The next step 54 (Fig. 4) of computing W and W_c will be discussed with reference to Fig. 7. In general, the performance of a multichannel noise control system can be described by the relation: r = Q + XW, where r is a vector describing the residual noise in the quiet zone, Q is a vector describing the noise source as measured in the quiet zone, X is a crosspoint matrix obtained by convolving microphone data with loudspeaker data, and W is a vector containing all forward filter weights. This relation states that the residue is the sum of the source output and the system output. The residue exists because the system output cannot perfectly match the source output. The objective of the filter weight computation is to make the residue as small as possible. The Wiener filter solution for W can be used to minimize the energy in the quiet zone:

W = - (XT )-i X Q where X^τ is the transpose of matrix X, and X and Q are as defined previously. The Wiener filter solution, as mentioned previously, is well known to those skilled in the art. Introduction to Random Signals by Robert Brown and Patrick Hwang, John Wiley & Sons Inc. Publishers, 1992 (2nd Ed.), includes a discussion on Wiener filters. To illustrate the computation of forward filter weights W, it is convenient to begin with the single input, single output, single channel case. Filter weights for multichannel systems can be calculated in the same fashion, but the computational complexity increases significantly as the numbers of channels is increased.

The following calibration data are available for weight calculations: π = [ lι l₂ l3 ... 1a ],

Qn = [ qi q2 Q3 - q_c ]>

where Li ₁ is the loudspeaker time series of length "a" corresponding to quiet zone location 1 and loudspeaker number 1; Mπ is the input microphone time series of length "b" corresponding to input microphone number 1 and time slot number 1; and Qπ is the quiet zone time series of length "c" corresponding to quiet zone location 1 and time slot number 1.

For the single input, single output, single channel case, W„ = - ( Xι,τχ₁₁)-ιχ₁₁τ_{Ql l =} [ w, W₂ W₃ ... _T ] where τ is the number of taps on the delay line representing the digital filter, and X„ = L„M„, which is set forth below in matrix form:

706-P5

The convolution matrices Lπ and Mn are calculated from time series Ln and Mn prior to computing Xn. A necessary condition for selecting the number of taps τ is that a + b + τ - 2 > 2.

In Fig. 7, the calculation of the filter weights begins with step 70, where the length of the forward filters (τ) is chosen. As mentioned earlier, the length of the forward filters is one of the determinants of the filter characteristics. In step 72, a weight calculation method is chosen. For small systems (less than 12 channels), a direct weight calculation method is preferred because of its relative simplicity and accuracy. For larger system, (12 channels or more), an iterative weight calculation method is preferred for speed. In step 74, a decision is made whether or not to use system constraints. If not, the forward filter weights W are calculated in step 76 (using the method selected in step 72), and the neutralization filter weights W_c are calculated in step 77. If system constraints are used, then the forward filters weights W are computed in step 78 (using the method selected in step 72), and a comparison is made in step 80 to determine whether or not the constraints have been satisfied. If the constraints are satisfied, step 77 is performed immediately. If the constraints are not satisfied, the calculation parameters for W are updated in step 82, and steps 78, 80, and 82 are repeated iteratively until the constraints are satisfied. Finally, the neutralization filter weights are computed in step 77. Step 70 selects the length of the forward filters, which is defined as the number of delay-line taps τ. Filter length is determined in part by the bandwidth requirement of the noise control system and in part by the need to treat delayed reflections in a multipath environment (e.g., a room with reflecting walls). However, as is evident from the above discussion, long filters are computationally very demanding. In the preferred embodiment of this invention, τ is in the range 100 - 500.

In step 72 it is determined whether a direct or an iterative method will be used to compute filter weights. As noted previously, the direct method is fast enough for noise control systems with fewer than 12 channels; the iterative method is needed for larger systems, which are computationally more demanding and would take too much time to calibrate with the direct method.

The option to use constraints in step 74 allows for the calculation of forward filter weights W while constraining certain quantities. Either a peak power constraint or an energy constraint can be applied to weight amplitude, loop gain, or loudspeaker input. A peak power constraint limits the maximum value of the quantity to be constrained. For example, a peak power constraint can be placed on the loudspeaker inputs so that the loudspeakers always operate in their linear range. An energy constraint limits the sum of the squares of the elements of the quantity to be constrained.

Steps 76 and 78 compute the forward filter weights using the method selected in step 72. The direct method will be discussed with reference to Fig.

8, and the iterative method will be discussed with reference to Fig. 9.

As indicated in Fig. 8, the direct method begins with step 84, the calculation of the crosspoint matrix

X = ML.

Next, in step 86, the normal equation

X^TXW = -X^TQ

is calculated. Then, in step 88, an efficient Block Toeplitz package is used to solve for W:

In Fig. 9, an iterative method 90 in accordance with the present invention begins with step 92, where the forward filter weight vector W is set to a preselected initial value. This initial value can be zero, but user or system experience may allow the selection of a non-zero initial value which converges more quickly on a proper filter weight. Next, in step 94, the crosspoint matrix

X = ML

is calculated. In step 96, an error score (i.e., a measure of how well the weights are converging) is computed from

error score = X^TXW + X^TQ,

and in step 98, a decision is made as to whether this error score exceeds a set threshold. If the error score exceeds the threshold, in step 100, the direction vector v is updated at each iteration by its previous value and the current iteration's error score. A new forward filter value W is calculated from

W = Woid + X^τXv.

Steps 96, 98, and 100 are repeated until the error score no longer exceeds the threshold.

As detailed in this invention, the iterative method for weight calculation is primarily an implementation of a conjugate gradient algorithm that achieves a significant speed-up by exploiting the structure of the crosspoint matrix. The problem can be written as

XW = -Q,

where X, W and Q are as previously defined. To solve for the weights in a least-squares sense, the normal equation

X^TXW = - X^TQ is formed, and a conjugate gradient algorithm is used to solve for W. This method for solving linear equations is well known to those skilled in the art. Matrix Computations, by Golub and Van Loan, The Johns Hopkins University Press, 1989 (2nd Ed.) is one of many references which discusses the algorithm.

Most of the computational time needed to solve for W is spent in the calculation of X^τXv. For purposes of this invention, v = v in the calculation of W in step 100 of Fig. 9, and v = W in the calculation of the error score in step 96 of Fig. 9. Both of these quantities are computed for every conjugate gradient iteration. The benefit of this invention is that it allows for the quick calculation of X^τXv by taking advantage of the structure of matrix X. In the frequency domain, matrix X has a block-diagonal structure, which allows factoring of matrices M and L. The normal equation can then be rewritten in a form which requires fewer computational steps in its solution.

The factorization referred to above will be illustrated for the single input, single output, single channel case. To simplify the writing of the matrices, assume a time series length of 2 for Mn, Ln. and Q - Also, assume a tap length (τ) of 2 for vector Wi i . Because X is the convolution of M_{\ \} and Li i , the equation

XW = -Q

can be rewritten as

Mιι*Lπ*W„ = -Qπ.

The convolutions can be written as the following matrix equation

) )

708-P11

or more succinctly as MLW = -Q.

Next, the Toeplitz structure of Mn and Ln is used to compute the product (Mj _\L_{\ \})^τM_{\ \}L_\ i W quickly for any W. To do this, M_{\ \} and Li i are extended to be circulant, and Wj i is padded with zeros:

Mud) 0 0 M (2)

M₁₁(2) Mud) 0 0

M 11 -

0 n(2) Mud) 0

0 0 Mn(2) Mud)

70B-P12

Notice that the extended matrices and vector still match the vector Qi _\ above because the zeros added to Wn will multiply all extended elements of Mj i and ^" _{\ \}. It is well known to those skilled in the art that once matrices have been extended to be circulant, they can be diagonalized using a simple Fast Fourier Transform (FFT) in C(xlog₂X) operations, where C is a constant. The matrices can be factored such that nO) o 0 0

0 Mn(2) 0 0

Mn = F^τ

0 0 Mn(3) 0 0 0 M11(4)

Lnd) o o 0

0 L₁₁(2) 0 0

Ln = FT

0 0 Ln(3) 0

Ln(4)

706-P13 where M] j and L] i are the Fourier transforms of the time series Mi i and Ln. and F is the Fourier matrix generated by

F_jk = n-l/2(ζ(j-O(k-l)),

where j and k are indices that go from 1 to n, n is the matrix order, and

ζ _ _ei2π/n

Because the matrix F is orthogonal, we can rewrite MLW = -Q as

M ^"

708-P14

The normal equation can now be rewritten as

708-P15

and solved using a conjugate gradient algorithm. The advantage of writing the normal equation in the above form is that the per-iteration cost of calculating XTXW is C(xlog2X) rather than C(x²).

The matrix calculations for the single input, single output, single channel case show the advantage of writing the normal equations in a form that offers a per-iteration calculation savings. Although mathematically more complicated, all results can be extended to the multichannel case. In practice, all calculations can be done much more simply and quickly by first calculating the standard FFT of Mn , Ln , Wn , and Qn . and then performing the convolutions in the frequency domain.

In Fig. 10, the simulate performance step 56 of Fig. 4 is illustrated in greater detail. The simulation takes the transfer functions S, M, C, Wc, W, and L and combines them with summers 102, 104, and 106 to create a simulated residual r. More specifically, M and C are combined in summer 102, and the output of summer 102 is combined with the output of Wc in summer 104. The output of summer 104 is applied to filter W, and the output of W is applied to Wc, C, and L. The output of L is combined with S in summer 106 to derive the simulated residual r. When this simulated residual is small, (i.e., the simulated system gives adequate performance for the noise control problem at hand), the filter weights W and Wc are considered able to meet system requirements and suitable for loading into a hardware signal processor (step 58 of Fig. 4).

Claims

1. A method of calculating filter weights for a compression wave cancellation system having multichannel digital signal processing means, said method comprising: specifying a set of system parameters (48); deriving calibration data (50) from a first set of measurements taken with a plurality of input transducers (18) in response to compression waves (12) from a noise source (16), and from a second set of measurements taken with said plurality of input transducers in response to an artificially generated compression waves (p) produced for the purpose of calibrating said compression wave cancellation system; and calculating (54) a plurality of forward filter weights from said system parameters and said calibration data.

2. A method of calculating filter weights as recited in claim 1 wherein said input transducers are coupled to inputs of a plurality of forward filters (F) of a multichannel digital signal processing means (20), and wherein said artificially generated compression waves are generated by a plurality of output transducers coupled to outputs of said plurality of forward filters, and further comprising a step of deactivating said multichannel signal processing means prior to deriving said calibration data.

3. A method of calculating filter weights as recited in claim 2 further comprising a step of preprocessing (52) said calibration data to control the bandwidth of said compression wave cancellation system.

4. A method of calculating filter weights as recited in claim 2 wherein said step of specifying a set of system parameters includes: specifying a size and shape of a quiet zone; specifying a desired attenuation level; specifying a number of input transducers; and specifying a number of output transducers.

5. A method of calculating filter weights as recited in claim 2 wherein said step of deriving calibration data comprises: deriving input transducer calibration data; deriving noise source calibration data; and deriving output transducer calibration data.

6. A method of calculating filter weights as recited in claim 5 wherein said step of deriving calibration data further comprises: deriving coupling calibration data to calculate weights for a plurality of neutralization filters provided by said multichannel digital signal processing means.

7. A method of calculating filter weights as recited in claim 5 wherein said step of deriving an output transducer calibration data is accomplished by exciting said output transducers with a wideband probe signal.

8. A method of calculating filter weights as recited in claim 6 wherein said step of preprocessing said calibration data comprises: windowing input transducer calibration data and noise source calibration data at selected time slots that characterize said noise source; averaging output transducer calibration data and coupling calibration data to improve the signal-to-noise ratio; and filtering a subset of the calibration data to select a frequency band for system operation.

9. A method of calculating filter weights as recited in claim 1 wherein said step of calculating forward filter weights comprises an iterative process.

10. A method of calculating filter weights as recited in claim 1 further comprising: simulating the performance (56) of said compression wave cancellation system on a digital computer, said simulation using said set of system parameters, said calibration data, and said forward filter weights to predict the simulated performance that would be achieved if the simulated system were implemented in hardware.

1 1. A method of calculating filter weights as recited in claim 10 further comprising: repeating said steps of specifying a set of system parameters, deriving calibration data, and calculating forward filter weights if the simulated performance does not meet system requirements.

12. A noise cancellation system comprising: a multichannel digital signal processing means (20) having a plurality of forward filters (F), each of which can be programmed with at least one weight; a plurality of output transducers (22) coupled to outputs of a plurality of said forward filters, said plurality of output transducers being arranged to create a quiet zone (66) in a defined region of space; a plurality of input transducers (18) coupled to inputs of a plurality of said forward filters, said plurality of input transducers being arranged between a noise source (16,S) and said output transducers; means for deriving calibration data (50) for said compression wave cancellation system; and means for calculating weights (54) for said forward filters from said calibration data.

13. A noise cancellation system as recited in claim 12 wherein said means for deriving calibration data includes digital computation means responsive to measurements of compression waves by said plurality of input transducers.

14. A noise cancellation system as recited in claim 13 wherein said means for calculating weights for said forward filters comprises digital computation means which utilizes said calibration data.

15. A noise cancellation system as recited in claim 14 wherein said weights are calculated by an interative method implemented on said digital computation means.

16. An iterative method of calculating filter weights for a compression wave cancellation system comprising the steps of: deactivating a multichannel digital signal processing means (20) having a plurality of input transducers coupled to inputs of a plurality of forward filters (F) and a plurality of output transducers (22) coupled to outputs of a plurality of forward filters; deriving calibration data (50) from a first set of measurements taken with said plurality of input transducers in response to compression waves from a noise source (16), and from a second set of measurements taken with said plurality of input transducers in response to artificially generated compression waves (p) produced by said plurality of output transducers for the purpose of calibrating said compression wave cancellation system; and setting the weights W of said forward filters to a preselected initial value; computing a crosspoint matrix X from said calibration data; computing an error score e from said filter weights and crosspoint matrix X; using said error score e to recalculate said filter weights W if said error score e exceeds a given threshold; repeating the previous two steps until said error score e is at or below a predetermined threshold; and loading resulting filter weights W into said forward filters.

17. An iterative method of calculating filter weights as recited in claim 16 wherein said step of computing a crosspoint matrix X is computed as

X = ML, where M comprises calibration data relating to said plurality of input transducers and L comprises calibration data relating to said plurality of output transducers.

18. An iterative method of calculating filter weights as recited in claim 17 wherein said step of computing an error score e is computed as:

e = X^TXW + X^TQ

where X^τ is the transpose of X, and Q is a vector describing the noise source.

19. An iterative method of calculating filter weights as recited in claim 18 wherein said step of using said error score to recalculate said filter weights produces a new filter weight W from an old filter weight W_cid according to:

where the v is a direction vector calculated using the error score e.