The Stationary Phase Approximation, Time-Frequency Decomposition and Auditory Processing

The principle of stationary phase (PSP) is re-examined in the context of linear time-frequency (TF) decomposition using Gaussian, gammatone and gammachirp filters at uniform, logarithmic and cochlear spacings in frequency. This necessitates consideration of the use the PSP on non-asymptotic integrals and leads to the introduction of a test for phase rate dominance. Regions of the TF plane that pass the test and don't contain stationary phase points contribute little or nothing to the final output. Analysis values that lie in these regions can thus be set to zero, i.e. sparsity. In regions of the TF plane that fail the test or are in the vicinity of stationary phase points, synthesis is performed in the usual way. A new interpretation of the location parameters associated with the synthesis filters leads to: (i) a new method for locating stationary phase points in the TF plane; (ii) a test for phase rate dominance in that plane. Together this is a TF stationary phase approximation (TFSFA) for both analysis and synthesis. The stationary phase regions of several elementary signals are identified theoretically and examples of reconstruction given. An analysis of the TF phase rate characteristics for the case of two simultaneous tones predicts and quantifies a form of simultaneous masking similar to that which characterizes the auditory system.


I. INTRODUCTION
The principle (or method) of stationary phase (PSP) [1] is a result from asymptotics that can provide closed-form approximations, in the limit as λ → ∞, to often intractable oscillatory integrals of the form where a, b, t, λ ∈ R.There are two steps involved: (i) recognition that in the limit the integral will be almost zero everywhere in the interval {t : −∞ ≤ t ≤ ∞} except near values of t where the derivative ḃ(t) is zero, the stationary phase points; (ii) the integral in the vicinity of these stationary phase points can be expressed in terms of the second derivative of the phase i.e. b(t).Perhaps the most successful application of the PSP in signal processing has been in the context of synthetic aperture radar (SAR), where it is the starting point in the development of many of the Fourier-based imaging algorithms, c.f. [2].Application of the PSP in not without its pitfalls.It is tempting to use the PSP in the non-asymptotic cases where λ = 1 to find closed form approximations to integrals such as Fourier transforms.The argument for this requires that the phase b(t) is changing much more rapidly than the amplitude a(t).However, as pointed out in [3], a degree of care must be exercised, particularly with step (ii).The primary interest here is its application to linear time-frequency (TF) decomposition [4].The motivation is the recent resurgence of interest in analogue filter banks both as part of a synthetic cochlea and as a means to provide power efficient implementations of analysis filter banks [5].The desire with both is to extract salient features from the TF decomposition using the limited functionality associated with analogue circuitry.This does not deny the considerable work that have been done on computational modelling of the auditory system, typified by papers such as [6] and the references therein.However the main purpose of such work is to model and predict the response of the auditory system to stimulus rather than to expose the signal processing principles that might be at work.
The PSP is a natural place to start because of the prevalence of oscillatory terms in TF decompositions and the hope that the stationary phase points may provide a means for identifying salient features as well as a focus for sparse decompositions without the need for the usual iterative re-synthesize steps [7].The PSP has been applied to linear TF decomposition for both analysis, [8] & [9], and synthesis [10], the latter leading to the method of reassignment.The approach here is to revisit [10] and to fundamentally re-interpret it to provide a PSP-based approximation to the TF synthesis integral.There is no attempt to either reassign [10] or relocate [11] components in the TF plane because of the limited functionality mentioned above.Another alternative would be to follow an amplitude-based approach such as ridgelets [12].However this leads to algorithms that are far from feasible with analogue circuitry and further, ridges in the TF plane may not be appropriate when dealing with auditory filters such as the gammatone [13] and gammachirp [14] which have asymmetrical impulse responses and, in the case of the latter, an asymmetrical frequency response.
Concentration on the synthesis rather than the analysis integral is advocated because: (i) most methods for sparse atomic compensation, c.f. [7], are based on re-synthesis of the original waveform; (ii) one of the main functions of the auditory system is to code the incident waveforms and coding requires at least some consideration of the potential for reconstruction (even when reconstruction is not a requirement); (iii) smooth variations of the magnitude and phase of the integrand are more readily satisfied for the synthesis integral than the analysis integral (because the signal has already been filtered in the analysis process).
In Section II a preliminary description of the filter banks to be considered is given.A framework is presented that accommodates Gaussian, gammatone and gammachirp filters at uniform, logarithmic and cochlear spacings in frequency.More details of the relevant properties are provided in Appendix A. Section III is an examination of the application of the PSP to non-asymptotic integrals and proposes the use of step (i), above, without step (ii).The concept of phase-rate dominance is introduced to partition the interval into sets where the PSP is and is not applied.The synthesis double-integral is addressed in section IV and the PSP is applied to it in a new way to facilitate the detection of stationary phase points in the TF plane.Extending the single-integral concepts to double integrals (Appendix B) leads to a test for phase-rate dominance in the TF plane.The stationary phase regions of several elementary signals are identified in section V and examples of reconstruction provided.Section VI examines the phase rate characteristics in the TF plane for the case of two simultaneous tones .It is shown that if these characteristics are used to test for the presence of a tone at a particular frequency, a form of simultaneous masking, similar to that which characterizes the auditory system [15], is observed.All of the above is addressed from the perspective of deterministic signals.Consideration of random processes is left for future work.The effects of noise on the related method of reassignment are covered in [16] and [17].Finally, in section VII, conclusions are drawn.

II. PRELIMINARIES
Consider a time-frequency analysis X ω (t) of a signal of interest x(t) of the form: where * denotes convolution and the impulse response of a single filter in the analysis filter bank is given by: Each filter is formed using a prototype filter h(t) and a nominal bandwidth β.The frequency response of the analysis filter is related to the frequency response of the prototype in a straight forward way, i.e.
Prototype filters that will be considered here are a gammatone pulse, a gammachirp pulse and a Gaussian pulse.While the Gaussian pulse is a common [4] and analytically convenient choice for time-frequency analysis, the gammatone [13] and gammachirp pulses [14] more closely model the cochlea in the ear.Details of these prototype filters and their properties are summarized in Appendix A. The prototype filters are normalised such that H(0) = 1.This property and multiplication by the nominal bandwidth β in (3) normalizes the maximum gain of each filter to unity at ω rad/s.A re-synthesis, x(t), of the signal of interest is performed using filters matched to h ω (t).Thus for real signals of interest and where C is a constant and {.} and {.} denote 'real part of' and 'imaginary part of' respectively.As a convenience, in order to deal with a number of possible filter bank spacings, define a filter bank variable µ that lies in the range [0, 1], A value µ = 0 indicates the lower edge of the filter bank and µ = 1 indicates the upper edge.The frequency ω at which the filters gain is a maximum is a function of µ and the bandwidth of the filter β is proportional to the derivative dω dµ , i.e. β ∝ dω dµ .Hence the total number of filters required to just cover the band of interest is given approximately by dω dµ 1 β .In the following ω min is nominally the lowest frequency covered by the filter bank and ω max is the maximum.For a uniformly spaced filter bank ω = {ω max − ω min }µ + ω min and hence u = {ω − ω min }/{ω max − ω min }, dω dµ = ω max − ω min , the nominal bandwidth β is a constant and dβ dω = 0.For logarithmically spaced filter banks, similar to wavelets [4], ω = ω min e bµ where b = ln( ωmax ωmin ) and hence µ = 1 b ln ω ωmin , dω dµ = bω, the nominal bandwidth β is proportional to ω and dβ dω is a constant.For a cochlear spaced filter banks, based on the approximation of [18] for low sound pressure levels, ω = {ω min + a}e bµ − a, where a = 2π×10 3 4.37 , b = ln ωmax+a ωmin+a and hence µ = 1 b ln ω+a ωmin+a , dω dµ = b{ω + a} and dβ dω is a constant.For all the above filter banks, the equivalent rectangular bandwidth (ERB) [15] and the 3 dB bandwidth are related to the nominal bandwidth β in a straightforward way.

III. THE PSP AND NON-ASYMPTOTIC INTEGRALS
Consider the integral (1) evaluated over an interval [t−∆/2, t+∆/2].Assuming that the function f (t) = a(t)e b(t) is well approximated by its Taylor series over this interval, gives The derivatives can expressed as and: In the asymptotic case, at a suitably large value of λ, ḟ (t) f (t) ≈ λ ḃ(t) and f (t) f (t) ≈ −λ 2 ḃ2 (t) + λ b(t) and hence the integral is dominated by the phase derivatives ḃ(t) and b(t).Thus, at a stationary point, where ḃ(t) = 0, the integral can be evaluated in terms of b(t) without reference to derivatives of a(t).For the non-asymptotic case where λ = 1, these the approximations are also dependent on the relationships between the derivatives of the magnitude and phase of the integrand.Thus ḟ (t) f (t) ≈  ḃ(t) provided, as in [8], that and If the amplitude and phase derivatives are available then at every value of t it is possible to test for what might be called first order or second order phase-rate dominance using (5) or, ( 5) and ( 6), respectively.The approach adopted here, where a closed form approximation to the integral is not required, is to circumvent the difficulties associated with step (ii) of the PSP by simply not using that step.Rather ( 5) is used to identify a set S 0 that is the union of intervals of t where step (i) is valid and the compliment to that set S0 where it is not.If the stationary phase points {t i } i where ḃ(t i ) = 0 and intervals such as S i ∈ {t : t i − δ 1 < t < t i + δ 2 , δ 1 ≥ 0, δ 2 ≥ 0} that contain them can be identified, the integral over the whole real line (1) can be replaced by a similar integral over the union of intervals S ∈ { i S i } ∪ S0 .It is also worth noting that for (5) to be satisfied at or near a stationary phase point, would also require that | ȧ(t)| a(t) | → 0 as t → t 0 .Thus the normalized amplitude rate must go to zero more rapidly than the phase rate, i.e. (6) must apply.Thus there are liable to be be intervals where S i ∈ S0 and for these intervals there is no need to identify δ 1 and δ 2 .

IV. A TIME-FREQUENCY STATIONARY PHASE APPROXIMATION
The PSP can be extended to double integrals such as (4) with the result defined in terms of the gradient and Hessian of the phase of the integrand, c.f. [1], pg.478.However, as in Section III the full form the PSP is not used here.Specifically, the gradient is used in step (i) to identify the stationary phase points but the Hessian of step (ii) is not used to form an approximation to the integral in the vicinity of these points.A test similar to (5) is developed for double integrals in Appendix B to identify regions of phase-rate dominance in the TF plane where the PSP is applied by numerical integration in the vicinity of the stationary phase points.As in Section III, it is convenient to define normalized time-and frequency-derivatives of the integrand Z(τ, µ), i.e.Z τ respectively, where and where Both Z τ and Z µ are additions of a signal dependent term, e.g.
The former is a function of the pair (ω, τ ) whereas the latter is a function of the pair (ω, τ − t).Note that (ω, τ − t) are themselves parameters of the filter, specifically, the frequency where the filter has maximum gain ω and the delay τ − t between the input and output of the filter.In contrast to the method of re-assignment derived in [10], the delay term {τ − t} is interpreted as the group delay of the filter (3) at ω, where the notation ∠H indicates the argument of the complex variable H. Simple expressions for the group delay are given in the Appendix A in terms of the group delays of the prototype filters at zero frequency.The justification for the use of ( 10) in ( 7), ( 8) and ( 9) proceeds as follows: {τ − t} is the delay between the input signal x(t) and the output of the analysis filter X ω (τ ); since each filter is tuned to have a maximum gain a particular frequency ω, the delay through the filter is also tuned to the rate of change of the phase response at the frequency where the gain is maximum.A particular filter is thus jointly labeled with both its frequency ω and the group delay at that frequency g(ω).Thus ( 7) and ( 8) can be written as and respectively, where Then because {η(ω)} = 0 for all three filter types including the complex gammachirp (c.f.Appendix A), the time derivative of the phase of the integrand is and the frequency derivative is Time and frequency derivatives of the analysis integral (2) are constructed using the derivative filters ∂ ∂τ h ω (τ ) and Alternatively, as suggested in [10], they can be derived from the output of (2) by direct differentiation of the signal X ω (τ ) to obtain (15) and by using neighbouring analysis filters to obtain an approximation to (16).Stationary phase points {(ω i , τ i )} i are solutions to: The derivative filters can be expressed in terms of the prototype filter h(t) and its derivative ḣ(t).Expressions for the derivative ḣ(t) of the Gaussian and gammachirp filters can be found in Appendix A. From ( 17) there are two conditions that must be satisfied simultaneously for a stationary phase point to occur at (ω i , τ i ), specifically: 1) the frequency, ∂∠Xω(τ , observed at the output of filter at time τ , is equal to the centre frequency ω of the filter; 2) the delay, ∂∠Xω(τ , observed at the output of the filter at frequency ω, is equal to the group delay g(ω) of the filter at that frequency.Together these define a signal matching condition: at (ω, τ ) the frequency and delay observed at the output of the filter must match the designed centre frequency and group delay of that filter.
Locating stationary phase point requires a grid search over ω for a bank of analogue filters or over both ω and τ , for a discrete-time filter bank.Such a grid search is not onerous since it is implicit in the implementation of the analysis integral.With a grid search there is always the risk of missing the pair (ω i , µ i ) that satisfy (17).This risk can be reduced by: (i) defining a phase gradient vector where the superscript T indicates matrix transpose; (ii) using the Euclidean norm of this vector to construct a test for stationary phase points, i.e.
φ(τ, µ) < C 1 (19) where the threshold C 1 is a small positive real constant.The Euclidean norm is used here for analytic convenience.
Other vector norms may be appropriate and may have desirable properties with respect to ease of implementation.
In addition to finding stationary phase points, equations ( 11) and ( 12) can also be used to test for phase-rate dominance in the TF plane.For phase-rate dominance the inequality must be satisfied, where the amplitude gradient vector is and the threshold C 2 is a positive real constant greater than or equal to one.The projection term is formed from the sum of the projection of the phase rate vector in the direction of the normalised amplitude rate, i.e. φ T a/ a , plus the projection in the orthogonal direction, i.e. φ T 0 1 −1 0 a/ a and hence W = 1 1 −1 1 .The test is derived in Appendix B. Thus the PSP divides the time frequency plane into two regions: a region S 0 where ( 20) is satisfied and the rest of the TF plane S0 where it is not.
Given the stationary phase points {(µ i , τ i )} i that are solutions to (17), the stationary phase approximation is invoked by replacing (4) by: where S is a subset of the TF plane defined as S = { i S i } ∪ S0 , S i is the neighbourhood of the ith stationary phase point such that (µ i , τ i ) ∈ S i and i S i contains all points in the TF plane that satisfy (19).Equation (23) promises sparsity directly from analysis without the computationally expensive re-synthesis step associated with most methods for sparse atomic decomposition.The atomic decomposition of (23) is sparse in the sense that ∀ (µ, τ ) S the coefficient X ω (τ ) is implicitly set to zero.However there are no guarantees about the degree of sparsity that can be achieved or the quality of reconstruction that might be expected apart from the usual ones that might be expected from a well-designed snug or tight frame [20].This will be explored in the following section.

V. ELEMENTARY SIGNALS
This Section is devoted to a consideration of the stationary phase regions of the TF plane associated with some important elementary signals and also in identifying regions of phase-rate dominance.Because of the elementary nature of the signals there is some hope that closed form solutions are possible.These elementary waveforms are: an impulse; a single tone; a linear chirp; a decaying phasor.Closed form solutions that define the stationary phase points are derived.Within the constraints of the available space results of the numerical evaluation of these regions are also provided to confirm and expand upon the theoretical results.
The results are for order n = 4 gammatone (c = 0) and gammachirp (c = 4) filter banks at the cochlear spacing described in Section II.Numerical evaluation considers a frequency band from ω min = 2π × 100 rad/s to ω max = 2π × 5000 rad/s at a sampling rate of 20 kHz.For the gammatone, dω .These values are commensurate with ERB figures for the cochlea, c.f. [18].The spacing of the filters in frequency is such that there are 4 filters with their maximum gain within the ERB of each filter [20].In this case 103 filters are used to cover the stated band.The combined frequency response of the analysis and synthesis filters banks has a linear phase and is flat to within a fraction of a dB over the band.Thus, over the band of interest, the combination of analysis and synthesis filter banks, ( 2) & (4), is effectively distortionless.All filters are approximated by simply truncating them at a suitable point to give finite impulse response (FIR) filters.Anticausal filters (e.g.synthesis filters) are simulated by incorporating suitable delays.
For the inequalities (( 19) and ( 20)), C 1 = 10 and C 2 = 1.The former was found experimentally to provide a good indication of the stationary phase regions for the filter bank described above and for all combinations of filter spacing and prototype filter described in Section II.The latter is an extreme value.The inequalities, (20) of Section III and (54) of Appendix B, suggest a larger value such as C 2 = 10.However a value of unity was chosen because: (i) it makes clear the boundary in the TF plane between regions where the projection p(τ, µ) is greater than the norm a(τ, µ) ; (ii) it reflects a desire to test the degree of sparsity that could be obtained.Thus (20) becomes Studies of ridgelets and skeletons as in [12] suggest that for asymptotic signals, very sparse representations are possible. A

. An impulse
The impulse is obviously important because it is an extreme example of a transient signal.Further it is not an asymptotic signal in the sense considered in [8] and as such illustrates the advantages of applying to PSP to the synthesis integral rather than the analysis integral.For an impulse x(t) = δ(t) the analysis of (2) yields X ω (τ ) = βh(βτ )e ωτ .From which, (11) and (12) give and For a Gaussian pulse g(ω) = 0, η = 0 and ḣ(βτ ) h(βτ ) = −βτ .For a gammachirp pulse g(ω) = n β , η = −1/n and Thus, for all the three pulse type at filter spacing considered, the stationary phase region is defined as τ = g(ω), ∀ω, i.e. a contour in the TF plane at the group delay.
As to phase-rate dominance, the most straightforward case to consider is a uniform filterbank constructed from Gaussian filters for which Z τ = −β 2 τ , Z µ =  dω dµ τ , Φ T = 0 dω dµ τ and a T = −β 2 τ 0 .Using (20) and ( 22), the test for phase dominance becomes Recalling that for uniform filter banks dω dµ = ω max − ω min and thus In the following sub-section this result will be considered again and contrasted with a related result involving the response to a single phasor.More generally it is worth noting that Z µ in (26) is a function of normalized time βτ alone and is not dependent on the filter frequency ω since both dω dµ 1 β and dβ dω are constants for a given filterbank.Z τ , on the other hand, is dependent on ω indirectly through β if the filterbank is non-unform.If Z τ varies with ω and Z µ does not then the relationship between p(τ, µ) and ||a(τ, µ)|| may also be dependent on ω.
The amplitude rate norm is thus which has a zero at βτ = { n−1 n+1 }n, c.f. Figure 1(a).Straightforward manipulation also gives: which, not surprisingly, has a zero at the stationary phase point.Consider two limiting cases of (24).First, as βτ → 0, it reduces to βc βτ > β n−1 βτ , which is satisfied if |c| > n − 1.Thus, for example, in Figure 1(a), for βτ < 3, the inequality is satisfied for c ≥ 4 and phase rate dominance can be achieved to the left of the peak in the impulse response.Second, as βτ → ∞ and c ≥ 0, the inequality reduces to: The LHS is always positive and linearly increasing with βτ , hence a value βτ 0 can always be found above which the inequality is satisfied.For {c : 0 ≤ c < n + 1}, increasing c reduces the RHS and hence reduces βτ 0 .For c ≥ n+1 the inequality is satisfied for all βτ .A degree of caution must be exercised in the application of this result as it is based on asymptotic arguments.However it does predict the trend.As c increases from 0, βτ 0 decreases, increasing the region {βτ : βτ > βτ 0 } where the inequality is satisfied.This trend is evident in Figure 1(a) for βτ > 4. In summary, increasing the value of c increases the region to left of the peak response where phase rate dominance can be achieved.To the right of the stationary phase point this region can also be increased by increasing and positive values of c.The effect of this increase will however, as for the Gaussian pulse, be influenced by the relative sizes of dω dµ 1 β and β, c.f. (33) Figure 2 shows both the complete TF response and results for TFSPA for n = 4 and c = 4.The stationary phased points that satisfy (19) are indicated in white on Figure 2 (b).This stationary phase contour is not co-incident with the peak response or ridge but rather lies at the group delay of each filter as shown earlier.Figure 2(b) is obtained by only plotting |X ω (τ )| at points in the TF plane where (24) is not satisfied.As expected from Figure 1, the leading edge of the response has been removed.Reconstructions of the input waveform are shown in Figure 3 using both full TF plane and TFSPA.In this case they are almost identical despite the removal of the leading edge of the response.

B. A single phasor
Apply a single tone of λ rad/s to the analysis filter bank i.e. x(t) = u(t)e jλt .In the steady state for t 0 this gives where Ω = λ−ω β , with Z τ = β {η + Ω} and Thus for all filters considered at the any of the filter spacings considered the stationary phase point occurs when λ = ω, ∀τ since, by definition, the group delay of the filter at ω is g H(0) .As to phase-rate dominance, the most straightforward case to consider is again a uniform filterbank constructed from Gaussian filters.For a Gaussian pulse, g(ω) = 0, Ḣ(Ω) = −ΩH(Ω) and η = 0. Thus Z τ = βΩ and Noting that Z τ is purely imaginary and For uniformly spaced filters, this will be satisfied if β > dω dµ 1 β .Note that this in direct contradiction to (27) and hence it is not possible to design a Gaussian filter bank that exhibits phase-rate dominance in response to both impulse and phasor inputs.For non-uniform filter banks the RHS is quadratic in Ω. However provided dβ dω |Ω| << 1, the same general conclusion can be drawn for a large range of Ω and hence λ.
For a gammachirp pulse, g(ω) = n/β, Ḣ(Ω)/H(Ω) = −n{Ω+j}−c{Ω+c/n} 1+{Ω+c/n} 2 and η = −1/n.Thus At the stationary phase point where, Ω = 0, As in the previous sub-section, for uniform filter banks or for non-uniform filter banks where |c|  24) is satisfied and phase rate dominance is achieved outside the nominal bandwidth β of the filters, i.e. for |Ω| > 1  2 .The implication is that the output of an analysis filter at ω rad/s, in response to a phasor at λ rad/s, contributes little to the output of the synthesis filter bank if |λ − ω| > β 2 .Further, (24) can be used to locate the vicinity of the peak in the response |X ω (τ )|.
Figure 5 shows both the complete TF response and the results for TFSPA.The stationary phase points are indicated in white on Figure 5(b).The remaining values of |X ω (τ )| that are plotted are at points in the TF plane where ( 24) is not satisfied.The steady state behaviour is well predicted from the theoretical considerations above as illustrated in Figure 4.The horizontal white line at 1000 Hz corresponds with the stationary phase point on Figure 4.There is only a small region around Ω = 0 in Figure 4 where (24)is not satisfied.This region is just visible on Fig. 5 as a horizontal band in red at 1000 Hz.The transient behaviour has similarities to that of the impulse response of Figure 2(b).Specifically the stationary phase points, on Figure 2(b) at the the group delay, are present at most frequencies on the LHS of Figure 5(b) apart from a band of frequencies around the applied one.Figure 6 shows the reconstruction achieved both with the whole TF plane and with TFSPA.The transient response of TFSPA introduces an amplitude modulation that dies away before matching the steady state response.This modulation is evident for τ < 7 ms, the group delay of the filter, after which the stationary phase points at 1000 Hz is established.

C. Linear chirp
Consider a linear chirp with a rate γ rad/s 2 , i.e. x(t) = u(t)e  γ 2 t 2 , where u(t) is the step function at the origin.The response is: for τ 1 β .This integral is intractable and is often approximated using the PSP.However it might appear imprudent to proceed with repeated application of that principle.Rather the approach here is to approximate the integral by more explicit means.First assume that the frequency γt of the quadratic phase term e  γ 2 t 2 is approximately constant over the temporal extent of h(βt), specifically γ β 2 .Then replace the quadratic phase term with its frequency at the filter group delay, i.e. e  γ 2 t 2 ≈ e γg(ω)t , to give The envelope |X ω (τ )| will have a peak when −ω + γτ − γg(ω) = 0, thus τ = ω γ + g(ω).The stationary phase points are defined as the solution to {Z τ } = 0, hence and Substitution for { Ḣ(.)/H(.)}gives the stationary phase points as Note that ġ(ω) is generally negative because g(ω) ∝ 1 β and β decreases with ω (in non-uniform filter banks).More precisely, for a Gaussian pulse g(ω) = 0 and hence so is ġ(ω).For log and cochlear spaced filter banks using gammatone and gammachirp pulses g(ω) = n β and hence ġ(ω) = − n β 2 dβ dω , where dβ dω is positive and dβ dω 1.Thus, for γ β 2 , the denominator term 1 + γ ġ(ω) lies in the interval [0, 1) and thus the stationary phase points lag behind the peak response.When the chirp rate is low, i.e. provided γ β 2 , the behaviour of the p(τ, µ) and a(τ, µ) can be inferred from Figures 4.There will be a band in the TF around the stationary phase line where the (24) is not satisfied.
Figure 7 illustrates these points.The applied signal contains a down-chirp (γ = −10 × 10 4 rad/s/s) which can be observed in Figure 7(a) as a continuous ridge from top left to bottom right.A second lower amplitude up-chirp (γ = 8 × 10 4 rad/s/s) produces a ridge from bottom left to upper right.The down chirp produces a series of stationary phase points (in white) from top left to bottom right.The position of these points is well predicted by (40), shown as a black line, apart from bottom right where γ > β 2 .As expected there is a region of the TF plane around the stationary phase lines where (24) in not satisfied.The lower amplitude up-chirp is depicted in a similar manner, including stationary phase points, apart from the region where the two chirps cross.In this region the larger amplitude chirps hides the trajectory of the lower amplitude chirp.This is a form of simultaneous masking that will be explored more thoroughly in Section VI.Finally there is a low amplitude artifact that lies between the two chirps.This artifact does not contain stationary phase points and could be ignored on that basis in applications where signal analysis was the primary objective.The outputs from the synthesis filter banks are shown in Figure 8 at the point where the two frequencies cross.TFSPA produces some distortion in the reconstruction when compared with that produced using the whole TF plane.

D. A decaying phasor
Voiced speech can be viewed as a signal of the form λi n A i u(t − nT )e λi{t−nT } where T is the pitch period and {λ i } i are the complex formant frequencies with centre frequencies { (λ i )} i and time constants {− (λ i )} i ; A ii are the complex amplitudes of this atomic decomposition.Because of the relationship to speech it is informative to consider the stationary phase points associated with a single atom of this decomposition, i.e. when x(t) = u(t)e λt with a single complex frequency λ, {} < 0. Given this: To proceed further consider the case where h(t) ∈ R and thus restrict consideration to Gaussian and gammatone pulses.The normalised time derivative is  Note that since h(t) is real, the first term on the right hand side is also real if (λ) = ω at which point: {Z τ } = 0, c.f. (11).In a similar manner e ωt e −λt dt which is purely imaginary for Gaussian and gammatone pulses.Thus setting c.f. (12).While it is unlikely that closed form solutions of (41) for τ can be found, nevertheless it is still possible to infer something about the nature of the solutions.Because h(βt)e − (λ)t > 0, ∀t and {t − g(ω)} < 0, ∀t < g(ω), no solution exists for τ < g(ω).For τ ≥ g(ω) the integral can be split into two terms, i.e.
dt.The first term is a negative constant.The second is is nonnegative monotonically increasing function of τ because h(βt)e − (λ)t {t − g(ω)} ≥ 0, ∀τ ≥ g(ω).Thus, if a solution exists it will be the only solution.Such a solution will obviously be dependent on the value of the bandwidth β of the analysis filter for which ω = (λ) and the time constant − (λ) of the applied pulse.Thus it may be possible to infer both the frequency and time constant of the applied pulse from the position of the stationary phase point.

VI. SIMULTANEOUS MASKING
Simultaneous audio masking is the well studied process whereby the presence of one tone prevents the detection of a second tone [15].In this section it is shown that the use of the test (19) leads to a form of simultaneous masking.In particular, if the test (19) is used to indicate the presence or otherwise of a tonal component at ω, the amount by which that test is de-sensitized by the presence of a second tone at ω 1 is dependent on: (i) the frequency separation of the tones {ω 1 − ω} ; (ii) the relative amplitude |A 1 | of the tones; (iii) the magnitude frequency response H {ω1−ω} β of the filter at ω. Consider a signal x(t) = e ωt + A 1 e ω1t made up of two tones at ω and ω 1 , where A, the relative amplitude, is positive real.The steady-state response of the analysis filter at ω is where H 1 = H(Ω 1 ) and Ω 1 = {ω 1 − ω}/β is the normalized frequency separation with derivative dΩ1 dω = − 1 β Ω 1 dβ dω + 1 .The time derivative follows: as does the frequency derivative: It is convenient to rewrite the time derivative (13) of the phase of the integrand as and the frequency derivative (14) as with denominator and numerator of (42) where θ = {ω − ω 1 }τ − φ and φ = ∠{A 1 H 1 }.The angle variable θ indicates the position in time with respect to one period of the beat frequency {ω − ω 1 }.In (43) which, unlike the numerator in (42), is dependent on the particular pulse used.Consider a gammatone pulse for which g(ω and H * 1 Ḣ0 = −nH * 1 from which the normalized frequency derivative is obtained using (46).
After some manipulation this yields equation (52) (see over), where The numerator of (52) is the sum of two oscillatory terms in θ.Assuming that the filterbank is such that β 2 , then the first oscillatory term, from {Z τ }, is larger than the second term, from {Z µ }.
To obtain an approximation to the norm, first consider |A 1 H 1 | 1, in which case the denominator of ( 52) is a constant.For approximately out of phase and thus the minima of the numerator can again be found at the zeros of the larger term, i.e. cos(θ) = −A 1 H 1 .Substitution in (52) yields For |A 1 H 1 | 1 the first term does not go to zero.Given that this time derivative term much larger than the frequency derivative term, an approximation to the norm is obtained by assuming that the latter term is negligible and hence This function of angle has minima at cos(θ) = 1, which are: Together ( 50) and (51) approximate the behaviour of the minima of the norm Φ(τ, µ) as a function of the frequency separation Ω 1 and the amplitude of the the second tone as observed at the output to the filter |A 1 H 1 |.When two tones are present the norm will not go to zero.However, despite that, the test (19) may be satisfied for a given threshold C 1 twice per period of the beat frequency at cos(θ) = −|A 1 H 1 | and thus indicate that a region of stationary phase is present.It is also evident from (50) and (51) that the behaviour of the norm changes significantly at the point where |A 1 H 1 | = 1. Figure 9 illustrates Φ(τ, µ) min as a function of |A 1 H 1 | for a particular filter and a frequency separation, Ω 1 = 1.The minimum value of the norm is evaluated both by numerical minimization of (52) and using the approximations of (50) and (51).There is a clear discontinuity at |A 1 H 1 | = 1 as might be expected from the development of the approximations.The approximation are generally a good fit to the minima evaluated numerically.For this example the most significant errors can be observed around |A 1 H 1 | = 1.It is also clear from the numerical evaluation that the minimum of the norm is a montotonically increasing function of |A 1 H 1 | and as such is invertible.Thus if a particular value of the threshold C 1 in ( 19) is chosen to detect stationary phase points a corresponding value of A 1 H 1 can be calculated that will just achieve a minimum of C 1 .Together (50) and (51) are not guaranteed to provide a montonically increasing function.However minor adjustments in the vicinity of the discontinuity can overcome this and together they provide a simple means of inverting the function.Figure 10 illustrates how adoption of the stationary phase test (19) leads to a simultaneous masking effect.A 0 dB masking tone at 1 kHz is introduced that corresponds to the interferer at ω 1 rad/s.For each filter frequency ω rad/s, (50) and ( 51) are used to calculate the amplitude of a tonal component at ω rad/s that would produce a value of minimum norm equal to the threshold value of C 1 .This amplitude is converted to deciBels to give the detection threshold.For reference, the response of each of the filters, i.e. |X ω (τ )|, to the masking tone alone is shown.When ω 1 = ω a single tone is present and the results of sub-section V-B apply.The norm will go to zero at ω 1 and thus the test will be satisfied.There is a region around ω 1 where the detection threshold increases as |ω − ω 1 | increases.However outwith that region the detection threshold follows the same trends as |X ω (τ )| and, as such, this masking effect is not symmetrical and affects filters above ω 1 more than below it.

VII. CONCLUSIONS
The starting point for this paper was an examination of the application of the PSP to non-asymptotic integrals in general and TF synthesis in particular.The conclusion was that only one aspect of the PSP, location of stationary phase points, is required.The second aspect, approximation of the integral through use of the second derivative of the phase of the integrand, is only needed when closed form expressions that approximate the integral are required.When this requirement is removed, the second aspect can be replaced with a test for phase rate dominance.Regions of the TF plane that pass the test and do not contain stationary phase points contribute little or nothing to the final output.Analysis values that lie in these regions can thus be set to zero.In regions of the TF plane that fail the test or in the vicinity of stationary phase points, synthesis is performed in the usual way.In re-examining the application of the PSP to the TF synthesis integral, a new interpretation of the location parameters associated with the synthesis filters leads to: (i) a test for locating stationary phase points in the TF plane; (ii) a test for phase rate dominance in that plane.With this formulation the stationary phase regions of several elementary signals have been identified theoretically and it has been shown that sparse reconstructions of tones and chirps are possible.An analysis of the TF phase rate characteristics for the case of two simultaneous tones predicts and quantifies a form of simultaneous masking similar to that which characterizes the auditory system.All of the above was addressed from the perspective of deterministic signals.Consideration of random processes is left for future work.

APPENDIX A
This appendix details the relevant characteristics of Gaussian, gammatone and gammachirp prototype filters.All have a maximum unit gain at zero frequency.These prototype filters are scaled in time and modulated in frequency to form uniformly-spaced, logarithmically-spaced and cochlear-spaced analysis and synthesis filter banks.A Gaussian pulse has a non-causal impulse response h(t) = 1 √ 2π e −t 2 /2 with a peak at t = 0, a group delay of 0 and a derivative dh(t) dt = −th(t).It's frequency response is H(Ω) = e − Ω 2 2 with derivative dH(Ω) dΩ = −ΩH(Ω).An order n gammachirp pulse, with chirp rate parameter c has a causal impulse response defined in terms of the complex gamma function Γ(.) as h(t) {1 +  c n } n+c Γ(n + c) t n−1 e −{1+ c n }t+c ln(t) , ∀t ≥ 0, h(t) = 0, otherwise.This is a minor modification to the gammachirp filter of [14] in order to decouple the dependency between the location of the peak gain in the frequency response and the chirp rate parameter c.The time derivative is: The gammatone pulse is a specific case of this when the chirp rate parameter is zero, i.e. c = 0.The gammachirp pulse has frequency response: The maximum gain of the frequency response is 0 dB which occurs at Ω = 0.The derivative of its frequency response is The peak in the magnitude of the impulse response occurs at t = n − 1 while the group delay at Ω = 0 is − ∂∠H(0) ∂Ω = n.The gammachirp pulse, like the Gaussian, is scaled (in time) and modulated to give the analysis filter of 3. The chirp rate parameter is unaffected by the time scaling because of the action of the logarithm function i.e. c ln(βt) = c ln(β) + c ln(t).The time scaling adds a phase shift c ln(β) to the impulse response.Thus values of c can be used interchangeably with [14].From (52), the phase derivative of the impulse response is: d∠h(t) dt = − c n + c t , which is zero at the group delay of t = n.Likewise dh(t) dt = − 1 n at the group delay.
The integral of the function around a point v 0 over an interval S = {v : v < r} with radius r is I S (v 0 ) = S f (v 0 + v)dv.Assuming that the function is well approximated by its first order Taylor series over this interval and making a change of variable v = Qu, where Q is a orthonormal rotation matrix chosen such that Q T ∇ a (v 0 ) = [ ∇ a (v 0 ) 0 ] T and Q T ∇ b (v 0 ) [ ∇ b1 (v 0 ) ∇ b2 (v 0 ) ] T , gives: − sin(θ) cos(θ) , where cos(θ) sin(θ) = ∇a(v0) ∇a(v0) .Thus the integral over the interval will be dominated by the phase variations if |∇ b1 (v 0 ) + ∇ b2 (v 0 )| ∇ a (v 0 ) , where and W = 1 1 −1 1 .Therefore a test for first order phase rate dominance is

Figure 1 (= 1 9 , dω dµ 1 β 1 .
a) illustrates the relationship between p(τ, µ) and ||a(τ, µ)|| for a single filter in a filterbank where dβ dω = 35 and β = 712 rad/s for 3 values of the filter parameter c.The stationary phase point is at βτ = n = 4 whereas the peak magnitude of the impulse response occurs is at βτ = n − 1 = 3, c.f. Figure 1(b).Insight into the relative behaviour of p(τ, µ) and ||a(τ, µ)||, and how that is affected by the value of c, can be obtained by considering either a uniform filter bank, i.e. dβ dω = 0, or, equivalently, a non-uniform filter bank where dβ dω 1 and dβ dω | c n | The latter leads to the approximation

Figure 4 , dω dµ 1 β
Figure 4(a) illustrates the relationship between p(τ, µ) and ||a(τ, µ)|| for a single filter in a filterbank where dβ dω = 1 9 , dω dµ 1 β = 35 and β = 712 rad/s for 3 values of the filter parameter c as in the previous sub-section.As expected from the asymptotic analysis, the norm ||a(τ, µ)|| is approximately constant and virtually independent of c.The projection p(τ, µ) is approximately linear in Ω, as suggested by (38), and much less dependent on c than the response to an impulse.For these examples, (24) is satisfied and phase rate dominance is achieved outside the nominal bandwidth β of the filters, i.e. for |Ω| > 1 2 .The implication is that the output of an analysis filter at ω rad/s, in response to a phasor at λ rad/s, contributes little to the output of the synthesis filter bank if |λ − ω| > β 2 .Further, (24) can be used to locate the vicinity of the peak in the response |X ω (τ )|.Figure5shows both the complete TF response and the results for TFSPA.The stationary phase points are indicated in white on Figure5(b).The remaining values of |X ω (τ )| that are plotted are at points in the TF plane where (24) is not satisfied.The steady state behaviour is well predicted from the theoretical considerations above as illustrated in Figure4.The horizontal white line at 1000 Hz corresponds with the stationary phase point on Figure4.There is only a small region around Ω = 0 in Figure4where (24)is not satisfied.This region is just visible on Fig.5as a horizontal band in red at 1000 Hz.The transient behaviour has similarities to that of the impulse response of Figure2(b).Specifically the stationary phase points, on Figure2(b) at the the group delay, are present at most frequencies on the LHS of Figure5(b) apart from a band of frequencies around the applied one.Figure6shows the reconstruction achieved both with the whole TF plane and with TFSPA.The transient response of TFSPA introduces an amplitude modulation that dies away before matching the steady state response.This modulation is evident for τ < 7 ms, the group delay of the filter, after which the stationary phase points at 1000 Hz is established.