Filtering in the joint time/chirp-rate domain for separation of quadratic and cubic phase chirp signals

This article investigates the possibility and convenience of a filtering operation in the joint time/chirp-rate (TCR) domain, and proposes a novel linear TCR filter for decomposing multicomponent signals into their quadratic and/ or cubic phase chirp components with monotonic instantaneous chirp-rate (ICR) laws only. The TCR domain mask of the filter is selected on a display of a TCR representation of an input signal to isolate the desired chirp component. Projecting the input signal onto the phase signal associated with the TCR mask and approximating the phase difference in this projection operation in terms of ICR values result in the proposed TCR filter that recovers the selected component. Simulations illustrate the proposed filtering in recovery of undersampled cubic phase signals and in resolving back-to-back objects from in-line holograms for which cases it is easier to design filter masks in the TCR domain than in the time-frequency domain.


Introduction
Multicomponent nonstationary signals are widely encountered in many applications including radar, sonar, communications and optics. Parametric methods mostly based on polynomial phase modeling may be used to analyze and estimate such signals and separate them into their components; such as nonlinear least squares techniques [1,2], a maximum likelihood algorithm [3], an expectation-maximization based method [4], an array processing approach based on state estimation via an extended Kalman filter [5], a cyclic moment based method for polynomial phase signals with independent random amplitudes [6], techniques using transforms like high-order ambiguity function [7][8][9] and time-frequency (TF) Hough transform [10][11][12], and an approach for chirplet approximation [13], among other such methods.
The above methods require the number of components in the analyzed signal and/or orders of their polynomial phases to be known or estimated beforehand, although some of them are able to estimate these parameters along the way [1]. Besides, for some applications it is sufficient to decompose the analyzed signal into its chirp components, as in object reconstruction from inline Fresnel holograms [14], without much need for signal model parameters. Nonparametric signal separation methods may be more suitable for such applications, such as a periodicity-based algebraic separation algorithm [15] and an automatic signal separation method based on difference equation representation of chirp signals [16]. The first method requires the number of signal components and relies on inequality of component periodicities [15]. The second one has been reported to give better performance in instantaneous frequency (IF)/ amplitude estimation when applied to monocomponent signals especially for low SNR cases, and, has been suggested to be used after signal component separation by TF filtering in such cases [16]. Hence, TF filtering is still indispensible for many signal separation applications, as reviewed in [17].
As maintained in [35,36], if TF support region of a signal component is nearly disjoint from those of other components and the background noise in an input signal, then, TF transfer function of the Wiener filter that optimally estimates that component reduces to the indicator function of its TF support region. Hence, TF transfer function or pass region of such a filter is selected on a display of a TF representation of the input signal to isolate the desired component [17].
Kozek and Hlawatsch [37] compares linear TF filters to nonlinear TF analysis-masking-synthesis methods based on the Wigner distribution (WD) and the smoothed WD, with prescribed TF pass regions, in TF signal separation problems, and finds that TF filters, in general, yield improved performance with reduced computational cost. Indeed, our simulations indicate that especially Weyl and TF projection filters separate chirps with excellent accuracy. They usually give several percents of error in the noiseless case, where percentage error is defined as the energy of the deviation of the filter output from the desired chirp component normalized by the energy of that chirp.
Despite the good performance and convenience of linear TF filters in chirp separation applications, it may still be more convenient to prepare the mask function of a time-varying separating filter in the joint time/ chirp-rate (TCR) domain, rather than in the TF domain, for some of those applications. One such application is reconstruction of back-to-back objects from in-line Fresnel holograms [14]. Each such object is represented by a pair of lines with opposite slopes in an associated space/ spatial-frequency (or TF) representation obtained from the hologram, intersecting at the object coordinate. Magnitude of the slopes is inversely proportional with the object depth, i.e., the distance of the object to the hologram plane [14]. Thus, linear tracks associated with back-to-back objects overlap in the space-frequency (SF) domain, making it tedious to design SF filter mask functions to resolve such objects. In the space/chirp-rate (SCR) plane, those objects are represented by distinct horizontal strips corresponding to different slopes or depths. Hence, if the mask function of a separating filter can be prepared in the SCR plane to isolate such strips, that would further ease the filter design task.
Motivated by the above application, we propose a novel linear time-varying filter in the TCR domain, reminiscent to TF filtering, for decomposing multicomponent signals to reconstruct their chirp components of the form a(t) exp(j2π(t)), where (t) is a quadratic or cubic phase with a monotonic instantaneous chirp-rate (ICR) law (given by its second derivative). Aside from the mentioned problem, it is also more beneficial to use the proposed TCR filter with its mask function prepared in the TCR plane to recover undersampled quadratic or cubic phase chirps if their ICR curves change more slowly than their IF curves. We present simulations illustrating separation and reconstruction of severely undersampled cubic phase signals with IF curves traversing the discrete TF plane many times within the signal duration while their ICR laws vary much more slowly and exhibit single linear tracks in the TCR plane. For such signals, it is almost impossible to design a TF mask function for a separating TF filter, but a TCR mask can be easily prepared for the proposed TCR filter.
The idea of filtering in the joint TCR domain is novel. Filtering schemes based on domains other than frequency and TF domains have been developed before; however, they are not directly related to the TCR domain. As a previous work of this kind, [38,39] have proposed an extended FT (EFT) matching a known IF function and have developed a time-varying filter for reconstruction of signals with that known IF, by masking in the EFT domain and then taking the inverse EFT. Similarly, a filtering operation in the frequency modulation (FM) rate parameter domain has been performed for suppressing interference chirp signals with nonlinear phase functions in direct sequence spread spectrum communication systems [40]. A linear transform with a kernel that matches these phase functions maps these signals to impulses in the FM rate domain. Then, undesired chirps can be masked out and desired chirp components or the spread spectrum sequence can be recovered after an inverse matched signal transform [40]. Both methods employ an analysis-masking-synthesis approach.
Although mask design in the parameter domain is time-invariant in these methods, they yield time-varying filters with suitable TF transfer functions for separation of selected signal components. The former method requires a positive IF function [38,39], whereas the later does not have this restriction [40]. Both of them require the IF function of the signal to be reconstructed to be known up to a scaling constant. Another filtering approach has been developed in [41] to filter dispersive guided wave signals, based on unitary operators matching this kind of signals; however, proposed TF filters are specifically tailored for and limited to these Pekeris guided waves [41].
Unlike the two methods above [38][39][40], our proposed TCR filter does not employ an analysis-masking-synthesis procedure. Instead, a TCR domain mask function H(t, a) is prepared to enclose the TCR signature of the desired signal component, where a denotes the chirp-rate (or frequency-rate) parameter. Then, this mask function is transformed to yield a time-varying impulse response, as TF transfer functions are transformed to obtain time-varying impulse responses of linear Weyl and Zadeh TF filters [17]. To enable such a TCR filtering operation, firstly, a TCR representation of an input signal should be computed and displayed so that a desired chirp component of it can be identified on the TCR display and can be selected by a TCR mask that encloses its linear ICR trace.
Our proposed filter is more flexible than the above ones in that it does not require the knowledge of the IF or ICR functions of the desired component up to a constant but works for any linear strip-shaped TCR pass region. However, it can separate only quadratic or cubic phase signals with monotonic ICR laws exhibiting single linear tracks, as will be verified in the article. Piecewise linear ICR components can be recovered with repeated use of the filter for each linear segment.
Several TCR representations that can facilitate a TCR filter can be found in the literature. A generalized WD that serves as a joint time-phase derivatives representation for monocomponent, constant-amplitude polynomial phase signals has been proposed in [42], based on decomposition of polynomial derivatives in terms of shifted versions of the involved polynomial. O'Neill and Flandrin [43] has presented a quartic, shift-invariant TCR representation. O'Shea [44,45] have proposed the cubic phase function (CPF) for estimating phase parameters of cubic phase signals. A product CPF has been proposed for multicomponent chirps, in [46], to eliminate spurious peaks appearing in the CPF when applied to such signals. Extended versions of the CPF have been developed [45,47,48] to estimate polynomial phase signals with higher order phases. Finally, a class of joint time-phase derivatives distributions highly concentrated along phase derivative curves has been derived in [49].
Among the above, [43,49], beyond estimation of phase parameters, have also used their transforms as joint TCR representations or distributions in the form of two-dimensional (2-D) images that display ICR curves of analyzed signals.
In our article, we employ the CPF [44,45], the quartic TCR distribution of [43], a bilinear TCR distribution in [49], and a shifted version of a quadratic local polynomial periodogram [50][51][52] to obtain our TCR displays on which desired signal components are identified and masked.
Section 2 derives the proposed TCR filter by approximating the phase difference in terms of the second derivative of the phase, i.e., ICR values, while projecting an input signal onto the phase signal associated with the TCR pass region of the filter. One of the terms in the derived time-varying impulse response requires an approximate knowledge of the IF value of the desired signal component at a reference time instant, in the form of an IF distribution at that time instant. It should be selected on a TF display of the input signal. Hence, the proposed TCR filter is based on joint use of a TCR representation with a TF representation displayed as images.
Section 3 derives the equivalent Weyl TF transfer function for the filter with an infinitesimally narrow linear pass region in the TCR domain, and verifies that it correctly recovers the corresponding quadratic or cubic phase signal. An expression for the noise power at the filter output is also presented in this section. Section 4 addresses discrete implementation of the proposed TCR filter and its computational cost. Section 5 presents simulations that illustrate this filtering scheme in separation or recovery of quadratic and cubic phase signals, including how to resolve back-to-back particles from in-line Fresnel holograms. Separation performance of the proposed filter is compared with those of Weyl and TF projection filters. Section 6 concludes the article.

Derivation of the proposed TCR filter
Let x(t) be an input signal involving amplitude modulated chirp (AM/FM) signal components and possibly a background noise component. Let s(t) = a(t) exp (j2π(t)) be the desired signal component with a narrow support region in the TCR plane that is nearly disjoint from those of other components and the noise in the input signal x(t). Then, can be viewed as the approximate TCR mask function of a separating filter. a denotes the chirp-rate (or frequency-rate) parameter, and, (2) (t) is the second derivative of the phase of s(t) yielding its ICR curve. In the above, we assume that the ICR curve of s(t) is correctly and accurately read on a TCR display of the input signal x(t) and is taken as the TCR domain mask of the filter.
Letŝ(t) be an estimate of s(t) obtained from the input signal x(t) aŝ where h(t, t') is the time-varying impulse response of the separating filter. Integrals are from −∞ to ∞ in this article unless indicated otherwise. Under the disjoint, narrow support region assumption, the optimal estimatê s(t) in the sense of minimizing mean-square error will be close to the projection of the input signal x(t) onto the phase signal exp(j2π(t)) [23,25]: Comparison of Equations (2) and (3) reveals in order to recover s(t). Hence, TCR domain filtering for signal separation consists of (i) displaying a TCR representation of the input signal x(t), (ii) selecting the TCR mask function of the filter so as to isolate the TCR signature of the desired component s(t) on this display, as in Equation (1), (iii) and obtaining the impulse response of the filter given in Equation (4) from the selected mask function in Equation (1). Then, obtaining Equation (4) from Equation (1) reduces to estimating the phase difference (t) − (t') associated with desired s(t) from the second derivative of its phase function (2) (t).
The Taylor series expansion of the phase (t) at the time point t' is By taking the first three terms in Equation (5) and substituting the trapezoidal approximation where t 0 is a reference time instant. We then seek a transform that maps a TCR mask function H(t, a) to the corresponding time-varying impulse response h(t, t'). When Equation (1) is substituted into this transform as the TCR mask, Equation (4) should be obtained as the impulse response with the exact phase difference replaced by its approximation given in Equation (7). Such a transform is given by where H f (t 0 , f ) accounts for an estimate of the IF value of the desired signal s(t) at the reference time t = t 0 , f (t 0 ) = '(t 0 ).
H f (t 0 , f ) in Equation (8) serves as a reference IF distribution around the given IF value. If it is taken as H f (t 0 , f ) = δ(f − '(t 0 )), then, substitution of it and Equation (1) into Equation (8) gives Equation (4) as the impulse response where the phase difference is replaced by its approximation given in Equation (7). This reference IF distribution of the desired component is indispensable in our proposed TCR filter given by Equations (8) and (2).
The proposed filtering procedure is given as follows: Step. 1. A TCR representation and a TF representation of the input signal x(t) are displayed as 2-D images.
Step. 2. The TCR mask function H(t, a) is prepared on the TCR display to isolate the ICR strip of the desired component s(t). This is idealized by Equation (1).
Step. 3. TF display of the input x(t) is examined and a convenient reference point (t 0 , f (t 0 )) is selected on the IF curve or in the TF support region of the desired s(t). Then, a reference IF distribution H f (t 0 , f) is prepared around the value f (t 0 ) at the reference time t 0 . This is idealized as Step. 4. The TCR mask H(t, a), its slice at t 0 , and the reference IF distribution H f (t 0 , f ) are substituted into Equation (8) to obtain the filter impulse response h(t, t').
Step. 5. Time-varying impulse response h(t, t') is applied to the input signal x(t) by Equation (2) to yield an estimate of the desired component s(t).
Higher order derivatives in Equation (5) could also be retained and approximated by differences of second derivatives evaluated at different time points. This leads to alternative forms of the TCR filter in place of Equation (8). For example, third derivative in Equation (5) can be approximated by a difference of second derivatives. The remaining terms can be discarded. Alternatively, the integral in Equation (6) can be approximated at three time points t 0 , (t 0 + t')/2 and t', instead of two. Both approaches lead to time varying impulse responses with four product terms in them.
These alternative filters can also successfully recover quadratic and cubic phase signals with monotonic ICR laws exhibiting single linear tracks, as the one proposed in Equation (8) does. This can be verified by showing that their equivalent Weyl TF transfer functions are also concentrated around IF curves of desired signals, as we show for the proposed TCR filter in the next section.
However, our simulations indicate that their performances in chirp signal recovery are worse than that of the proposed one, since their equivalent TF transfer functions exhibit more severe peaks near the origin of the TF plane. Moreover, their discrete implementations require more than one discrete TCR mask functions to be prepared and used, each for a different product term in the filter impulse response. This further complicates their discrete implementations. Our proposed filter in Equation (8) has the best separation performance and is easiest to implement, among them.

Equivalent Weyl TF transfer function
Equivalent Weyl TF transfer function of the proposed TCR filter in Equation (8) will be derived below for a linear TCR pass region approximated as a line impulse as given by Equation (1) and a rectangular pulse-shaped IF distribution at a reference time. For a more realistic case of a linear strip-shaped TCR pass region, we could not evaluate the resulting complicated integral to obtain an analytical expression for the TF transfer function.
We assume that the TCR mask is selected on a TCR display to follow the ICR curve of the desired component accurately and it is approximated by Equation (1). We also assume that the IF value of the desired component at time t = t 0 is read from an accompanying TF display. Then, substitution of Equation (1) and an initial impulsive approximation for the reference IF distribu- is assumed to be the cubic phase of the desired signal s(t) or that of the phase signal underlying the filtering operation. We have to verify that the TF transfer function of the filter given by Equations (9) and (10) is concentrated along the IF curve f = '(t) of s(t), in the TF plane, so that this filter will recover the desired s(t).
The Weyl TF transfer function of a linear, time-varying filter is given by [17,19,20] in terms of its impulse response. By substituting Equations (9) and (10) into Equation (11), we obtain which is the FT of the above CPF. It is concentrated around the IF curve f = 3at 2 +2bt+c of the desired signal, as required for its recovery.
The above integral can be expressed in terms of Bessel functions [53] and can be related to an Airy function [54] to roughly characterize its TF pass region along the IF curve: for f > 0, a > 0 or f < 0, a < 0. K 1/3 (·) and Ai(·) above denote the modified Bessel function of the second kind with order 1/3 and the related Airy function, respectively. For f > 0, a < 0 or f < 0, a > 0, we have where J 1/3 and J −1/3 are Bessel functions of the first kind with given orders.
By combining Equations (12), (13) and (14), we obtain (15) as the Weyl TF transfer function for the case of an impulsive reference IF distribution.
We now take a rectangular reference IF distribution: where rect(x) = 1 for |x| ≤ 1/2 and zero otherwise. When it is substituted into Equation (8), together with Equation (1), the filter impulse response becomes where sinc(x) = sin(πx)/(πx) and the phase term above is as given by Equation (9) together with Equation (10).
Substitution of Equation (16) into (11) yields that can be evaluated by convolving the right side of Equation (15) with the FT of B f sinc(B f τ), i.e., rect(f / B f ) in the frequency direction. Then, where the profile of the TF pass region of the filter around the IF f = 3at 2 + 2bt + c, at a fixed time, is obtained as with sign(a) denoting the sign of a.
The integral above has been evaluated by using [54] f where c 1 = 0.355028053887817, c 2 = 0.258819403792807, (20) and Time-frequency pass region profile H(f ) of the proposed TCR filter is plotted in Figure 1a,b for the scale factor taken as (4π/3) 2/3 /a 1/3 = 1 and −1, respectively. B f = 4 Hz and 62 terms are included in power series expansions given in Equation (21), for both cases. These plots reveal that the profile function is concentrated around f = 0; hence the Weyl TF transfer function H W (t, f ) given by Equations (18)-(21) is concentrated around the IF curve f = 3at 2 + 2bt + c correctly.
Time-frequency pass region can be determined from first zeros of H(f ) given by Equation (19) for a > 0, and for a < 0. Equations (22) and (23) determine the resolution limit of the TCR filter for separation of cubic phase signals with respect to slopes of their ICR lines in the TCR plane. Suppose that two such signals, s(t) = exp(j2πat 3 ) ands = (t) = exp(j2πāt 3 ), -T/2 ≤ t ≤ T/2, are to be resolved, where both a,ā > 0. If the TCR mask is selected to isolate a = 6at in the TCR plane to reconstruct the first signal s(t), then the segment of the second signals(t) around the point (t, 3āt 2 ) in TF plane can not be resolved from the desired s(t) provided that The slope range above is obtained from Equation (22) by substituting ' '(t) = 3at 2 and f =φ (t) = 3āt 2 into it. B f ≥ 2/T should be maintained above.
If a = 0, corresponding to a quadratic phase desired signal or reference signal onto which the input signal is projected in our filtering scheme, then Equation (12) reduces to the FT of the unity signal: a line impulse along the linear IF law of the desired quadratic phase s(t), in the case of an impulsive refer-  as the Weyl TF transfer function for this more realistic assumption.
Then, portions of a signals(t) = exp(j2πbt 2 ), −T/2 ≤ t ≤ T/2, around the point (t, 2bt) in TF plane can not be resolved from a desired signal following from Equation (26). B f ≥ 2/T should hold for recovery of s(t), for the signal duration −T/2 ≤ t ≤ T/2.
For higher order polynomial phase signals, where Equation (6) is neither exact nor a good approximation, equivalent TF transfer function of the proposed filter does not capture the correct IF curve. This is also the case for segmented quadratic or cubic phase signals for which Equation (6) is again not valid for the whole signal duration. Components of such a signal should be recovered one by one by repeated use of our filter with a different TCR mask each time. Exactness of the trapezoidal approximation in Equation (6) is the key to our proposed TCR filter. It is valid only for a single, linear pass region in the TCR plane corresponding to a quadratic or cubic phase signal with a monotone ICR curve.

Output noise power
If we take a uniform strip-shaped TCR pass region and a pulse-shaped reference IF distribution, then the TCR mask function and the IF distribution of the proposed filter are given as (28) respectively. Substitution of Equation (28) into Equation (8) gives the filter impulse response as where the phase term above is as given in Equation (9). Let x(t) = s(t) + w(t) be a noisy input signal for the proposed filter with the impulse response given in Equation (29), where s(t) is the desired signal component and w(t) is additive, zero-mean, white noise with power spectral density S w (f) = h.
An estimateŝ(t) of the desired component s(t) is obtained at the filter output as given by Equation (2). The noise component at the filter output, denoted as n (t), corrupting this estimate is given by The variance, i.e., the average power of the noise at the filter output can be obtained as The first line above is obtained by substituting Equation (30) in the expectation, and, the second line is obtained by substituting Equation (29) into the first line.
The output noise power given in Equation (31) does not depend on the phase, (t), of the desired component or that of the phase signal associated with the TCR pass region of the filter. It is determined by TCR and reference IF bandwidths, B a and B f , respectively, current and reference time values t and t 0 , and input noise power only, regardless of the phase being quadratic or cubic.

Discrete implementation
Discrete implementation of the proposed filter is described below. The three integrals in Equation (8) are discretized by considering time samples h(nT, mT) of the impulse response h(t, t') and those of its three components with a common sampling period T. Taking T = 1, samples of the first filter component can be written as The second expression above is obtained by writing the first one as a sum of integrals over intervals of length 2, changing the variable in each integral and using the fact that the discrete complex exponential above is periodic in the chirp-rate variable α with period 2. n 0 denotes the discrete reference time of the filter.
The term in square brackets in Equation (32) is sampled in the time variable and is periodic in the chirp-rate variable a with period 2. Hence, it can be viewed as the discrete-time TCR mask of the filter, denoted as H d (n, a).
Discrete version of the first integral in Equation (8) is, then, given by where H d (n, a), 0 ≤ n ≤ N − 1, 0 ≤ a ≤ 2 is the discrete-time TCR mask function of the filter, associated with the term in Equation (32). a denotes the discrete chirp-rate parameter above with an overuse of notation. N is the length of the discrete-time input signal x(n), 0 ≤ n ≤ N − 1, of the filter.
M-point Riemann sum approximation of Equation (33) gives the first filter component: Similarly, the second integral in Equation (8) can be discretized as Discrete version of the third integral in Equation (8) is obtained by similar steps, as with H f,d (n 0 , f ) being the discrete reference IF distribution. N-point Riemann sum for Equation (36) is computed via N-point inverse discrete Fourier transform (IDFT), as  (34), (35) and (37). Finally, the output signal of the filter, y(n), is computed by 4N, 8N, 16N, etc., give good results in the simulations presented below.
Computational cost: The proposed TCR filter can be efficiently implemented by means of inverse fast Fourier transform (IFFT) algorithms. Equation (34) can be evaluated by an M-point IFFT for each m value and by a subsequent index-finding among the stored values using the periodicity of the complex exponential kernel with period M. Equations (35) and (37) N + N log N) per output sample [17].
However, the main computational expense of our filtering scheme results from computing a TCR representation on which the filter TCR mask is selected. Such a TCR representation either has a quadratic phase kernel function, as the CPF [44,45], and thus can not be computed by fast algorithms or it requires interpolation by irrational factors, as in [43,49]. Hence, its computa-

Reconstruction of cubic phase signals
The use of the proposed TCR filter is illustrated for the following noisy multicomponent input signal with three cubic phase and one quadratic phase components: for 0 ≤ n ≤ N − 1, where the signal length is taken as N = 128 samples. w(n) denotes additive, zero-mean, white Gaussian noise with variance s 2 , above. The desired signal component to be estimated at the filter output is the first cubic phase component: 0 ≤ n ≤ N -1.

Noiseless input case
We first explain steps of the proposed TCR filtering scheme when there is no noise in Equation (39).
(i) TF display of the input signal: Figure 2a,b display spectrograms of the noiseless input signal in Equation (39) with Hann windows of widths 23 and 11 samples, respectively.
A quadratic IF curve that starts at zero frequency and increases with time can be identified on the left part in Figure 2a. That IF curve belongs to the first cubic phase component in Equation (39), which is also given in Equation (40) as the desired component to be reconstructed at the filter output. Two quadratic IF curves that start at frequencies of π and 2π radians and decrease with time can also be identified on the left part in Figure 2a. They belong to the second cubic phase component in Equation (39). The quadratic IF curve that starts at π/2 radians and increases belongs to the third component. The quadratic phase component is represented by the IF line starting at zero and ending at π/2 values in Figure 2a.
Each quadratic IF curve traverses the discretized TF plane several times, in this simulated signal scenario; hence, they become difficult to identify as time increases, on the right part of Figure 2a. Figure 2b reveals these quadratic curves more clearly, on the right part.
In particular, it is difficult to identify the IF curve of the desired component in Equation (40) and prepare a TF mask to isolate it. TF filtering is difficult to use for this signal separation task involving rapidly changing IF curves of undersampled signal components. However, proposed TCR filter can handle it more easily.
(ii) TCR display of the input signal: TCR patterns of the input signal x(n) are obtained by computing and displaying a bilinear TCR distribution derived in [49], which can also be viewed as a modified version of the CPF: Discrete radian chirprate range is taken to be [0, 4π) above to match that of the discrete TCR mask H d (n, 2k/M), 0 ≤ k ≤ M − 1, of the proposed filter, since the mask is prepared based on a display of the TCR distribution given in Equation (41). Thus, two periods of the discrete bilinear TCR distribution in the chirp-rate variable are computed and displayed. M = 8N is taken. Figure 2c-e show segments of the absolute value of the bilinear TCR distribution in Equation (41) computed for the input signal in radian chirp-rate ranges [0, π/2], [π/2, π] and [2π, 3π], respectively. Horizontal lines are ICR lines of the quadratic phase chirp in Equation (39), and, oblique lines with positive and negative slopes are ICR lines of cubic phase input components in Equation (39), in two periods of the modulus of the TCR distribution.
Instantaneous chirp-rate lines with a larger positive slope, in these figures, belong to the desired component given by Equation (40). Those with a smaller positive slope represent the third cubic phase component. ICR lines with a negative slope that start at chirp-rate values of π and 3π radians correspond to the second cubic phase component. Figure 2e shows all these ICR traces together, in the second period.
(iii) Preparing the TCR mask: Figure 2f shows a segment of the prepared filter TCR mask H d (n, 2k/M), for 0 ≤ n, k ≤ N − 1, corresponding to the radian chirp-rate range [0, π/2], that isolates the ICR line with the larger positive slope in Figure 2c belonging to the desired signal in Equation (40). This linear mask is chosen to be 1 sample wide vertically.
The replica of this line in Figure 2e, located in the range [2π, 4π], should not be included in the TCR mask, since its inclusion would result in an additional, undesired pass region in the TF plane for the equivalent TF transfer function, in addition to the desired TF pass region. If the second cubic phase component were desired at the filter output, then the ICR line with negative slope in Figure 2e, located in the range [2π, 4π], would be selected by the filter TCR mask, but its replica in Figure 2d located in the range [0, 2π] would be left out.
(iv) Selecting the reference IF distribution: Figure 2a indicates that the quadratic IF curve of the desired component in Equation (40) starts from zero frequency at time zero. Hence, we choose n 0 = 0 as the reference time point and H f,d (0, k/N) = 1 for 0 ≤ k ≤ 5, and zero otherwise, as the discrete reference IF distribution of the filter around the zero frequency value. The width of the distribution is determined by a search to maximize the separation performance.
(v) Computing the filter output: The filter output signal y(n) is computed from the reference IF distribution in part (iv) and the TCR mask in Figure 2f via Equations (34), (35), (37) and (38). Figure 3a displays the reassigned spectrogram of the filter output y(n) with a Gaussian window of width 9 samples in the form w(n) = exp(−n 2 /20), −4 ≤ n ≤ 4, showing only the quadratic IF curve of the desired component given in Equation (40). Figure 3b plots real part of this desired signal. Figure 3c plots real part of the output signal y(n) of the proposed TCR filter, after it is scaled by a number chosen to minimize the mean-square error between the desired and scaled output signals. Comparisons of Figure 2a with Figure  3a, and, Figure 3b with Figure 3c indicate that desired component is captured and reconstructed by this TCR filter.
The scale factor that minimizes the mean-square error, mentioned above, is calculated as where the desired signal s(n) is assumed to be known, as in Equation (40) for this simulation example. Its estimateŝ(t), real part of which is given in Figure 3c, is obtained from the filter output signal y(n) bŷ s(n) = βy(n).
Then, the normalized mean-square error (NMSE) between the desired and scaled output signals, used as a measure of the separation performance, is given as The percentage NMSE of our proposed TCR filter according to Equation (42) is 15.7655% for this noiseless simulation example, as the best result achieved by our filter. Best NMSE figures achieved by interpolated halfband Weyl filter [17] and full-band suboptimal TF projection filter [24] are 25.9042 and 20.8597%, respectively, for this case. They give worse results than our proposed filter for this simulation example, because these filters with TF masks that enclose IF traces of the desired component traversing TF plane eight times as shown in Figure 3a are no longer underspread systems. Their performances deteriorate when they become overspread [17,22], as in this case.

Noisy input case
The noise component w(n) is included in Equation (39) and its variance is varied to give signal-to-noise ratio (SNR) values listed in Table 1. The first cubic phase component given in Equation (40) is separated from the noisy input signal by the proposed TCR, Weyl and TF projection filters. Average percentage NMSE figures over 1,000 independent trials achieved by them are listed for these SNR values, in Table 1.
The proposed TCR filter gives significantly better results for this signal separation example except for low SNR values, since the desired signal and the filters designed to recover it are overspread. However, NMSE figures for low SNR cases indicate that the TCR filter is less robust against noise than others. At and below SNR = 0 dB, the TCR filter seems to break down, and is surpassed by the TF projection filter. This may be due to the first and second terms in the integral in Equation (31) that might have caused an effect of amplifying the output noise power especially pronounced when the reference IF bandwidth B f is kept larger to obtain a broader TF pass region for capturing the rapidly varying cubic phase component in Equation (40).

Resolution for cubic phase signals
To check validity of Equation (24) that indicates resolution limit of the proposed filter with respect to slope parameter a in the TCR plane, we take a noiseless input signal, where the first component is the desired signal s(n) given by Equation (40). When s(n) is passed through the TCR filter described above only, the NMSE is 1.9841%. Its slope parameter is a = 1/(48N).
When the signal above is given as input to the filter including the second component with the slope difference obtained from the increment in the upper bound in Equation (24) by substituting t = N = 128 there and by taking B f /2 = 6/N for the designed TCR filter, the NMSE becomes 30.4978%. Since the second component falls inside the TF pass region of the filter given by Equation (22), for the whole signal duration, the TCR filter could not separate it from the desired component.
When t = N/2 = 64 is taken in the slope difference Δ, the NMSE reduces to 18.2503%, since later parts of these two components are resolved in this case. These results suggest that Equation (24) accurately gives the resolution limit of the proposed filter in the slope parameter.

Recovery of a more slowly varying component
Suppose that the first cubic phase component, given explicitly in Equation (40), is replaced by the following half-band cubic phase signal in the noisy input signal given by Equation (39): Table 1 Percentage NMSEs in recovery of a cubic phase component from a noisy input signal: Percentage NMSEs averaged over 1,000 independent trials that are achieved by the proposed TCR, interpolated half-band Weyl and fullband suboptimal TF projection filters in recovery of the cubic phase component in Equation (40 0 ≤ n ≤ N − 1. The above signal is now separated from the noisy input signal by the proposed TCR, Weyl and TF projection filters for this case where all filters become underspread systems designed to recover this more slowly varying desired component. Best average percentage NMSE figures over 1,000 independent trials achieved by them are listed versus given input SNR values, in Table 2. The TCR mask is prepared to be 1 sample wide along the ICR line of the desired component given by Equation (43), with the number of chirp-rate bins taken as M = 512N for this example to fully expose this ICR line on the underlying TCR pattern, and H f,d (0, k/N) = 1 for k = 0, 1, N − 1, and zero otherwise, as the discrete reference IF distribution of the TCR filter around the zero frequency value, to yield best results achieved by the TCR filter as listed in Table 2.
For high and medium input SNR values, NMSE figures of the TCR filter are slightly worse than those of the Weyl filter, and slightly better than those of the TF projection filter, as seen in Table 2. But differences are not significant since all filters are underspread systems for this case of a more slowly varying desired component.
At and below SNR = 0 dB, results of the TCR filter deteriorate again due to output noise amplification caused by the first two terms in the integral of Equation (31). However, this deterioration, i.e., the increase in percentage error as input SNR is increased, is about half of the increase observed in Table 1 for the rapidly varying desired signal in Equation (40). This is because the reference IF bandwidth B f is chosen to be 3 samples wide for the desired signal in Equation (43), whereas it was 6 samples wide for the one in Equation (40).

Noiseless input case
To recover the quadratic phase component in Equation (39), the TCR mask of the proposed filter is selected to isolate the horizontal ICR line shown in Figure 2c, in the radian chirp-rate range [0, 2π], corresponding to this signal component. The replica of this line in Figure  2e, in the range [2π, 4π], is again left out of the TCR mask. The mask is chosen to be 1 sample wide vertically for the best separation performance. Figure 2a indicates that the linear IF curve of the desired component also starts from zero frequency at time zero, hence we choose n 0 = 0 as the reference time point again, and H f,d (0, k/N) = 1 for k = 0 and zero otherwise, as the discrete reference IF distribution of the filter as determined by a search to achieve the maximum separation performance. Figure 3d displays the alias-free WD (AF-WD) [55] of the output signal of our filter with the above TCR mask and reference IF distribution, demonstrating that the quadratic phase component is separated by the proposed filter. Figure 3e,f plot real parts of the desired quadratic phase component in Equation (39) and the normalized filter output, respectively.
The percentage NMSE of the proposed TCR filter for this separation problem, defined by Equation (42), is 0.2793%, whereas best NMSE results achieved by interpolated half-band Weyl filter and full-band suboptimal TF projection filter are 0.2795 and 0.9992%, respectively.
A noteworthy observation is that the filter TCR mask should be kept as narrow as possible for a good separation performance. While recovering the quadratic phase component in Equation (39), if a neighboring chirp-rate bin is also included in the horizontal TCR mask making it two samples wide, then the NMSE falls down to 14.5101%. Inclusion of the other neighbor for a threesample TCR mask results in a NMSE value of 17.6307%.
Hence, the chirp-rate value of the desired quadratic phase signal component should be determined accurately from the underlying TCR representation of the input signal, and, the filter TCR mask should match this value in order to recover that component with high accuracy. Quadratic phase signals with chirp-rates differing less than 4π/M radians are not resolvable by the filter.

Noisy input case
When the input signal given in Equation (39) is noisy, percentage NMSE values averaged over 1,000 independent trials that are achieved by the proposed TCR, Weyl and TF projection filters in recovery of the quadratic phase component are listed in Table 3 for various SNR values. The proposed TCR filter gives slightly better results for this example except for SNR = 5, 0 and −3 dB cases, but differences are not significant since all filters are underspread for this case. TF projection filter gives slightly worse results here only because its pass region could not be chosen as narrow as those of others, since it could not capture the desired component for such a narrow pass region due to implementation matters.
The TCR filter does not break down for low SNR cases here, since the reference IF bandwidth B f is only one sample wide in this case, limiting noise leakage.

Recovery of back-to-back objects from in-line Fresnel holograms
Formation of in-line Fresnel holograms can be modelled by a linear, shift-invariant system with a quadratic phase impulse response [56]. Hence, such a hologram can be roughly viewed as a linear combination of chirp signals centered at object plane locations of objects encoded in it.
For a discretized hologram, the normalized depth parameterᾱ of an encoded object is defined by [56] where X is the spatial sampling period in hologram plane coordinates, N is the discrete hologram image size in samples taken to be N = 256 in our presented simulations, z is the actual object depth, and l is the illumination wavelength. Figure 4a shows a simulated hologram of two back-toback objects having one-dimensional (1-D) variation like a long thin wire. The objects have the same object plane coordinates, but are located at different normalized depths α = 1 and 1/ √ 2, respectively. Horizontal profile of each object, when viewed as a 1-D discrete signal, is a binary rectangular pulse five samples wide. In Figure 4b, a row of this hologram is plotted after eliminating its DC level. Figure 4c displays the AF-WD of the signal in Figure  4b, as a SF pattern. The horizontal axis is the discrete space axis; x = 0, 1, ..., 255. The vertical axis is the discrete spatial radian frequency axis; w m = (2π/256)m, m [−128, 127]. Each object is represented by a pair of lines with opposite slopes that intersect on the spatial axis at the object coordinate. Magnitudes of the slopes are inversely proportional with real object depths. As can be seen in the figure, these pairs of linear tracks associated with the two objects overlap in the SF plane, making it tedious to design SF filter mask functions to resolve them [14]. However, objects are represented by distinct horizontal strips corresponding to different slopes or depths in the SCR plane. Hence, it is easier to design SCR mask functions isolating such strips and to use our proposed filter to separate and recover these overlapping objects.
The underlying SCR distribution for the proposed filtering scheme is chosen to be the CPF [44,45]: N denotes the input signal length again, and, M = 16N is chosen for our SCR patterns in this subsection. Discrete radian chirprate range is, again, taken to be [0, 4π) above, encompassing two periods, to match that of discrete SCR masks of the proposed filter that will be used in hologram component separation. Figure 4d-f display modulus of the CPF for the hologram signal in Figure 4b in chirp-rate ranges taken near the zero value, around 2π, and near 4π, respectively. Longer horizontal lines in these figures correspond to slope values of the line pair in Figure 4c with smaller slope magnitude, for the object with normalized depth α = 1/ √ 2. Shorter, less visible horizontal lines located symmetrically above or below longer ones correspond to the object with normalized depthᾱ = 1. We want to recover the object pulse located at the normalized depth α = 1/ √ 2 by separating its hologram component from the two-component hologram signal in Figure 4b by our proposed filter. For this purpose, filter SCR masks should be designed to isolate the longer ICR lines in these SCR patterns. Figure 4c indicates that IF lines of the desired object pass through the zero-frequency level at the intersection point with spatial coordinate n 0 = 101. Hence, reference spatial point for our filter is chosen as n 0 = 101, and the discrete reference IF distribution for our filtering operations is taken as H f,d (101, k/N) = 1 for 0 ≤ k ≤ 3 and N − 4 ≤ k ≤ N − 1 around the zero frequency value, and zero elsewhere.  Figure 5a shows a segment of the first SCR mask that selects the longer line in Figure 4d close to the zero chirp-rate value. Equivalent Zadeh SF transfer function [18] of the proposed filter, with this SCR mask and the reference IF distribution described above, is shown in Figure 5b. It isolates the IF line in  Figure 4b is, firstly, passed through this SCR filter. Figure 5c shows a segment of the second SCR mask that selects the longer line in Figure 4f close to the chirp-rate value of 4π. Equivalent Zadeh transfer function of the proposed filter, with this SCR mask and the same reference IF distribution, is shown in Figure 5d. It isolates the IF line in Figure 4c with the larger negative slope. The input hologram signal in Figure 4b is, then, passed through this SCR filter.
The hologram component encoding only the desired object is estimated as the sum of these two filter output signals obtained with different SCR masks given in Figure 5a,c, respectively. Real part of the sum signal is plotted in Figure 5e, as this separated component. The desired object pulse with normalized depthᾱ = 1/ √ 2 is finally estimated by applying a Fourier synthesis algorithm [14] to the separated hologram component. The recovered object is plotted in Figure 5f. Figure 6a,b show two simulated, planar objects overlapping in the object plane. Figure 6c displays their simulated hologram, where square shaped and triangular objects are located at normalized depthsᾱ = 1/ √ 2 and 1, respectively. A row that passes through the diffraction pattern in this discrete hologram image is selected and its DC value is eliminated. Figure 6d shows the AF-WD of the DC-levelled hologram row signal, as its SF pattern.
In Figure 6d, two back-to-back objects are represented by two pairs of lines with opposite slopes that intersect on the spatial axis (zero-frequency axis) at the common object coordinate n 0 = 130. Hence, this intersection point is chosen as the reference spatial point for our filter and its reference IF distribution at this point is, again, chosen as a rectangular pulse centered at the zero-frequency value: H f,d (130, k/N) = 1 for 0 ≤ k ≤ 7 and N − 8 ≤ k ≤ N − 1, and zero elsewhere.
We use a shifted version of the local quadratic periodogram (LQP) [50][51][52] to obtain SCR patterns for the hologram row: 16N again, and the shift amount is taken to be equal to the spatial coordinate of the intersection point in Figure 6d, n 0 = 130, since input chirp components are centered at this point. Discrete radian chirp-rate range is, this time, taken to be [0, 2π), corresponding to a single period. By this way, the computed transform peaks at chirp-rate index values matching the correct chirp-rate values of the input components, which are twice the coefficient values of the quadratic phases of input components. Figure 6e,f display modulus of the shifted LQP in Equation (45) for the DC-levelled hologram row, with a rectangular window of width 141 samples, displayed in chirp-rate ranges taken near the zero value and 2π, respectively. The lower smeared line in Figure 6e and the upper smeared one in Figure 6f correspond to the square shaped object. The upper horizontal line in Figure 6e and the lower one in Figure 6f are ICR lines corresponding to the triangular object, in the SCR plane. Figure 7a,b show filter SCR masks that isolate horizontal traces in Figure 6e,f, respectively, for the square shaped object. The hologram image in Figure 6c is DC levelled, then, its rows are passed through the proposed filter with the described IF distribution and these different SCR masks, separately. The two filter output images are summed, and then a Fourier synthesis algorithm [14] is applied to the resulting image, to recover the square shaped object as shown in Figure 7c.
Similarly, Figure 7d,e show filter SCR masks that isolate horizontal traces in Figure 6e,f, respectively, for the triangular object. Rows of the DC-levelled hologram image are passed through the proposed filter with the described IF distribution and these different SCR masks given in Figure 7d,e, separately. The two filter output images obtained with these two masks are summed, and the triangular object is recovered from the sum image, as shown in Figure 7f.
Object reconstruction procedure described in this subsection does not require prior knowledge of object depths, or a manual search for them, in order to recover the encoded objects from their holograms. It automates this task.

Conclusion
We propose a novel linear time-varying filtering scheme in the joint TCR domain for decomposition of multicomponent signals into their quadratic and/or cubic phase components, in this article. It is valid only when the approximation involved in Equation (6) becomes exact, i.e., it can only separate quadratic or cubic phase signal components with single linear pass regions in the TCR plane corresponding to monotonic ICR laws. This is verified by deriving approximate expressions for equivalent Weyl TF transfer functions of the proposed filter and by checking whether they contain IF curves of such signal components or not, in Section 3.
Multicomponent quadratic or cubic phase signals can be recovered by repeated use of our filter with a different TCR mask, each time, designed to separate a particular signal component. Simulation results presented in Table 1 suggest that our filter decomposes (e) (f) Figure 7 Reconstruction of back-to-back objects from a two-dimensional hologram: (a) and (b) Filter SCR masks that select horizontal traces in Figure 6e,f, respectively, to recover the square shaped object in Figure 6a. (c) Recovered object to compare with the original one in Figure 6a. (d) and (e) SCR masks that select horizontal traces in Figure 6e,f, respectively, to recover the triangular object in Figure 6b. (f) Recovered triangular object for comparison with Figure 6b.
Özgen EURASIP Journal on Advances in Signal Processing 2012, 2012:122 http://asp.eurasipjournals.com/content/2012/1/122 multicomponent quadratic and cubic phase chirps with rapidly varying IF curves with significantly improved accuracies over Weyl and TF projection filters for high and medium input SNR values, albeit with increased computational cost. The reason for this improvement is that filters designed to capture such signals can no longer be underspread systems and hence, Weyl and TF projection filters can not perform so well in recovery of such signals. However, for low SNR cases, the proposed TCR filter is less robust against noise in recovery of rapidly varying signals since broader TF pass regions required to capture such signals leak more noise and the first and second terms in the integral in Equation (31) further amplify the output noise power.
Average percentage errors listed in Tables 2 and 3 show that the proposed TCR, Weyl and TF projection filters give comparable accuracies in separation of more slowly varying signal components, at least for high and medium input SNR values, since all filters are underspread in such cases.
The idea of filtering in the joint TCR domain is novel. However, since the proposed filter is linear with a timevarying impulse response; it has an equivalent TF transfer function and it can also be designed in the TF domain in theory. The proposed TCR filter intends to ease the filter design task for some applications where it is more convenient to design the filter mask function in the TCR domain than in the TF domain; such as reconstruction of back-to-back objects from in-line holograms for which SF masks have to overlap, whereas SCR masks do not.
Another application involves separation and recovery of undersampled chirps with rapidly varying IF curves that are difficult to identify and enclose in the TF plane. If ICR lines vary more slowly and hence are more readable in the TCR plane, then it is more convenient to perform the filtering operation in the TCR plane by our proposed filter. For example, received radar signals from rapidly maneuvering targets can be modeled as multicomponent cubic phase signals [57]. They can be decomposed by the proposed TCR filter first, for improved estimation of phase parameters of each component to obtain focused SAR or ISAR images of targets as presented in [57].