This disclosure is related to: (1) U.S. patent application Ser. No. 09/927,784 filed on even date herewith and entitled “SOUND PROCESSING SYSTEM INCLUDING FORWARD FILTER THAT EXHIBITS ARBITRARY DIRECTIVITY AND GRADIENT RESPONSE IN MULTIPLE WAVE SOUND ENVIRONMENT” in the name of Erik W. Rasmussen, and commonly assigned herewith; (2) U.S. patent application Ser. No. 09/927,783 filed on even date herewith and entitled “SOUND PROCESSING SYSTEM INCLUDING WAVE GENERATOR THAT EXHIBITS ARBITRARY DIRECTIVITY AND GRADIENT RESPONSE” in the name of Erik W. Rasmussen, and commonly assigned herewith; and (3) U.S. patent application Ser. No. 09/928,229 filed on even date herewith and entitled “SOUND PROCESSING SYSTEM THAT EXHIBITS ARBITRARY GRADIENT RESPONSE” in the name of Erik W. Rasmussen, and commonly assigned herewith.
The present invention relates generally to audio signal processing. More specifically, the present invention relates to an audio processing system that exhibits an arbitrary directivity and gradient response.
There are many instances where it is desirable to have a system capable of receiving information from a particular signal source where the environment includes sources of interference signals at locations different from that of the information signal source. For discussion purposes, the specific instances will be generalized to the extent possible. Turning first to
In general, the sound field into which the system 10 is placed contains not only the sounds to be picked up, referred to as a utility signal, but also unwanted sounds, referred to as noise or noise signals. In these situations, it is desirable to process the signals picked up by the microphones 12 in order to reduce the noise contents electronically. There are several conventional methods for reducing the noise electronically through the audio processor 14. One is known as static beamforming where the signals from two or more microphones are passed through filters and combined to form a single signal. The resulting signal will show a sensitivity to sounds that depends upon the direction of the sound incidence as compared to the direction of the microphone assembly. The directivity response will take the form of one or more beams. Due to the fact that the lobes of the directivity response have different magnitudes, the beamformer will show a signal to noise improvement when the beam is oriented so that the utility signal falls within the main lobe and the main part of the noise falls outside the main lobe. Static beamformers have the disadvantage that, in order to provide substantial noise reduction under general noise conditions, a large number of microphones are required.
Another conventional method is known as adaptive beamforming which is achieved when the filters of a beamformer are variable and controlled by an adaptation process. Normally such an adaptation process works to minimize the output signal power. An adaptive beamformer can track noise sources and dynamically adjust the directivity response such that the sensitivity at the direction of the noise incidence is minimized while keeping the sensitivity at the utility direction high. Currently known adaptive beamformers show the disadvantage that they are only capable of tracking a limited number of noise sources, mostly only a single. Furthermore adaptive beamformers work with a fairly large time constant in the adaptation process. Therefore they are only able to track quasi-static noise sources.
Yet other conventional methods apply only to a single microphone multiband noise reduction situation, that is, spectral subtraction. When only a single microphone signal is available, it is not possible to obtain information as to the direction of sound incidence from the microphone signal and it is therefore not possible to perform beamforming as above. Still a reduction of the noise contents can be achieved under these circumstances. Such methods all rely on dividing the signal into a number of frequency bands. In each band the signal is analyzed statistically to derive measures of the utility signal and noise content. Based on these measures a band gain is derived and applied that amplifies bands with utility signal contents while attenuating bands with noise contents. Unfortunately, the statistical analysis requires long time constants. Therefore the single microphone methods are limited to a sound field with stationary noise signals and non-stationary utility signals.
Further conventional methods are known that use two microphones and analyze the microphone signal contents to derive and apply a gain in frequency bands. The gain is a function of how the microphone signals relate to each other. The known methods have the disadvantage that they only work when the signal in each frequency band consists of either utility signal or noise and not when a combination of utility signal and noise is present.
By contrast to the conventional methods above, the present invention uses a different approach to the problem. It uses the general equations for sound fields to analyze the microphone signals and find required properties of one or more components or waves contained in the input signals. The desired properties can for example be the direction of sound incidence or the pressure gradient of the impinging waves. The incoming waves are amplified with a gain function based on these properties, that is, the directivity or the gradient. Based on the amplified waves an output signal is generated either by synthesizing the amplified waves or by applying filtering to an input signal combination. The present invention can operate in a number of applications including hearing aids, directional microphones, microphone arrays, silicon microphone assemblies, headsets, hearing protectors, cordless phones, mobile phones, camcorders, personal computers, laptops, palmtops, and personal digital assistants, among others. In some embodiments, the present invention is especially suited to work with head worn microphones that pick up the speech signal of the wearer. In this application, the present invention offers a substantially improved noise reduction when compared to conventional solutions with comparable sound quality.
A sound processing system including at least one microphone, an audio processor, and at least one output device is disclosed. The audio processor includes an analog beamformer, a microphone equalizer, and an apparent incidence processor. Two different embodiments of the apparent incidence processor are disclosed, that is, a wave generation method and a forward filtering method. Both embodiments use the same principles to estimate the properties of the individual waves of the sound field. With the present invention, it is possible to implement arbitrary directivity or gradient responses using a small number of microphones only, that is, two or three microphones. The present invention offers improved noise reduction also for environments with many independent noise sources. Furthermore, the present invention works for signals and noises with arbitrary statistics.
The accompanying drawings, which are incorporated into and constitute a part of this specification, illustrate one or more embodiments of the present invention and, together with the detailed description, serve to explain the principles and implementations of the invention.
In the drawings:
Embodiments of the present invention are described herein in the context of a sound processing system, including a forward filter, that exhibits an arbitrary directivity and gradient response in a single wave sound environment. Those of ordinary skill in the art will realize that the following detailed description of the present invention is illustrative only and is not intended to be in any way limiting. Other embodiments of the present invention will readily suggest themselves to such skilled persons having the benefit of this disclosure. Reference will now be made in detail to implementations of the present invention as illustrated in the accompanying drawings. The same reference indicators will be used throughout the drawings and the following detailed description to refer to the same or like parts.
In the interest of clarity, not all of the routine features of the implementations described herein are shown and described. It will, of course, be appreciated that in the development of any such actual implementation, numerous implementation-specific decisions must be made in order to achieve the specific goals of the developer, such as compliance with application- and business-related constraints, and that these specific goals will vary from one implementation to another and from one developer to another. Moreover, it will be appreciated that such a development effort might be complex and time-consuming, but would nevertheless be a routine undertaking of engineering for those of ordinary skill in the art having the benefit of this disclosure.
In accordance with the present invention, the components, process steps, and/or data structures may be implemented using various types of operating systems, computing platforms, computer programs, and/or general purpose machines. In addition, those of ordinary skill in the art will recognize that devices of a less general purpose nature, such as hardwired devices, field programmable gate arrays (FPGAs), application specific integrated circuits (ASICs), or the like, may also be used without departing from the scope and spirit of the inventive concepts disclosed herein.
The figures and equations in this document will contain signals and variables that are vectors or matrices. Unless otherwise noted all operations on these vectors and matrices are to be interpreted as being performed element-by-element. For example, a multiplication of two vectors should, unless otherwise noted, be interpreted as an element-by-element multiplication and not as a vector multiplication. For the sake of clarity, the figures have been reduced to a minimum number of elements, however, as one of ordinary skill in the art will recognize, the number of elements will vary with the particular application.
Turning now to
The method of apparent incidence processing involves complex signal operations. It is therefore possible that the processing will be performed with digital techniques. Thus the microphone signals will be converted to digital signals with at least one analog to digital (A/D) converter 24 and the output signal will be converted back to an analog signal, if needed, with a digital to analog (D/A) converter 26.
The analog beamformer 18 provides analog preprocessing according to conventional techniques of the microphone signals that enables the reduction of the resolution of all but one of the A/D converters 24. This can save size and reduce the power consumption. For hearing aids, for example, these properties are highly desirable. Depending on the circumstances however, the analog beamformer 18 may be deleted as unnecessary or too costly.
The method of apparent incidence processing requires, as generally do conventional beamformers, that the microphones 12 of
Turning now to
In a embodiment, the beam amic(1), for example, is formed as the sum of all microphone inputs and the other beams are formed as the difference of a specific microphone signal and a reference microphone signal (Microphone(1)). The filter transfer functions in the laplace domain that implements such a beamforming are as shown in (1) below. In this way, as many analog beam outputs will be processed as there are microphone units, N1=N2. This analog beamforming is well suited for use with microphone arrays implemented as silicon transducers on a single piece of silicon with small transducer spacing.
With this beamforming, the generalized beam amic(i) will relate to the microphone signals as follows:
In another embodiment, all microphone signals are passed directly to the A/D converters without beamforming, amic(i)=Microphone(i), N2=N1.
In yet another embodiment, each of the filters, Filter(j,i), approximates different time delays of the microphone signals, that is, either inverting or non-inverting.
In a further embodiment, only two of the filters, Filters(j,i), for each beam i are present.
In yet a further embodiment, each of the beams amic(i) implements the same directivity using different microphone combinations.
In a still further embodiment, a variation of the last described, the microphones are placed equidistant along a common axis. The analog beamformer is defined according to (3) below. The numbering j of the microphones follows their placement along the common axis with number one being closest to the sound source. NB-1 is commonly referred to as the order of the directivity.
Turning now to
To perform the equalization, one of the microphone signals is chosen as a reference. For convenience, the reference signal is referred to with index one, mic(1). This reference signal is by definition equalized and is thus passed through the processing of the microphone equalizer 20 unaltered. To provide for the necessary equalization, the equalization functions generated by the microphone equalization updater 34 should follow the definition of (4) below.
In (4), S(1) is the sensitivity defined as the digital value at the input of the microphone equalizer 20 of
Turning now to
To derive a magnitude update, the phase of the quotient is analyzed in a zero phase condition detector 46. If the phase indicates that the sound incidence is from a direction and a distance for which the current microphone signal should have the same magnitude sensitivity as the reference microphone signal, then the zero phase condition detector 46 will output a logic one as a ZeroPhase output signal. If the phase condition does not hold, then the analyzer will output a logic zero. When the quotient magnitude minus one is gated with this ZeroPhase switch signal and scaled with a MagnitudeCoef coefficient, a magnitude update value is obtained that is suitable to use to update the magnitude accumulator 44 for the equalizing function.
To derive a phase update, both the phase and the magnitude of the quotient are analyzed. Under normal conditions, it will generally not be possible to derive any information regarding a misfit of the phase of the microphone equalization response. But once in a while, triggered by a specific input from a specific direction, the microphone signals will relate to each other in a way that can only be possible if the equalization response is incorrect. In a compute excess phase monitor 48, the signals are monitored and if such an “impossible condition” is found to exist the amount by which the phase is un-natural will be output as an excess phase signal. If the phase conditions are natural, the compute excess phase monitor 48 will output zero. As the excess phase conditions depend upon signal frequency, the compute excess phase monitor 48 estimates the frequency for its use. The phase update signal is scaled with the PhaseCoef coefficient.
To further improve the quality of the magnitude and phase update signals it is gated with an Inband signal output of a signal inband detector 50. The signal inband detector 50 outputs a logic one if the power in the current frequency band is contributed mainly by input contents of frequencies falling within the band. Conversely, the signal inband detector 50 outputs a logic zero if the contents are due mainly to input at frequencies outside of the band. It is widely known that most time-to-frequency transforms “spill” energy from the source band to neighboring bands due to windowing effects or similar mechanisms. However, only signals within the frequency band should be allowed to influence the equalizing value for a band. This differentiation is possible through the use of the Inband signal.
With the equalization processing described here, the equalization response will be updated dynamically. The phase and the magnitude of the equalization response will be regulated independently. Updating follow through statistical processes that rely on the noise-like nature of the signals that are most likely to be encountered as inputs to the audio processor, such as speech, machine noise, wind noise, and the like. The update signals will contain large noise components, therefore they are scaled with small coefficients, that is, PhaseCoef and MagnitudeCoef, respectively, such that the adaptation times are slow. The use of the coefficients prevents noise from entering the forward signals through the equalization processes.
In an embodiment, the phase and magnitude accumulators 42 and 44 are divided into static and dynamic parts, where the updates only influence the dynamic parts. The effective equalization response is the product of the static and dynamic parts.
In yet another embodiment, the static part of the equalization response is measured with standard measuring techniques once or regularly at the time of production test or at some other convenient times and saved.
In a further embodiment, a forgetting factor is included with the dynamic part of the accumulator. The forgetting factor causes the dynamic response to converge towards zero when no updates are received.
In yet a further embodiment, means are provided that can save the accumulated equalization response when the audio processor is powered down and read the saved response again when the processor is powered up the next time.
In a still further embodiment, the signals mic(i) are all omni directional and the zero phase condition detector 46 is implemented so as to compare the magnitude of the phase with a constant value. If the phase magnitude is smaller than the constant, then a logic one ZeroPhase signal is generated.
In another embodiment, the signals mic(i) are all omni directional and the compute excess phase monitor 48 generates the phase update signal according to (5) below. In (5),f(i) is the frequency as estimated by the compute excess phase monitor 48. In (5), a,f, and ExcessPhase are all vectors covering the frequency range of the frequency transformation used. The * operation in (5) denotes an element-by-element multiplication and not a vector multiply. d(i) is the physical spacing between Microphone(i) and the reference microphone, Microphone(1). c is the speed of sound. ∈ is a small positive constant.
In still another embodiment, the center frequency, f(i), for the equalization is estimated as the center frequency of the band i.
In yet another embodiment, the frequency, f(i), for the equalization is estimated as in (6) below. In (6), k is the frequency band index, fc(k) is the center frequency of band k, BW(k) is the bandwidth of band k, and b is a positive constant.
In a further embodiment, the signal inband detector 50 for each frequency band evaluates the absolute value of its input signal in the current band and the two nearest neighboring bands. If the current band carries the highest absolute value, then the Inband signal for the current band is generated as a logic one and otherwise it is generated as a logic zero.
Turning now to
Turning now to
In an embodiment, analysis beamforming is not performed and thus the signals cmic(i) are directly copied to abeam(i) with N2=N.
Another embodiment is described with the set of filters described in (7) below.
In yet another embodiment, the analysis beamforming of (7) is combined with the analog beamforming of (1) above. In such a combination, it can then be shown that the output of the analysis beamformer 52 will be estimates of the microphone outputs as described in (8) below.
abeam(i)=cADMicrophone(i) (8)
In (8), cAD is the A/D converter conversion gain. With this embodiment the sequence of the processing is preferably changed so that the analysis beamformer 52 precedes the microphone equalizer 20 of
In another embodiment, each of the filters 58, Filter(j,i), approximates different time delays of the microphone signals, that is, either inverting or non-inverting.
In yet another embodiment, the latest embodiment described is changed so that only two of the Filters(j,i) for each beam i are present.
In a further embodiment, each of the beams, abeam(i), implements the same directivity using different microphone combinations.
In yet a further embodiment, a variation of the last described, the microphones are placed equidistant along a common axis. In this case, the analysis beamformer 52 is defined according to (9) below. The numbering j of the microphones follows their placement along the common axis with number one being closest to the sound source. NB-1 is commonly referred to as the order of the directivity.
In still another embodiment, the analysis beamforming is performed in frequency bands.
Turning now to
The equation solver 68 utilizes mathematical functions. Such functions can be included either through table look-up, Taylor-series approximation, or the like. In any case, the dynamic range of the functions may be limited due to hardware constraints. In order to gain maximal use of a limited mathematical dynamic range, the input signals are normalized in the normalizer 66. However normalization may not be necessary or desirable, in which case the normalizer 66 may be omitted. When implemented, the output from the normalizer 66 carries not only the normalized frequency domain signals but also information about the amount by which the signals have been normalized, an exponent, collected in the output signal BeamExp. This exponent enables the recovery of the absolute values from the normalized values. Each beam input and frequency may be equalized independently but the same exponent may also be used across all beams, frequencies, or both.
In an embodiment, the analysis filters 62 perform differentiation with respect to time. Note that differentiation with respect to time can be expressed in the frequency domain with the transfer function of (10) below. In (10), j is the imaginary unit.
Hdiff(ω)=(−j·ω)D (10)
where D is the order of differentiation.
In another embodiment, the analysis filters 62 use the difference equation of (11) below to approximate first order differentiation with respect to time. In (11), n is the sample index and FS is the sampling frequency.
y(n)=(x(n)−x(n−1))·FS (11)
In yet another embodiment, only a single analysis filter 62 per abeam is included.
In still another embodiment, a single analysis filters 62 is included to filter the first abeam signals only, that is, abeam(1).
In an embodiment, the forward transformers 64 are FFT based.
In another embodiment, the forward transformers 64 are performed with a time domain filterbank.
In yet another embodiment, the forward transformers 64 are performed with a time domain filterbank that delivers quadrature outputs from which phase information can be extracted.
In an embodiment, the microphone equalizer 20 of
In an embodiment, the normalizer 66, independently for each frequency band, finds the complex component, real or imaginary, of the ABEAM(i) signals, with the largest magnitude. A common exponent for all beams is found using this largest component. All beams are then normalized with the common exponent.
In another embodiment, the equation solver 68 uses floating point arithmetic. In such a case, the normalizer 66 also converts each of the beams to floating point notation.
Turning now to
In (11) and (12), a is the analysis filter index and i the abeam index.
The equation solver 68 may include a time domain integrator 76 which is optional. When implemented, it integrates product factors P(i1,a1)*P(i2,a2) over time. Through the use of the time domain integrator 76, it may be feasible to enhance the analysis especially for any embodiment using time domain filterbanks as the forward transformers 64 of
The equation solver 68 most importantly includes a core solver 78 which solves the sound field equations. A sound field can be described in several ways. The description can be in the time domain or in the frequency domain, among others. Furthermore, the description can involve a potential field describing both pressure and velocity with the same function or the description can have distinct pressure and velocity equations.
In this case, it is preferred that the sound field be described primarily in the frequency domain and simplifications will be used when judged feasible. Initially, the sound field will be described with two functions. The first function is a complex scalar function. It gives the sound pressure as a function of frequency and position. The second function is a complex vector function. It gives the sound particle velocity as a function of frequency and position.
To move forward towards a solution, it is necessary to make some assumptions regarding the sound field. First, the sound consists of M waves. The waves are not required to be plane waves. Second, only the sound pressure and the direction of the sound particle velocity are of interest. The value of the sound particle velocity is not important. Both the sound pressure and the sound particle velocity of the individual waves will be functions of time. Third, each of the waves in the description may in reality consist of the sum of several waves. The principle of superposition holds for sounds up until very high sound pressure levels. Fourth, one is only interested in the sound pressure levels at the places of the microphones and only as seen through the initial beamforming, A/D conversion, and so forth. Fifth, each of the waves is quasi-sinusoidal in the sense that within each frequency band of the analysis the energy is mainly due to a single sinusoid only. Thus each wave may consist of several sinusoids, spread over the frequency range.
Given these assumptions, the sound field can be described with a set of equations (13) below that gives the sound pressures at the locations of the microphone sound inlets.
In (13), the small dots denote element by element multiplication and the large dots denote inner products. c is the speed of sound. k is the frequency band index. ωm is the angular frequency of the wave m. m is the wave index. x(i) is a vector and is the position of microphone i. vm is a unit vector in the direction of the sound particle velocity. vm is thus the direction of sound incidence of the wave. δm is the damping factor along vm. Note that all wave parameters generally will depend upon k, the frequency band index. This dependency stems in part from the fact that it is assumed that each wave in the description can be the sum of more than one factual waves in the sound field. However, it also accounts for windowing effects and other non-idealities associated with the frequency transformation used. Note also that the wave parameters will be functions of time. In (13) and the equations to come, the dependencies upon the frequency band index and time are implicit except when otherwise noted.
For convenience, the sound field equations (13) above will be rewritten in (14) below as it can be observed through the input signals that are supplied to the core solver 78.
In (14), the effects of the microphone and A/D converter sensitivities, among others, have been included in the values Am. a indexes the Analysis Filters Ha. G(i,vm) is a function that collects all variations due to wave damping, microphone directivity, analog beamforming, and analysis beamforming.
In some applications it will be feasible to place the microphones along a common axis. In this case, the method only allows for imperfect detection of the direction of the sound incidence. Only the angle αm, between the microphone axis and the wave direction in the plane in space that contains both vectors, can be detected. For such embodiments, (14) can be rewritten in terms of this angle as in (15) below.
In (15), d(i) is the physical distance from the reference microphone, microphone 1, to microphone i.
As noted above, it is preferred that the sound field be described primarily in the frequency domain. The frequency domain is generally the most advantageous domain in terms of complexity and computing costs. Nevertheless in some cases the processing of the present invention is most feasible performed together with other audio processing applications. If such audio processing runs in the time domain it may prove efficient to implement the apparent incidence audio processing in the time domain as well. In the following, the sound field equations will therefore be stated in the time domain. Below, (16) states the sound field equations in the time domain under the assumptions as used in the formulation of (14) above. In (16), p is used to describe the time domain version of P and n is the sample index.
(16) can generally not be solved for a single sample measurement only. It is necessary to measure over a number of samples and perform some form of averaging. To enable the solution, new measurement signals are defined in (17) below. In (17), indexes a1 and a2 may have a value of zero. B is the number of samples to average over and win( ) is an optional window to weight the measurements with.
The product terms in the right hand side of (17) will according to (16) be the product of constants and two cosine terms dependent on time through the sample index n. The product of two cosine results in a sum of two different cosines, namely at the sum angle and the difference angle. The sum and difference angles may be DC terms or AC terms. When integrating over a large number of samples the AC terms will diminish from (17). Furthermore, as is known from Fourier transform theory, the integral of the product of two cosines goes to zero for large integration times if the cosines are of a different frequencies. Below, (18) shows what remains of (17) when B is sufficiently large.
A subset of (18) will normally suffice to solve for all parameters except the phase, φ. To solve for φ, a subset of (16) combined with a subset of (18) is needed.
In general all parameters of the waves of the sound field will a priori be unknown and the sound field equations must be solved for all parameters. However, some special cases may exist where, for example, a wave is known to come from a specific direction. In other cases, it may prove useful to make certain assumptions regarding wave parameters. When certain parameters are known or assumed a priori, the task of solving the sound field equations becomes less complex. In still other cases, it may be desirable to imply constraints upon the parameters of the waves. This will generally increase the complexity of the task of solving the equations for the wave parameters, but it can simplify the task of mapping the parameters to wave gains.
In an embodiment, the wave frequencies are assumed such that ωc=2πfc, where fc is the center frequency of the individual bands of the frequency transform.
In another embodiment, the waves are assumed to be plane waves with δm equal to zero.
In yet another embodiment, the first wave is assumed to impinge from a target direction.
A further embodiment solves for two waves. The first, signal, wave is constrained to have a direction from within a certain tolerance around a target direction. The second, noise, wave is constrained to have a direction from outside the tolerance field.
A still further embodiment solves for M waves. The solutions are not directly constrained but they are ordered so that the waves impinge further away from the target direction the higher the wave number is.
There exists several different techniques that allow for the solution of (14), (16), or(18) with respect to parameters of the waves that the sound field consists of. The parameters to solve for may include the following:
Am—the wave amplitude (note that Am is a complex number incorporating phase information)
vm—the direction of sound incidence of the wave
ωm—the angular frequency of the wave (ω=2πf)
δm—the wave damping factor
Next, five different techniques will be described for equation solving. The first technique is referred to as direct solving. Under special conditions it is possible to solve (14) above directly using arithmetical methods. Such direct solving yields the wanted parameters as mathematical functions of the input signals to the core solver 78.
In a specific application, a system with two microphones 12 of
In (20) below, (19) above can be solved to yield the wave parameters in terms of the input signals to the core solver 78.
In (21) below the result of (20) has been rewritten to take advantage of the ratio inputs to the core solver 78. Furthermore, the result has been simplified in order to always return values for the wave frequency and angle of incidence that are real valued.
In (21) above a special technique has been introduced. The wave direction has been solved for while disregarding input amplitude information. In some cases it is possible to estimate some wave parameters using mainly the phase information of the input signals. To disregard the amplitude information has the advantage that the estimation will not be vulnerable to sensitivity changes of the microphones and thus a microphone equalization circuit may not be needed.
The second technique for equation solving is referred to as iteration. Unfortunately it is not generally the case that a solution for the sound field equations in (14) above can be found directly. When a solution cannot be found explicitly as a mathematical expression, it is still possible to use numerical methods. A group of numerical methods can collectively be called iterative methods.
Iteration can roughly be said to include the following steps:
Formulate an initial guess.
Compute an error.
Compare computed error to a predetermined limit.
If the error is less than the predetermined limit, then a solution has been found.
If the error is greater than or equal to the predetermined limit, then map the current solution to a subsequent solution.
Repeat computation and comparison steps until a solution is found.
The error is found by subtracting the actual measurement from the results that are obtained by inserting the current solution in the right hand side of (14). The error can be expressed in a mean square sense but it can also be taken as maximal of any of the measurements. However, the error can also be expressed as the relative difference between two successive solutions.
In an embodiment, the wanted directivity response is symmetric around a given axis. The wave parameters are solved using Newton-Rhapson iteration. Parameter error functions are defined as in (22) below.
In (22), the superscript I is the iteration index. Below, (23) gives the equations to find the next step of the iteration.
The iteration stops when the mean square error nerr is smaller than a given value.
In another embodiment, the error cost function, nerr, is defined as the maximal relative difference between the last two iterative solutions as in (24) below.
In a further embodiment, an extra level of iteration is introduced, that is, wave iteration. Wave iteration may be said to include the following steps:
Solve, using iteration or another method, a subset of the sound field equations for a limited number, for example two, of the waves.
Evaluate a cost function werr, for example as defined in (25) below.
If werr is smaller than a predefined constant, then continue with the next step below otherwise increase the number of waves and return to the solution step above.
Set the amplitudes for the resting M−1 waves to a value of zero and the rest of the parameters for these waves to predefined values.
It should be noted that iteration for a solution can have limitations. For one, the process may not converge. Adding appropriate means to the iteration process can solve this problem. However, iteration might also find a “local extreme” where the stop conditions for the iteration are fulfilled even when the solution is not the correct “global solution.”
The third technique for equation solving is referred to as a full parameter scan. Unlike the iteration solution, this method will always find the global solution when properly set up. The drawback is a greater amount of computation. With the full parameter scan, the possible solution ranges are defined and within these ranges grids are set up. The grids are set to implement the wanted resolution for the respective parameters. Then the error cost function, for example as in (22) above, is evaluated for all possible combinations of wave parameters on the grids. All parameter combinations for which the error cost function is smaller than a given threshold are possible solutions.
In a specific embodiment, the microphone configuration is axis-symmetric and the error cost function is defined as in (22). The parameter ranges are set to [0 . . . 2*BeamExp] for |A|, [0 . . . 2*π] for arg(A), [0 . . . π] for α, and [−3 . . . 3] for δ. Parameter scan is performed and the wave parameters are chosen as the set that gives the lowest value for nerr.
In another embodiment, full parameter scan is combined with iteration. Parameter ranges and grids are set up using a coarse grid. All parameter combinations on the coarse grid are used as starting guesses for an iteration leading to a “local solution” for the wave parameters. The local solution with the lowest error cost function is chosen as the correct solution.
The fourth technique for equation solving is referred to as solution screening/optimizing for minimal power. The methods described above may yield more than one possible solution for the sound field equations for a given a set of measurements. Measurement noise from sources such as A/D conversion, microphones, etc. can be sources of ambiguity but the system may also be underestimated. The sources of underestimation are that the number of microphones used is not large enough to solve for the number of waves assumed or that the sound in reality consists of more waves than solved for.
Even when the system is known to be underestimated it may due to cost reasons be attractive to use solution screening/optimization to work around the problem instead of making the system unambiguous. With solution screening the solution with the lowest cost function is not simply chosen. In stead, a threshold for the error cost function is defined. All solutions for which the error cost function is lower than the threshold are deemed possible solutions. From the set of possible solutions the solution is chosen for which a power estimate is minimal.
With this strategy the chosen solution may not be the correct one, but even if the correct solution is not chosen the strategy results in a system gain for noise components equal to or lower than the noise gain that would have been the result if the correct solution was to be used.
A specific embodiment uses the full parameter scan method. All parameter sets for which the error cost function is equal to or lower than the minimal cost function encountered plus a threshold are deemed possible solutions. The solution with the lowest Ptot, (26), is chosen.
The fifth technique for equation solving is referred to as solving for a subset. In some applications, knowledge of the full set of wave parameters is not necessary in order to control the wave gain. A subset only will suffice. This may simplify the task of solving the sound field equations.
In an embodiment, only the wave damping is of interest. Two microphones are used in this embodiment and it is assumed that the sound field consists of a single wave coming from a direction on the microphone axis. In this case, the sound field equations can be simplified as shown in (27) below.
Above, five techniques for solving the sound field equations have been described with respect to the parameters for the waves in the sound field. Two implementations for these techniques will now be described. A first implementation involves the use of conventional computer software. Several brands of application software exist that are capable of solving the sound field equations with either symbolical or numerical methods. These include mathematical programs such as Mathematica, Matlab, Maple, and the like. In addition circuit simulators may be used under certain conditions.
An embodiment of the system includes a standard computer architecture to do parts of the acoustical sound processing. The sound field equations are defined and solved for the wave parameters within a conventional software package on a conventional computer architecture.
A second implementation involves the use of a conventional table look up. A table can be pre-computed that contains the optimal solutions with a certain resolution. The table can be computed using any of the solving techniques described. Once the table has been computed, one simply looks the solution up in the table to solve the sound field equations. Adding an iteration or approximation process to the look up process can enhance it by minimizing the size of the storage used for the table.
Turning now to
Turning now to
Turning now to
The IsSignal and IsNoise switch signals are multiplied with the wave powers, that is, the squared wave amplitudes. The wave power are summed over all waves in a signal summer 104 and a noise summer 106. The summed signal and noise powers are low pass filtered in a NarrowBand filter 108 to yield narrow band estimates of the signal and noise powers. The effective integration time of the filter 108 controls the speed of the measurement. It must be set large enough that inaccuracies in the wave parameter estimates are filtered out. The narrow band measurement may thus be relatively slow.
To obtain a faster measurement of the signal and noise powers may also be made with a coarser frequency resolution. The WideBandPowers output provides the same measurements as the NarrowBandPowers output with the exception that the measurement has been integrated over wide bands in sum over bands integrators 114 before being low pass filtered in WideBand filter 110. Due to the wide bandwidth the measurement may be performed at a faster rate, that is, a shorter integration time, and with a smaller delay than the narrow band measurement. Note that the dynamic characteristics of filters 108 and 110 control the update speed of the power signals. Therefore the filters will in general have different characteristics.
In an embodiment, the signal or noise analysis 102 is based on the measured direction of sound incidence. If this is within a given tolerance equal to a target direction, then the wave/frequency pair is judged to be signal and otherwise it is judged to be noise.
In another embodiment, the signal or noise analysis 102 is based on the measured direction of sound incidence. The IsNoise signal is generated with the help of a directivity function as shown, for example, in the polar plots of
In yet another embodiment, the signal or noise analysis 102 is based on the measured wave damping. If this is greater than a given threshold, then the wave/frequency pair is judged to be signal and otherwise it is judged to be noise.
In an embodiment, an additional path, generating a second NarrowBandPowers signal, is provided. The two NarrowBandPowers are generated with two different update rates.
Turning now to
The gain function appliers 120 can be thought of as the heart of the processing. They directly control the gain of each wave as a function of some or all of the wave parameters including the direction of the sound incidence, the wave damping, and the frequency and amplitude. It is thus here that the directivity of the processing is implemented. The gain function chooser 118 selects the gain function that serves best for the current signal to noise ratio in view of the strategy input that has been chosen. The output of the gain function chooser 118 will typically control the width of the main lobe of the directional response.
In an embodiment, two gain strategies are implemented, that is, both omni directional and directional. The strategy chooser 116 compares the wideband signal power with the wideband noise power. The omni directional strategy is chosen in all narrow frequency bands covered by wide bands where signal power is greater than a predefined constant times the noise power. In all other bands, the directional strategy is chosen.
In another embodiment, only the directional strategy is implemented.
In an embodiment, the gain function appliers 120 operate by comparing the direction of incidence with a target direction. If the direction of incidence is within a predefined tolerance, the cut-off angle, from the target direction, then the raw gain is set to a predefined maximal gain and otherwise the raw gain is set to a predefined minimal gain. This results in a directivity as shown, for example, in the polar plot of
In an enhancement of the embodiment just described above, the gain function chooser 118 outputs the cut-off angle as a GainSelector signal. The cut-off angle is controlled as a function of the narrowband band signal to noise ratio.
In another embodiment, the cut-off angle is controlled so as to produce a wide mainlobe of the beam for poor signal to noise ratios and a narrow mainlobe for good signal to noise ratios.
Conversely, in yet another embodiment, the cut-off angle is controlled so as to produce a narrow mainlobe of the beam for poor signal to noise ratios and a wide mainlobe for good signal to noise ratios.
In a further embodiment, the gain function appliers 120 operate by table look-up and optional approximation/interpolation in a similar way as the table look up process described with respect to
In a still further embodiment, the GainSelector signal switches between a finite number of different gain functions, for example, those implemented by different tables. An example of a set of gain versus direction functions is shown in the polar plot of
In a far field embodiment, the gain function is chosen to attenuate waves where the absolute value of the wave damping is greater than a predefined threshold. Thus waves are attenuated that are not far field waves.
In a near field embodiment, the gain function is chosen to attenuate waves where the value of the wave damping is lower than a given threshold. Thus far field waves are attenuated.
In another embodiment, only the wave(s) with the largest relative absolute amplitude(s) are amplified and the rest of the waves are attenuated.
The wave estimation and gain control processes will normally be performed on a block of samples. The duration of the blocks will be so large that it is possibly that the raw gain for a specific frequency band in consecutive blocks will differ significantly. Unfortunately, an abrupt gain change may cause unwanted audible effects. Therefore it will generally be desirable to prevent abrupt gain changes. This is the purpose of the gain smoother 96 of
In an embodiment, the gain smoother 96 of
In another embodiment, the gain is smoothed in the gain smoother 96 of
In yet another embodiment, the gain is smoothed with exponential averaging of the raw gains of successive blocks in the gain smoother 96 of
Turning now to
In an embodiment, the wave estimation is performed in the frequency domain. The reverse analysis transformer 122 is thus FFT-based.
In another embodiment, the wave estimation is performed in time domain filter bands. The reverse analysis transformer 122 is then implemented as a reconstruction filter bank.
In yet another embodiment, the wave estimation is performed in the time domain. The reverse analysis transformer 122 is thus omitted.
In an embodiment, the output of the gain mapper 98 is in the frequency domain. The forward gain transformer 126 is thus FFT-based.
In another embodiment, the output of the gain mapper 98 is in time domain filter bands. The forward gain transformer 126 is thus implemented as an analysis filter bank.
In yet another embodiment, the output of the gain mapper 98 is in the time domain. The forward gain transformer 126 is thus omitted.
Turning now to
In an embodiment, the signals are generated in the frequency domain. The wave generator 128 thus merely has to output the complex amplitude, Am or abs(Am)*exp(j*φm). Then the gain is applied by multiplying with the frequency domain gain. Finally, the reverse signal transformer 134 performs an inverse frequency transform.
In another embodiment, the signals are generated in the time domain with sine wave generators for the wave generators 128. The reverse signal transformer 134 is omitted.
In yet another embodiment, the signals are generated in the time domain in narrow frequency bands by the wave generators 128. Then the gain is applied by multiplying in the bands. Finally, the reverse signal transformer 134 is implemented with a reconstruction filter bank.
Recall from above that two different embodiments of the apparent incidence processor 22 of
Turning next to
Turning now to
The forward beamformer 138 is optional and may be deleted leaving the input signals to be directly connected to the signal filter 142. When implemented, the forward beamformer 138 serves to remove noise from the signal thus enhancing the noise reduction performance achieved by the wave parameter controlled gain of the FF gain controller 140. The processing is similar to that of the analysis beamformer 52 of
In an embodiment, the forward beamformer 138 is identical to the analysis beamformer 52 of
In another embodiment, the forward beamformer 138 is in principle identical to the analysis beamformer 52 of
In a further embodiment, two fbeam signals are generated, that is, an omni directional signal and a narrow beam, for example a supercardioid.
Turning now to
Turning now to
In an embodiment, the fbeam and rbeam signals are implemented in the frequency domain. The forward transforms 154 are thus not implemented.
In another embodiment, the correlator 156, the power filter 160, and the integrator and limiter 162 are implemented in the time domain. The rbeam and fbeam signals are likewise also implemented as time domain signals and the forward transformers 154 are omitted. As a result, the gain mapper 164 merely windows the filter response.
In yet another embodiment, the integrator and limiter 162 includes a forgetting factor causing the integrated response to tend towards zero during periods of no signal activity.
In a further embodiment, two microphones are positioned along the target axis. The forward beams include an adaptive beam fbeam(i1). pbeam(i1) is implemented with a beamformer 146 of
In (29) and (30), e is constant in the range 0.5 to 1 and d(2) is the microphone spacing.
In a still further embodiment using at least one adaptive forward beam, in
Turning now to
The signal filtering in the forward filtering embodiments can be based on already beamformed input signals. The FF gain controller 140 therefore has to compensate for the directivity and near field characteristics of the forward beams. The BeamPar signal carries enough information about the forward beam that a plurality of FF gain function appliers 166 can compute a gain that implements the target directivity as shown, for example, in the polar plots of
Turning now to
In an embodiment, all of the forward beams are static. The beam dependency has been included with the pre-computed tables for the gain function chooser 118 of
In another embodiment, the FORWARDGAIN(i) signal is taken as the maximal of the wave gains GAINRAW(j,i).
In yet another embodiment, the FORWARDGAIN(i) signal is taken as the minimal of the wave gains GAINRAW(j,i).
In still another embodiment, the FORWARDGAIN(i) signal is the power weighted average of the individual wave gains as defined in (31) below.
In a further embodiment, at least one static forward beam signal and at least one adaptive beam signal are implemented. The FF gain function applier 166 monitors and analyses the BeamPar signal from the adaptive beam over time. When the BeamPar signal is stable, indicating a significant noise signal from a constant direction, then the GainSelector signal is switched to use the adaptive beam mainly to build the output. When the BeamPar signal resembles random noise, then the GainSelector signal is switched to remove the adaptive beam from the output.
Returning to
In addition to the single output embodiments of the apparent incidence processors 22 of
In an embodiment of either method, two outputs are generated. The first output contains the sum of the near field waves while the second output contains the sum of the far field waves.
In another embodiment of either method, there are N5 different outputs generated.
Each output contains the sum of the waves originating from a specific range of directions.
In yet another embodiment of either method, again there are N5 different outputs generated. Each output contains the sum of the waves originating from a specific range of directions. The wide band power of the waves in the sound field is measured. The individual output generation blocks are controlled in a way such that the first output is always generated using a range of directions centered around the origin of the wave with the largest power.
To this point of the discussion, embodiments have been described that employed either the wave generation method or the forward filtering method. It is also possible to combine the two methods. One simply replaces the forward filter 136 of
Turning now to
In a specific embodiment utilizing both the wave generation method and the forward filtering method, the WGFF gain controller 176 chooses the gain function so that: at high signal to noise ratios, forward filtering is the primary contributor to the output; and
at low signal to noise ratios, wave generation is the primary contributor to the output.
In the embodiments so far the process of finding gains for the waves of the sound has included two main steps, that is, finding the parameters of the waves and deriving a gain based on the parameters found. Both these main processes can be described by mathematical transforms, as depicted in (32) below, and in many cases they are best implemented using techniques known from mathematical pocket calculators Mathematically the gain may be described as a transform directly of the inputs as described in (33) below.
Turning now to
Assume that a specific example exists that uses two microphones and that uses forward filtering of the first microphone signal without any beamforming. A single table lookup implements the combined core solving of the sound field equations and the gain control. No statistical evaluation is performed. The inputs to the table lookup are the magnitude and the phase of the ratio Q(2) of the quotient signal obtained when, in the frequency domain, the second microphone signal is divided by the first microphone signal. The phase of Q(2) is quantized to one of thirty two possible phases covering the total complex phase range. The magnitude of Q(2) is quantized to one of 512 possible magnitudes covering the range from 0.01 to 100. The gain is stored as a binary value of either one or zero. The table thus implemented requires 16 Kb of storage capacity.
For applications where speech is to be picked up from the wearer of a headset, as in a mobile phone, a hearing protector, or the like, two main directions regarding noise suppression have until now been used. The most effective method has been the use of so-called noise-canceling microphones. These microphones amplify near field signals while attenuating far field signals. Unfortunately, noise-canceling microphones have to be placed no farther than two to three centimeters away from the speech source in order to be effective. This may not always be possible or convenient. Another method has been to use directional microphones pointing towards the mouth of the wearer. Unfortunately, a directional microphone can make no distinction between near and far field signals and thus it will not offer as large of a noise reduction as is possible with a properly placed noise-canceling microphone. A preferred near field embodiment of the present invention enables signal processing methods with which it is possible to produce sound pick-up with near field characteristics. It is possible to obtain noise reduction better than that possible with noise-canceling microphones. Furthermore, it is possible to maintain the near field characteristic with its noise reducing virtues at a distance further away from the speech source than is possible with conventional noise-canceling microphones.
The near field method works by dividing the input signal into a number of frequency bands. In each band, the input signals are analyzed to see whether the activity in that band is due to near field sources or to far field sources. If the activity is from near field sources, then that band is replicated in the output with a high gain and otherwise it is replicated with a low gain.
Turning now to
Turning now to
Turning now to
In an embodiment, the accumulator 210 is divided into static and dynamic parts, where the updates only influence the dynamic parts. The effective equalization response is the product of the static and dynamic parts.
In another embodiment, the static part of the equalization response is measured with standard measuring techniques once at the time of production test or at some other convenient time and saved.
In yet another embodiment, a forgetting factor is included with the dynamic part of the accumulator 210. The forgetting factor causes the dynamic response to converge towards zero when no updates are received.
In a further embodiment, means are provided that can save the accumulated equalization response when the near field audio processor is powered down and read the saved response again when the processor is powered up the next time.
In yet a further embodiment, the microphones used have the same directivity and frequency response except for the small tolerances that the microphone equalizer 200 of
In still another embodiment, the signal inband detector 212 for each frequency band evaluates the absolute value of its input signal in the current band and the two nearest neighboring bands. If the current band carries the highest absolute value, then the Inband signal for the current band is generated as a logic one and otherwise it is generated as a logic zero.
Turning now to
In an embodiment, M microphones are placed along a common axis. The microphone signals are pair wise beamformed to yield N=M1 beams with identical directivity response.
In a variation of the embodiment above, the beams are supercardioids.
In another variation, the beams are figures of eight.
In another embodiment, the beamforming is performed in the time domain with FIR or IIR filters.
In yet another embodiment, the beamforming is performed in the frequency domain.
Turning now to
In an embodiment, the gain control processing is performed on blocks of samples. For each block, a single complex signal value per frequency band is computed. The power and phase filters 220, 222 only use the values from the current block to compute their respective outputs.
In another embodiment, the power and phase filters 220, 222 averages the signal powers and phases over consecutive blocks.
In yet another embodiment, the phase averaging is power weighted.
In still another embodiment, no phase information is utilized.
In a further embodiment, the forward transformer 218 is implemented with a time domain filterbank, no phase information is generated or used, and the signal powers are measured with a finite time constant.
In yet a further embodiment, the forward transformer 218 is FFT based.
Turning now to
The IsSignal and IsNoise switch signals are multiplied with the wave powers, that is, the squared wave amplitudes. The weighted signal and noise powers are low pass filtered in NarrowBand filter 230 to yield narrow band estimates of the signal and noise powers. The effective integration time of the filter 230 controls the speed of the measurement. It must be set large enough that inaccuracies in the wave parameter estimates are filtered out. The narrow band measurement may thus be relatively slow.
To obtain a faster measurement of the signal and noise powers may also be made with a coarser frequency resolution. The WideBandPowers output provides the same measurements as the NarrowBandPowers output with the exception that the measurement has been integrated over wide bands in sum over bands integrators 240 before being low pass filtered in WideBand filter 232. Due to the wide bandwidth the measurement may be performed at a faster rate, that is, a shorter integration time, and with a smaller delay than the narrow band measurement. Note that the dynamic characteristics of filters 230 and 232 control the update speed of the power signals. Therefore the filters will in general have different characteristics.
In an embodiment, the signal or noise analyzer 228 is based on the R(2) signal. If this signal is less than a predefined threshold, then the signal is judged to be utility signal and otherwise it is judged to be noise.
In another embodiment, there is an additional path that generates a second NarrowBandPowers signal. The two NarrowBandPowers are generated at two different update rates.
The near field gain function applier 226 of
In an embodiment, two microphones 12 of
Turning now to
In an embodiment, two microphones are used, no statistical evaluation is performed and neither is the signal phase evaluated, and a single magnitude threshold is implemented. In this case, the near field gain function can be written as in (36) below.
In a variation of the previous embodiment, the phase is evaluated as well. If the phase of the two microphone signals differs too much, then the band will probably contain energy from more than one source and thus be noisy, in which case, a small gain is assigned. Below, (37) shows the gain function for this situation.
In another variation of the embodiment described above, the narrow band powers are evaluated. In this case, the gain function can be described with (38) below.
Turning now to
In a specific embodiment, the WidebandPowers are monitored. At good signal to noise ratios, all of the signals are passed through without attenuation. At poor signal to noise ratios, a near field characteristic is used.
In another embodiment, the NarrowBandPowers are monitored. Depending on the value of the signal to noise ratio, gain functions of different widths are chosen. For example, the gain function of different widths are shown in the graph of
In an embodiment, the wideband signal power is compared with the wideband noise power and two gain strategies are implemented, that is, both omni directional and directional. The omni directional strategy is chosen in all narrow frequency bands covered by wide bands where signal power is greater than a predefined constant times the noise power. In all other bands, the directional strategy is chosen.
Returning to
While embodiments and applications of this invention have been shown and described, it would be apparent to those skilled in the art having the benefit of this disclosure that many more modifications than mentioned above are possible without departing from the inventive concepts herein. The invention, therefore, is not to be restricted except in the spirit of the appended claims.
Number | Name | Date | Kind |
---|---|---|---|
3057960 | Kaiser | Oct 1962 | A |
4589137 | Miller | May 1986 | A |
4628529 | Borth et al. | Dec 1986 | A |
4750207 | Gebert et al. | Jun 1988 | A |
4759071 | Heide | Jul 1988 | A |
4802227 | Elko et al. | Jan 1989 | A |
4837832 | Fanshel | Jun 1989 | A |
5097510 | Graupe | Mar 1992 | A |
5170434 | Anderson | Dec 1992 | A |
5263019 | Chu | Nov 1993 | A |
5285502 | Walton et al. | Feb 1994 | A |
5289544 | Franklin | Feb 1994 | A |
5305307 | Chu | Apr 1994 | A |
5347496 | Clark et al. | Sep 1994 | A |
5412735 | Engebretson et al. | May 1995 | A |
5473684 | Bartlett et al. | Dec 1995 | A |
5473701 | Cezanne et al. | Dec 1995 | A |
5511128 | Lindemann | Apr 1996 | A |
5553134 | Allen et al. | Sep 1996 | A |
5586191 | Elko et al. | Dec 1996 | A |
5602962 | Kellermann | Feb 1997 | A |
5648936 | Clark et al. | Jul 1997 | A |
5651071 | Lindemann et al. | Jul 1997 | A |
5677987 | Seki et al. | Oct 1997 | A |
5680467 | Hansen | Oct 1997 | A |
5689572 | Ohki et al. | Nov 1997 | A |
5694474 | Ngo et al. | Dec 1997 | A |
5706394 | Wynn | Jan 1998 | A |
5825898 | Marash | Oct 1998 | A |
5836790 | Barnett | Nov 1998 | A |
5978490 | Choi et al. | Nov 1999 | A |
6023514 | Strandberg | Feb 2000 | A |
6157403 | Nagata | Dec 2000 | A |
6160757 | Tager | Dec 2000 | A |
6272229 | Backgaard | Aug 2001 | B1 |
6317703 | Linsker | Nov 2001 | B1 |
6351529 | Holeva | Feb 2002 | B1 |
6549630 | Bobisuthi | Apr 2003 | B1 |
6738481 | Krasny et al. | May 2004 | B2 |
7046812 | Kochanski et al. | May 2006 | B1 |
20010038699 | Hou | Nov 2001 | A1 |
Number | Date | Country |
---|---|---|
195 40 795 | May 1997 | DE |
0 679 044 | Oct 1995 | EP |
0 795 851 | Sep 1997 | EP |
0 831 458 | Mar 1998 | EP |
0 883 325 | Dec 1998 | EP |
0 942 628 | Sep 1999 | EP |
1 091 615 | Apr 2001 | EP |
1 065 909 | Jan 2002 | EP |
2 202 942 | Oct 1988 | GB |
9-146443 | Jun 1997 | JP |
11-18191 | Jan 1999 | JP |
2000-47699 | Feb 2000 | JP |
WO 9508248 | Mar 1995 | WO |
9534983 | Dec 1995 | WO |
WO 9729614 | Aug 1997 | WO |
WO 9750186 | Dec 1997 | WO |
WO 9847227 | Oct 1998 | WO |
WO 9904598 | Jan 1999 | WO |
WO 9909786 | Feb 1999 | WO |
WO 9939497 | Aug 1999 | WO |
9945741 | Sep 1999 | WO |
WO 9952211 | Oct 1999 | WO |
WO 9953336 | Oct 1999 | WO |
WO 0019770 | Apr 2000 | WO |
WO 0033634 | Jun 2000 | WO |
WO 0041436 | Jul 2000 | WO |
0110169 | Feb 2001 | WO |
WO 0169968 | Sep 2001 | WO |
WO 0172085 | Sep 2001 | WO |
03015459 | Feb 2003 | WO |
03015460 | Feb 2003 | WO |