An Alternative Derivation of the Method of Stationary ... - IEEE Xplore
2710
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
An Alternative Derivation of the Method of Stationary Phase for Multivariate Signals Otmar Loffeld, Senior Member, IEEE
Abstract—The method of stationary phase (msp) is a commonly established technique to asymptotically determine integrals (such as Fourier integrals) with strongly oscillating integrands. Basically, the method states that an integral over a strongly oscillating integrand only delivers contributions at those (“stationary”) points or regions, where the integrand does not change, i.e., where the phase derivative of its oscillation vanishes. In the engineering literature, various approaches for the derivation exist for univariate (1-D) signals, which, however, find their limitations for multivariate (multidimensional) signals, which cannot be factorized. On the other hand, from the field of applied mathematics, there are derivations, formulating and solving the problem in a rigorous and sophisticated way. This paper tries to close what the author interprets as a gap between these two worlds, aiming to substitute intuitive insights by some basic Dirac impulse modeling concepts, which are familiar to signal processing engineers. We first start with the general multivariate (multidimensional) case, introducing vectorial time and frequency representations. Here, we replace 1-D Taylor-series-based quadratic phase signals by multivariate signals with phase histories, being given by more general quadratic forms. The exponential, containing the quadratic form, is then related to an n-variate Dirac impulse, which provides the sifting property needed to establish the basic msp property that the integration over the entire integration range yields only one value of the integrand, evaluated in the point of stationary phase. This sifting capacity in connection with (time limiting) window functions in the integrand migrates time limitations of the signals to frequency limitations of the spectra—a multiplicative time window is transformed in a multiplicative frequency window of the same form, rather than into a convolution of the Fourier transforms. In this respect, this paper provides a closed-form comprehensive treatment of the subject. Index Terms—Chirp signals, monostatic point target reference spectrum, multidimensional method of stationary phase (msp), n-variate msp, quadratic phase history.
I. I NTRODUCTION
T
HE calculation of Fourier spectra S(f ) of multivariate signals s(t), where t = [t1 , t2 , . . . , tn ]T is some vector, denoting the time, and/or space dependence of a
Manuscript received August 17, 2010; revised December 20, 2010, March 7, 2011, and July 29, 2011; accepted October 15, 2011. Date of publication December 23, 2011; date of current version June 20, 2012. This work was supported in part by the German Science Foundation (DFG) under Grant DFG Lo 455/7-1, 7-2. The author is with the Center for Sensor Systems (ZESS), University of Siegen, 57068 Siegen, Germany (e-mail: [email protected]). Digital Object Identifier 10.1109/TGRS.2011.2174245
scalar-valued signal s(·) involves the solution of integrals of the form ∞ S(f ) = −∞ ∞
=
s(t) · e−j2πf
T
·t
dt
|s(t)| · ejφs (t) · e−j2πf
T
·t
dt
−∞ ∞
a(t) · ejφ(t,f ) · dt.
=
(1)
−∞
In those cases, where the solution of such integrals is not analytically possible, asymptotic solutions, like the method of stationary phase (msp), come into the focus of interest. This method, which states that the stationary points give the main contributions to such integrals, has been and is extensively used in, e.g., monostatic and bistatic synthetic aperture radar (SAR) systems [1]–[6], geophysics, antenna theory, and wave theory [7], [8] or seismics when multidimensional formulations are included. a(t) is some envelope or window function, limiting the time extent of the integrand in a way that a(t) sufficiently fast decays for large values of t. This function may even be discontinuous at the ends of the interval of integration but is supposed to contain neither singularities nor heavy oscillations. The phasor function φ(t, f ) is supposed to possess time derivatives up to second order and a nonvanishing and finite secondorder derivative in the interval of integration, particularly at the point of stationary phase, e.g., where the first-order derivative with respect to time vanishes. This phasor introduces an oscillating behavior (depending on f ) of the integrand, and the msp basically states that the integral only delivers a contribution at those positions t, where the phase is stationary, e.g., where its first-order derivative with respect to time t vanishes. In ray or diffraction theory, such points of stationary phase (stationary points) are often called critical points of the first kind [9]. They correspond to rays leaving the reflector normal to the tangent plane. In addition, there are critical points of the second and third kinds, due to diffraction at an edge and a corner, respectively, where cancellations in the exponential phase factor may become ineffective in the neighborhood of the limits of integration. They are considered in [10] and [11] in connection with the evaluation of infinite integrals, arising in optical diffraction problems. Readers are also referred to [12] for a further discussion. While the msp is readily formulated for univariate signals s(t), the extension to the general case,
0196-2892/$26.00 © 2011 IEEE
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
particularly where the signal s(t) cannot be factorized in a product of 1-D signals s(t) = ni=1 Si (ti ), needs some special consideration. In the multivariate case, this point of stationary phase becomes multidimensional and, due to the dependence of the phase on the integration variable t, may lead to cases where the multidimensional integral cannot be decoupled into individual integrals with a decoupled determination of the individual points of stationary phase for each integral. On the one hand, simple (and possibly too simple) extensions from the 1-D (univariate) case to the 2-D case have been published in [13] and [14]; on the other hand, this method has received ample treatment in the mathematical literature, and several rigorous expositions are available, e.g., [15], [16], and references therein. The first citation treats the subject using integration by parts, and the second one treats the subject using generalized Dirac functions and their derivatives. A multiple integral is reduced to a single line integral, using Dirac line functions, and then solved, where a proof, contained in the latter reference, has been further simplified in [17]. The advantage of such mathematic rigor is a precise formulation and understanding of preconditions and constraints, under which the approximation holds, as well as error bounds of the approximation. On the other hand, these works do not seem common reading, neither in signal processing nor in geoscience and remote sensing. The multidimensional stationary phase method has also appeared in book form, for example, [18] and [19], treating the subject in a comprehensive and thorough (excellent) way, but here, the price is paid in form of the length. This paper seeks a compromise; it substitutes an oversimplistic [13], [14] approach by an engineering mathematical approach, without, however, claiming the mathematic rigor in [15] and [16]. We start with the general multivariate (multidimensional) case, introducing vectorial time and frequency representations and replacing 1-D Taylor-series-based quadratic phase signals by multivariate signals with phase histories being given by more general quadratic forms. Summarizing the problem, we then look at integrals of the form
spectrum of an infinitely extended chirp and compare it with the analytical solution, which can be obtained. The second example considers the determination of the monostatic SAR point target reference spectrum, for which the results are commonly known, and the third example looks at an arbitrary chirp with some general quadratic form phase history and some other elliptical or hyperbolic envelope, to demonstrate the potential of the framework formulation. II. n-VARIATE MSP Recalling (2), we want to obtain the solution of integrals of the form ∞ A(x) · ejΦ(x) dx
I= −∞
=
···
A(x1 , . . . , xn ) · ejΦ(x1 ,...,xn ) dx1 , . . . , dxn (3)
n
where ejΦ(x) constitutes a scalar and quickly oscillating term with a large second-order derivative matrix (with all eigenvalues being sufficiently large). Furthermore, we assume that the assumptions concerning the window function A(x) and phase argument Φ(x), as discussed in the introduction, are met. Using the results in Appendix A, we suppose in the first step that we may approximate with sufficient exactness 1 x) + (x − x ˜)T · P0 (˜ x) · (x − x ˜) Φ(x) = Φ(˜ 2
A(x) · ejΦ(x) dx.
(2)
−∞
The dependence of the integral and the phasor φ(t, f ) on the parameter f has been omitted for simplicity. In the next section, the phase function Φ(x) will be approximated by a secondorder Taylor series approximation, expressed by a quadratic form. Such exponentials, containing quadratic forms, constitute multivariate Gaussians (cf. [20] for the scalar case). In the Appendix, we further show how the multivariate Dirac impulses can be interpreted as limiting cases of Gaussians, bringing a Dirac function into our integral. Conceptually, these Dirac functions, used in the integrand and being only nonzero, where their argument vanishes, give rise to line integration concepts, as discussed in [12]. The idea of expressing the implicit msp sifting property by Dirac impulses has already been used in [21] but here with the restriction to scalar signals. Having derived the formula framework, we apply this technique in the final chapter to three different applications. First, we calculate the n-variate
(4)
x ˜ is the point of stationary phase defined in Appendix A, given x) = (d/dx)[[(d/dx)Φ(x)]T ]x˜ by x ˜ : (dΦ(x)/dx)|x˜ = 0T . P0 (˜ is the second-order derivative of the multivariate phase history in the point of stationary phase, defined in Appendix A. We further assume that all eigenvalues of the matrix P0 are sufficiently large. Then, we have from (4) I∼ = I˜ =
∞ I=
2711
∞
1
A(x) · eΦ(˜x) · ej 2 (x−˜x) −∞ ∞
1
A(x) · ej 2 (x−˜x)
= eΦ(˜x) −∞
T
T
·P0 (˜ x)·(x−˜ x)
·P0 (˜ x)·(x−˜ x)
dx
dx .
(5)
I1
Looking at the integral I1 and then substituting x − x ˜ = u, where du = dx, we may write ∞
1
A(x) · ej 2 (x−˜x)
I1 = −∞ ∞
1
T
A(u + x ˜ ) · ej 2 u
=
·P0 (˜ x)·(x−˜ x)
T
·P0 (˜ x)·u
dx
du.
(6)
−∞
We now introduce an invertible orthogonal transformation of the form v = T · u => u = S · v where S −1 = T and T T = T −1 .
(7)
2712
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
Then, we have for the integration volume element
Substituting (15) into (5) and further using (3), the multivariate msp becomes
dv = |det(T )| · du
(8)
∞
−1
T
and due to T = T , we have det(T ) = 1. Substituting this into (6), we find that ∞ A(S · v + x ˜) · e
I1 =
j 12 (S·v)T ·P0 (˜ x)·S·v
−∞
dv
i=1
1
T
·S T ·P0 (˜ x)·S·v
1
T
·D(˜ x)·v
dv
−∞
∞
A(S · v + x ˜ ) · ej 2 v
dv.
(9)
−∞
Now, we require the matrix in the exponential to take on a diagonal form S T · P0 (˜ x) · S = D(˜ x) = diag (dii (˜ x)) .
(10)
i=1,n
It can be shown that the solution of (10) is a column vector matrix S, where the column vectors si , i = 1, . . . , n, obey the following equation: x) · si = dii (˜ x) · si , P0 (˜
i = 1, . . . , n
(11)
x) and si denotes saying that dii denotes the eigenvalues of P0 (˜ the eigenvectors corresponding to dii . Furthermore, due to the orthogonality of T and S, we also have: det(D(˜ x)) = x)). For the integral (9), we have then det(P0 (˜ ∞ I1 = −∞
1 T A(S · v + x ˜) · exp j v · D(˜ x) · v dv. 2
(12)
In Appendix B, we show that, provided that the eigenvalues of the diagonal matrix D are sufficiently large, the complex exponential with the quadratic form, given in (12), can be approximated by a multivariate Dirac impulse. Using this equivalence as specified in (57) in Appendix B, we obtain due to the diagonal form of D ∞
n
−∞ n n ∼ = (2π) 2 · (j) 2 ·
1 1 2
1
n
A(˜ x + S · v) · (2π) 2 · (j) 2 ·
I1 =
(16)
1
A(S · v + x ˜ ) · ej 2 v
=
−1
dii 2 (˜ x)
=det− 2 (P0 (˜ x)).
∞
!
n
n ∼ x) · ejΦ(˜x) · (j2π) 2 · = A(˜
−∞
=
A(x) · ejΦ(x) dx
I=
det D(˜ x)
1 2
det D(˜ x)
· δ(v)dv
· A(˜ x)
(13)
x ˜ : (dΦ(x)/dx)|x˜ = 0T is the multidimensional point of stationary phase, defined in (44), which, in the case of a Fourier transformation integral, depends on the frequency. x) denotes the eigenvalues of P0 (˜ x), where P0 (˜ x) = dii (˜ (d/dx)[[(d/dx)Φ(x)]T ]x˜ ) is the second-order derivative of the multivariate phase history in the point of stationary phase defined in (34). Two major items are worth noting in (16): First, we note that the product term of eigenvalues is equivalent to the x). Second, we note that rather than having determinant of P0 (˜ to deal with a convolution of Fourier transforms, we have a multiplication of the window function A(˜ x) with the phase history term, evaluated in the point of stationary phase. This essentially says that a multiplication with a strongly oscillating phase term carries over into the multiplication of the same parts in the frequency domain, indicating that multiplications do not necessarily always transform in convolutions in the frequency domain, if strongly oscillating terms are involved. This feature is implicitly introduced by the sifting property of the msp integral, being expressed by the Dirac impulse formalism in our derivation. It is essentially the deeper understanding of this feature which constitutes the implicit worth of the formal derivation. In the next part, we will look at application examples, further demonstrating the worth of this closed-form solution. III. E XAMPLES AND A PPLICATIONS A. Spectrum of the Infinitely Extended Chirp The first example is meant to provide a simple check of reasonableness with some infinitely extended chirp signal s(t) of elliptical geometry, for which its Fourier transform can be analytically determined. In order to ensure bounded total variation and, thus, existence of the involved Fourier integrals, we first assume s(t) = exp{jπ·tT ·Q·t} ⇒ S(f ) = F exp{jπ·tT ·Q·t} (17) where Q = Q1 + jQ2
where 1 2
D(˜ x) = diag (dii (˜ x)) → det (D(˜ x)) = i=1,n
n
1
dii (˜ x) 2 .
(14)
i=1
Hence, we obtain for the integral (13) n I1 ∼ = (j2π) 2 ·
n
i=1
−1
dii 2 (˜ x) · A(˜ x).
(15)
with Q1 being real, symmetric, and positive definite, and Q2 being real, symmetric, and nonnegative definite, with the following symmetries: QH = QT1 − jQT2 = Q1 − jQ2 = Q∗ . Conceptually, the matrix Q2 in the exponential ensures an exponential decay of the chirp envelope. Later on, we will obtain the spectrum of the infinite chirp as the limiting case, letting Q2 (more precisely its eigenvalues) go to zero.
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
Straightforward Solution: Rather than trying to evaluate the Fourier integral, we derive a differential equation describing the signal, which will then be transformed into the frequency domain, showing that the frequency spectrum basically obeys the same differential equation and thus will be again a chirp d s(t) = jπ·exp{jπ·tT ·Q·t}·tT ·Q = j2π·s(t)·tT ·Q. dt
Then, the corresponding phase at the stationary point is obtained as T φ(t˜) = − [2πf T − π t˜ Q1 ] · t˜ −1 = − [2πf T − πf T Q−1 1 Q1 ] · Q1 ·f = − πf T Q−1 1 f.
(18)
Transforming (18) into the frequency domain, we obtain, using (d/dt)s(t) = F −1 {j2πf T · S(f )} (19) j2πf T · S(f ) = F F j2π · s(t) · tT · Q .
2713
(26)
The eigenvalues of the matrix Q1 are given by 2 q11 + q22 q11 · q22 − q12 · 1± 1−4· (27) λ1,2 = 2 [q11 + q22 ]2 where
Now, using the dual relationship F {−j2π · s(t)} = (d/df )S(f ) for the spectrum, we obtain from (19) j2πf T · S(f ) = −
d S(f ) · Q df
d S(f ) = − j2π · S(f )·f T · Q−1 df
(20)
which is exactly the same differential equation as (18), except for the minus sign and the inverse matrix Q−1 . Hence, the wanted spectrum must be given by S(f ) = C · exp{−jπ · f T · Q−1 · f }
(21)
where the constant C can be determined by ∞ S(0) = C =
∞
n
Now, substituting the results of (25), (26), and (27) into (16), we obtain the final result for the spectrum −1 S(f ) = j2π · (2π)2 λ1 · λ2 · exp −π · f T · Q−1T Q2 · Q−1 1 1 f (28) · exp −jπ · f T · Q−1 1 ·f . Now, letting Q2 go to zero, we obtain the final result of (23) −1 S(f ) = j2π · (2π)2 λ1 · λ2 · exp −jπ · f T · Q−1 1 ·f 1 · exp −jπ · f T · Q−1 =j 1 · f . (29) 2 q11 · q22 − q12 C
−∞
B. Monostatic SAR Point Target Reference Spectrum
− 12
= j 2 · det(Q).
(22)
For n = 2, and letting Q2 go to zero, we obtain the spectrum of the infinite chirp S(f ) = j · det(Q1 ) · exp −jπ · f T · Q−1 1 ·f . − 12
(23)
We now determine the spectrum using the msp, given in (16) ∧ and t = x; A(x) = exp[−πxT · Q2 · x]. We have
φ(t) = − π 2f T · t − tT Q1 · t d φ(t) = − 2π[f T − tT Q1 ] dt T d d P0 (t) = φ(t) = 2π · Q1 . dt dt
(24)
The point of stationary phase t˜ is given by d ˜ T φ(t) = 0T ⇒ −2π[f T − t˜ Q1 ] = 0T dt −1 ˜ ⇒ t˜ = f T · Q−1 1 ⇒ t = Q1 ·f . T
d22 = 2πλ2
exp{jπ · tT · Q · t}dt
s(t)dt = −∞
d11 = 2πλ1
(25)
Our next application is the determination of the 2-D point target reference spectrum of a monostatic SAR sensor. This is an example, where the outcoming result is well known and well reported in the literature. Conventionally, the monostatic point target spectrum is calculated by first transforming the point target response into the range frequency domain and then subsequently into the azimuth frequency (Doppler) domain, using the 1-D msp. Furthermore, implicitly present in the classical understanding, the derivation must be done exactly in this sequence of 1-D transformations, since the other way round (first, an azimuth Fourier transformation, followed by the range Fourier transformation) would not work, since no stationary point could be identified. Our example will show, however, that, by addressing the problem in a 2-D way, we implicitly identify a 2-D stationary point. In doing so, the sequence of transform operations is not an issue anymore. The whole problem reduces to the solution of a system of equations, where we have two equations for the determination of a 2-D point of stationary phase. Considering the prerequisites of applying the msp (requiring a second-order derivative matrix with sufficiently large eigenvalues), there are conditions formulated for the 1-D case, such as in [22, pp. 272] or in the analysis given in [23] which, however, still need to be reformulated for the multivariate case. In some publications, like [24], this is reflected by requiring signals with sufficiently large time-bandwidth product, which
2714
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
seems to be necessary but not completely sufficient condition. We omit a deeper analysis here for simplicity and space. For the given example, the sensor moves along the azimuth axis, denoted by azimuth time τ , the fast time (range time) coordinate is denoted by t, the sensor’s azimuth velocity is v, and the point target position is described by the azimuth time instant of closest approach τ0 and by the (orthogonal) slant range vector R0 at that time, where R0 = R(τ, R0 , τ0 )|τ =τ0 = R02 + ν 2 (τ − τ02 )τ =τ0 . After quadrature down conversion to zero offset frequency, the point target response in the range azimuth time domain is then given by t − t0 (τ ) sT (t, τ ) = σ(R0 , τ0 ) · rect TR · exp jπkr (t − t0 (τ ))2 τ − τc · rect · exp (−j2πf0 · t0 (τ )) . (30) Taz kr is the chirp rate of the transmitted signal, and TR and Taz are the time extensions, specifying the time apertures of the signal, which are described by rectangular windows. τc is the azimuth time, when the point target is seen in the center of the (azimuth) antenna footprint (due to a squinted geometry), f0 is the carrier frequency, σ(R0 , τ0 ) describes the radar brightness of the point target, and the two-way propagation delay is given by 2 2 t0 (τ ) = R(τ, R0 , τ0 ) = R02 + v 2 (τ − τ02 ). (31) c c Implicitly, the point target response would be 4-D, since the signal clearly depends on the location of the point target, but this dependence is neglected for convenience. From (30), the Fourier spectrum of the signal normalized to its radar brightness is obtained by ST (f, fτ ) = ST (f ) ∞ ∞ t − t0 (τ ) τ − τc rect · rect = TR Taz −∞ −∞
∞
· exp [jφ(t, τ )] dtdτ
w(t) · exp [jφ(t)] dt
=
kr is the range chirp rate of the transmitted signal. In order to identify the point t˜ of stationary phase (16), the two partial derivatives of (33) must vanish d ! φ(t)|t˜ = − π · 2 [f − kr · [t − t0 (τ )]]|t˜ = 0 dt d d φ(t)|t˜ = − π · 2 fτ + [f0 + kr · [t − t0 (τ )]] · t0 (τ ) dτ dτ |t˜ !
= 0. From the first equation in (34), we obtain
f f ! τ)= ⇒ t˜ = t0 (˜ τ )+ . [f −kr · [t−t0 (τ )]]|t˜ = 0 ⇒ t˜−t0 (˜ kr kr (35) Substituting this result into the second part of (34), we get fτ + (f0 + f ) ·
fτ +(f0 +f ) ·
(36)
τ −τ0 ) ! τ) · c 2 v 2 · (˜ fτ · R(˜ · = 0 ⇒ τ˜ −τ0 = − 2 c R(˜ τ) 2v · (f0 +f ) (37)
which can be solved straightforward, yielding the following two results for the point of stationary phase and the slant range at this point: |f0 + f | R(˜ τ ) = R0 · (f0 + f )2 − fτ2 ·
c2 4v 2
R0 · sign(f0 + f ) · 2vc 2 τ˜ − τ0 = − fτ · . c2 (f0 + f )2 − fτ2 · 4v 2
(38)
Substituting these two results into (33), we readily obtain
φ(t˜) = − π · 2f · t˜ + 2fτ · τ˜ + 2f0 · t0 (˜ τ) −kr · [t − t0 (˜ τ )]2 R0 f2 − 2πfτ τ0 − 4π = −π· kr c 2 c · (f0 + f )2 − fτ2 · 2 4v
−∞
where
Equation (32) is in the form that we can readily apply (16). For convenience and shortness of the example, we will restrict our interest to the phasor and neglect the window function. The phasor in (32) is given by φ(t) = −π· 2f ·t+2fτ ·τ +2f0 ·t0 (τ )−kr ·[t−t0 (τ )]2 . (33)
d ! t0 (τ )|˜τ = 0. dτ
Differentiating (31) with respect to and substituting the result back into (36), we get
(32)
t = [t, τ ]T f = [f, fτ ]T w(t) = w(t, τ ) t − t0 (τ ) τ − τc = rect · rect TR Taz
(34)
(39)
which is the well-known result for the 2-D point target reference phasor. Summarizing, we obtain, applying (16) ∞ w(t) · ejΦ(t) dt
ST (f, f )τ = −∞
∼ = (j2π) ·
2
i=1
−1
dii 2 (t˜) ·w(t˜)
A(f,fτ )
·e
−jπ·
f2 kr
+2fτ τ0 +4·
R0 c
·
(f0 +f )2 −fτ2 ·
c2 4v 2
. (40)
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
2715
The first amplitude term (usually ignored in the literature) has not been evaluated in detail here. The window term w(˜t) with its stationary point dependence on the range and Doppler frequency provides the spectral limitation of the result, particularly that the analysis of rect((˜ τ − τc )/Taz ) gives the commonly known results for Doppler centroid frequency and Doppler bandwidth of the point target reference spectrum. C. Arbitrary Chirp Spectrum In our last application example, we look at a quickly oscillating signal, where the phase is again directly given by a quadratic form, and the rectangular envelope signal with its argument being given by some other quadratic form. Although appearing rather artificial, this example in 1-D describes the monostatic case in azimuth direction. The phase history of any point target is centered on the point of closest approach. If we approximate this phase history by a quadratic function (chirp approximation), we may think of the squinted antenna beam as cutting out a “windowed” (and due to the squint, shifted) part of the infinitely extended quickly oscillating azimuth chirp, where the slope of the window function (envelope) is given by the two-way antenna beam. It is readily observed from (16) and has been discussed previously that this envelope is reproduced in the frequency domain. Shifting the envelope in the time domain, cutting out a different “portion” of a chirp reproduces a spectrum with the envelope shifted in the frequency domain (introducing the azimuth Doppler centroid in SAR). For 2-D problems, we may think of more “fancy” signals where the center of the 2-D envelope is not aligned with the phase minimum (maximum) of the oscillating part, such as in azimuth time-variant bistatic SAR processing. Here, the transmitter and receiver move at different velocities on nonparallel tracks, each one contributing individually to the overall phase history. An example for this case is described in [25]. Further applications are in the design and analysis of aspheric lenses by means of Fourier optics. In such applications, one would be interested in the transfer function of some noncentric and/or noncircular cutout of some 2-D lens with arbitrary higher order phase response. In applying the method successfully to such signals, we demonstrate its further potential and compare the results with the numerical solution obtained by fast Fourier transform (FFT) techniques. Again, for simplicity, we look at a 2-D signal
s(t) = rect (t − t0 )T ·A · (t−t0 ) ·exp jπ · tT · Q · t
(41)
w(t−t0 )
where
2558.8354 −1687.8297 −1687.8297 3154.0552 114.25781 0 0.8 · tdr A= t0 = . 0 10.1562 0.5 · tdaz
Q=
The sampling time in range time direction is TR = 0.000256410, and in azimuth direction, the sampling time is Taz = 0.000384615. Q models an elliptical chirp with its phase
Fig. 1. Elliptical infinite chirp.
Fig. 2. Envelope w(t − t0 ).
minimum at zero, which is counterclockwise rotated by 40◦ . The matrix A in the second window models an elliptical support. The whole envelope window is shifted to the point t0 , where tdr = 0.262564 and tdaz = 0.393846 are the time extents of the support window. All the signals [elliptical chirp, envelope, and outcoming signal s(t)] are shown in Figs. 1–3. Conceptually, the shifted envelope sifts out a certain frequency support from the infinite chirp in Fig. 1, which leads to a spectrum with bandpass properties. The point of stationary phase is the same as that given in (25): ˜t = Q−1 · f . Directly applying (16), we obtain, using the simplifications already made in (28) S(f ) = j
1 det(Q) C
·w(Q−1 ·f − t0 ) · exp{−jπ·f T · Q−1 ·f }
2716
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
Fig. 3. Outcoming signal s(t). Fig. 5.
Modulus of FFT spectrum.
Fig. 6.
Modulus of msp spectrum.
Fig. 4. Envelope of msp spectrum.
where ⎞T
⎡⎛
⎤
⎟ ⎢⎜ ⎥ w(Q−1 ·f −t0 ) =rect⎣⎝f −Q·t0⎠ · Q−1TAQ−1 ·(f −Q · t0 )⎦ f0
˜ A
(42) with
0.00060400077 0.00032321896 −0.0040622155 ˜ A= −0.0020697351 205.11408 f0 = 266.57544 C = 0.00043760684 · j.
Q−1 =
0.00032321896 0.00049001633 −0.0020697351 −0.00094978919
In order to verify the result, we first calculate the FFT spectrum of s(t). First, we plot the envelope given in (42) in Fig. 4. The modulus of that spectrum is shown in Fig. 5. Finally, Fig. 6 shows the modulus of the msp spectrum, calculated by (42). Referring to the visual comparison, the results are in excellent agreement. In a next step, we calculate the interferogram of the FFT and msp spectra, shown in Fig. 7. No fringe patterns are visible, indicating no phase errors. The inverse transformation of the spectral interferogram gives the matched filter result of the given chirp, compressed in the frequency domain with the msp spectrum. This result is shown in Fig. 8. Comparing this result with the compression result of the chirp compressed with its FFT spectrum, which is shown in Fig. 9, we see a further excellent agreement.
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
Fig. 7.
FFT spectrum multiplied with conjugate msp spectrum.
2717
Fig. 9. Matched filter result of chirp compressed with FFT spectrum.
A PPENDIX A Quadratic Term Approximation We want to approximate the phase history by a second-order quadratic form around the stationary point, where the first-order derivative is zero 1 Φ(x) ∼ x) + (x − x ˜)T P (˜ x) · (x − x ˜). = Φ(˜ 2
(43)
P (˜ x) is a symmetric matrix, explained hereinafter, and x ˜ is the point of stationary phase being characterized by dΦ(x) dΦ(x) dΦ(x) dΦ(x) dΦ(x) = , , ,. . . , = 0T . (44) dx |˜x dx1 dx2 dx3 dxn |˜x Now, letting x − x ˜ = u, we obtain 1 x) + uT · P (˜ x) · u Φ(u) ∼ = Φ(˜ 2 Fig. 8.
Matched filter result of chirp compressed with msp spectrum.
IV. C ONCLUSION In this paper, we have presented the derivation and conceptual framework of the n-variate (multidimensional) msp. Rather than aiming at fundamentally new results, the focus is on the derivation, comprehensive formulation, and treatment of the problem, yielding a closed-form solution for the general case. The solution has been successfully tested in three examples of varying complexity, starting with a simple case, where an analytic solution is available, continuing with a monostatic SAR example, and ending up with some numerical problem, where the asymptotic solution was successfully verified with FFT results, which were numerically obtained. However, more than demonstrating the validity of the results, it was again the demonstration of the simple applicability of the derived framework which motivated the examples.
(45)
where Φ(u) is a scalar, depending on the vector u. For the second-order derivative, we have T d d d 1 P (˜ x) + P (˜ Φ(u) = x )T · u du du du 2 =
d P0 (˜ x) · u = P0 (˜ x) du
(46)
x) = (1/2)[p(˜ x) + P (˜ x)T ]. where P0 (˜ x), we obtain from (45) the Using the symmetric matrix P0 (˜ second-order phase history approximation as 1 x) + uT ·0 (˜ x) · u. Φ(u) ∼ = Φ(˜ 2
(47)
Resubstituting u, we may write 1 x) + [x − x ˜]T · P0 (˜ x) · [x − x ˜] Φ(x) ∼ = Φ(˜ 2
(48)
2718
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
where the matrix P0 and point of stationary phase are given by ' T '' d dΦ(x) '' d ' P0 (˜ x) = Φ(x) = 0T . (49) ' ' dx dx dx 'x=˜x x=˜ x
Conceptually, our assumption of a nonvanishing secondorder derivative in the point of stationary phase (formulated in the introduction) is expressed by requiring nonvanishing eigenvalues of P0 (˜ x), and as the next section will show, the larger x) are, the more exact the approximate these eigenvalues of P0 (˜ solution will be.
Both sides of (54) become equal if 1 x) = −P −1 => P = − P0 (˜ x)−1 = jP0 (˜ x)−1 . jP0 (˜ j
(55)
We note that the eigenvalues of P0 (˜ x) are equal in magnitude to the inverse of the eigenvalues of P . Hence, requiring that the eigenvalues of P be sufficiently small would force the x) to be sufficiently large. Normally, P0 (x) eigenvalues of P0 (˜ will not necessarily be of diagonal form, but by introducing a principal component transformation, we can assure the desired diagonal form. Hence, the inverse becomes very simple: x)−1 = diagi=1,n [d−1 P0 (˜ ii ]. For the determinant, we then obtain 1 1 n − 12
2 2 n n −1 det P = det jP0 (˜ x)−1 = (j) 2 · det P0 (˜ x) = (j) 2 · dii 2 .
A PPENDIX B Dirac Approximations by Gaussians
i=1
We further need a description for multivariate Dirac impulses in terms of multivariate Gaussian signals. This interpretation is very much inspired by [20] and extends the idea of Papoulis [21], who used this approach in the 1-D case. We start with a multivariate Gaussian with unit area −1 1 2 −1 n f (x) = (2π) det[P ] ·exp [x−x0 ]T ·P −1 [x−x0 ] (50) 2 where x0 = E{x} could be considered as some mean value if f (x) were interpreted as a probability distribution density. By Fourier transformation arguments, it can be shown that f (x) approximates a Dirac impulse in the sense that −1 1
(56) with dii denoting the eigenvalues of P0 (and D), which are nonzero by assumption. n gives the number of vector elements. Summarizing, we have the equivalence − 12 n n 1 T ∼ x) · x = (2π) 2 · (j) 2 · det P0 (˜ x) · δ(x) . exp j x · P0 (˜ 2 (57) ACKNOWLEDGMENT The author would like to thank the anonymous reviewers for their valuable comments and criticism.
2
lim f (x) = lim (2π)n det[P ]
P →0
R EFERENCES
P →0
1 T −1 · exp − [x − x0 ] · P [x − x0 ] 2 = δ(x − x0 )
(51)
where P → 0 means that all eigenvalues of P go to zero. For sufficiently small eigenvalues, we may then approximately write −1 1 2 n −1 T −1 2 [x − x0 ] · P [x − x0 ] (2π) det P · exp 2 ∼ = δ(x − x0 ). Conversely, we would have then 1 2 n 1 T exp − x · P −1 · x ∼ = (2π) 2 det P · δ(x). 2
(52)
(53)
Now, we intend to substitute the real exponential by a complex quadratic form, to express the complex integrand in (2) or (12) in terms of a Dirac impulse. Hence, we require 1 T 1 T −1 x) · x = exp − x · P · x . (54) exp j x · P0 (˜ 2 2
[1] J. Fortuny and A. J. Sieber, “Fast algorithm for a near-field synthetic aperture radar processor,” IEEE Trans. Antennas Propag., vol. 42, no. 10, pp. 1458–1460, Oct. 1994. [2] G. Franceschetti and R. Lanari, Synthetic Aperture Radar Processing. Boca Raton, FL: CRC Press, 1999. [3] O. Loffeld and A. Hein, “SAR processing by inverse scaled Fourier transformation,” in Proc. EUSAR, Königswinter, Germany, 1996, pp. 143–146. [4] J. Mittermayer, A. Moreira, and O. Loffeld, “Spotlight SAR data processing using the frequency scaling algorithm,” IEEE Trans. Geosci. Remote Sens., vol. 37, no. 5, pp. 2198–2214, Sep. 1999. [5] K. Natroshvili, O. Loffeld, H. Nies, A. M. Ortiz, and S. Knedlik, “Focusing of general bistatic SAR configuration data with 2-D inverse scaled FFT,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 10, pp. 2718–2727, Oct. 2006. [6] O. Loffeld, H. Nies, V. Peters, and S. Knedlik, “Models and useful relations for bistatic SAR processing,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 10, pp. 2031–2038, Oct. 2004. [7] S. A. Greenhalgh and L. Marescot, “Modeling and migration of 2-D georadar data: A stationary phase approach,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 9, pp. 2421–2429, Sep. 2006. [8] R. Mager and N. Bleistein, “An examination of the limited aperture problem of physical optics inverse scattering,” IEEE Trans. Antennas Propag., vol. AP-26, no. 5, pp. 695–699, Sep. 1978. [9] G. L. James, Geometrical Theory of Diffraction for Electromagnetic Waves, 3 Sub. London, U.K.: IET, 1986. [10] N. G. Van Kampen, “An asymptotic treatment of diffraction problems,” Physica, vol. 14, no. 9, pp. 575–589, Jan. 1949. [11] D. S. Jones, “Asymptotic expansion of multiple integrals and the method of stationary phase,” J. Math. Phys., vol. 37, pp. 1–28, 1956. [12] N. G. Van Kampen, “The method of stationary phase and the method of Fresnel zones,” Physica, vol. 24, no. 1–5, pp. 437–444, 1958. [13] B. York, “Stationary phase method,” in ECE 201C Antenna Theory, University of California in Santa Barbara, Electrical and Computer Engineering Department, 2010. [Online]. Available: http://my.ece.ucsb. edu/bobsclass/201C/Handouts/StationaryPhase.pdf
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
[14] J. Lega, Method of stationary phase, vol. Math 583 B, 1999. [Online]. Available: http://math.arizona.edu/~lega/583/Spring99/lectnotes/AE4. html [15] N. Bleistein and R. A. Handelsman, “Multidimensional stationary phase an alternative derivation,” SIAM J. Math. Anal., vol. 6, no. 3, pp. 480–487, 1975. [16] R. Wong and J. McClure, “On a method of asymptotic evaluation of multiple integrals,” Math. Comput., vol. 37, no. 156, pp. 509–521, Oct. 1981. [17] G. Meinardus, “Remark on a lemma by R. Wong and J. P. McClure,” Math. Comput., vol. 45, pp. 197–198, Jul. 1985. [18] N. Bleistein, “Mathematical methods for wave phenomena,” Amsterdam, The Netherlands, Elsevier, 1984. [Online]. Available: http://knovel.com/ web/portal/browse/display?_EXT_KNOVEL_DISPLAY_bookid=2844& VerticalID=0 [19] R. Wong, Asymptotic Approximations of Integrals. Philadelphia, PA: SIAM, 2001. [20] M. J. Lighthill, An Introduction to Fourier Analysis and Generalised Functions. Cambridge, U.K.: Cambridge Univ. Press, 2003. [21] A. Papoulis, Systems and Transforms With Applications in Optics. New York: McGraw-Hill, 1968. [22] A. Papoulis, Signal Analysis. New York: McGraw-Hill, 1987. [23] G. Franceschetti and G. Schirinzi, “A SAR processor based on twodimensional FFT codes,” IEEE Trans. Aerosp. Electron. Syst., vol. 26, no. 2, pp. 356–366, Mar. 1990. [24] R. C. Thor, “A large time-bandwidth product pulse-compression technique,”IRE Trans. Mil. Electron., vol. MIL-6, no. 2, pp. 169–173, Apr. 1962. [25] R. Y. Wang, O. Loffeld, H. Nies, J. H. G. Ender, I. Walterscheid, and T. Espeter, “Processing the azimuth-variant bistatic SAR data by using monostatic imaging algorithm based on 2D principle of stationary phase,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 10, pp. 3504–3520, Oct. 2011.
2719
Otmar Loffeld (M’05–SM’06) received the Diploma degree in electrical engineering from the Technical University of Aachen, Aachen, Germany, in 1982, and the Eng. Dr. degree and the Habilitation in the field of digital signal processing and estimation theory in 1986 and 1989, respectively. Since 1991, he has been a Professor for digital signal processing and estimation theory with the University of Siegen. In 1999, he was a Principal Investigator (PI) on baseline estimation for the X-band part of the Shuttle Radar Topography Mission (SRTM). He lectures on general communication theory, digital signal processing, stochastic models, estimation theory, and synthetic aperture radar. Since 2005, he has been the Chairman of the Center for Sensorsystems (ZESS), University of Siegen. He is a PI for interferometric techniques in the German TerraSAR-X mission and a PI for a bistatic spaceborne airborne experiment, where TerraSAR-X serves as the bistatic illuminator, while the airborne “Phased Array Multifunctional Imaging Radar” system of the Fraunhofer Institute for High Frequency Physics and Radar Techniques (FHR) is used as a bistatic receiver. He is currently a PI in the “TerraSAR-X add-on for Digital Elevation Measurement mission” and in bistatic “Hitchhiker” experiments with stationary receiver and TerraSAR-X illuminator. He founded the International Postgraduate Program “Multi Sensorics” at the University of Siegen in 2002, and in 2008, he established the “North Rhine-Westphalia Research School on Multi Modal Sensor Systems for Environmental Exploration and Safety.” Prof. Loffeld is a member of the Information Technology Society (ITG/VDE) and a senior member of the IEEE Geoscience and Remote Sensing Society.
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
An Alternative Derivation of the Method of Stationary Phase for Multivariate Signals Otmar Loffeld, Senior Member, IEEE
Abstract—The method of stationary phase (msp) is a commonly established technique to asymptotically determine integrals (such as Fourier integrals) with strongly oscillating integrands. Basically, the method states that an integral over a strongly oscillating integrand only delivers contributions at those (“stationary”) points or regions, where the integrand does not change, i.e., where the phase derivative of its oscillation vanishes. In the engineering literature, various approaches for the derivation exist for univariate (1-D) signals, which, however, find their limitations for multivariate (multidimensional) signals, which cannot be factorized. On the other hand, from the field of applied mathematics, there are derivations, formulating and solving the problem in a rigorous and sophisticated way. This paper tries to close what the author interprets as a gap between these two worlds, aiming to substitute intuitive insights by some basic Dirac impulse modeling concepts, which are familiar to signal processing engineers. We first start with the general multivariate (multidimensional) case, introducing vectorial time and frequency representations. Here, we replace 1-D Taylor-series-based quadratic phase signals by multivariate signals with phase histories, being given by more general quadratic forms. The exponential, containing the quadratic form, is then related to an n-variate Dirac impulse, which provides the sifting property needed to establish the basic msp property that the integration over the entire integration range yields only one value of the integrand, evaluated in the point of stationary phase. This sifting capacity in connection with (time limiting) window functions in the integrand migrates time limitations of the signals to frequency limitations of the spectra—a multiplicative time window is transformed in a multiplicative frequency window of the same form, rather than into a convolution of the Fourier transforms. In this respect, this paper provides a closed-form comprehensive treatment of the subject. Index Terms—Chirp signals, monostatic point target reference spectrum, multidimensional method of stationary phase (msp), n-variate msp, quadratic phase history.
I. I NTRODUCTION
T
HE calculation of Fourier spectra S(f ) of multivariate signals s(t), where t = [t1 , t2 , . . . , tn ]T is some vector, denoting the time, and/or space dependence of a
Manuscript received August 17, 2010; revised December 20, 2010, March 7, 2011, and July 29, 2011; accepted October 15, 2011. Date of publication December 23, 2011; date of current version June 20, 2012. This work was supported in part by the German Science Foundation (DFG) under Grant DFG Lo 455/7-1, 7-2. The author is with the Center for Sensor Systems (ZESS), University of Siegen, 57068 Siegen, Germany (e-mail: [email protected]). Digital Object Identifier 10.1109/TGRS.2011.2174245
scalar-valued signal s(·) involves the solution of integrals of the form ∞ S(f ) = −∞ ∞
=
s(t) · e−j2πf
T
·t
dt
|s(t)| · ejφs (t) · e−j2πf
T
·t
dt
−∞ ∞
a(t) · ejφ(t,f ) · dt.
=
(1)
−∞
In those cases, where the solution of such integrals is not analytically possible, asymptotic solutions, like the method of stationary phase (msp), come into the focus of interest. This method, which states that the stationary points give the main contributions to such integrals, has been and is extensively used in, e.g., monostatic and bistatic synthetic aperture radar (SAR) systems [1]–[6], geophysics, antenna theory, and wave theory [7], [8] or seismics when multidimensional formulations are included. a(t) is some envelope or window function, limiting the time extent of the integrand in a way that a(t) sufficiently fast decays for large values of t. This function may even be discontinuous at the ends of the interval of integration but is supposed to contain neither singularities nor heavy oscillations. The phasor function φ(t, f ) is supposed to possess time derivatives up to second order and a nonvanishing and finite secondorder derivative in the interval of integration, particularly at the point of stationary phase, e.g., where the first-order derivative with respect to time vanishes. This phasor introduces an oscillating behavior (depending on f ) of the integrand, and the msp basically states that the integral only delivers a contribution at those positions t, where the phase is stationary, e.g., where its first-order derivative with respect to time t vanishes. In ray or diffraction theory, such points of stationary phase (stationary points) are often called critical points of the first kind [9]. They correspond to rays leaving the reflector normal to the tangent plane. In addition, there are critical points of the second and third kinds, due to diffraction at an edge and a corner, respectively, where cancellations in the exponential phase factor may become ineffective in the neighborhood of the limits of integration. They are considered in [10] and [11] in connection with the evaluation of infinite integrals, arising in optical diffraction problems. Readers are also referred to [12] for a further discussion. While the msp is readily formulated for univariate signals s(t), the extension to the general case,
0196-2892/$26.00 © 2011 IEEE
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
particularly where the signal s(t) cannot be factorized in a product of 1-D signals s(t) = ni=1 Si (ti ), needs some special consideration. In the multivariate case, this point of stationary phase becomes multidimensional and, due to the dependence of the phase on the integration variable t, may lead to cases where the multidimensional integral cannot be decoupled into individual integrals with a decoupled determination of the individual points of stationary phase for each integral. On the one hand, simple (and possibly too simple) extensions from the 1-D (univariate) case to the 2-D case have been published in [13] and [14]; on the other hand, this method has received ample treatment in the mathematical literature, and several rigorous expositions are available, e.g., [15], [16], and references therein. The first citation treats the subject using integration by parts, and the second one treats the subject using generalized Dirac functions and their derivatives. A multiple integral is reduced to a single line integral, using Dirac line functions, and then solved, where a proof, contained in the latter reference, has been further simplified in [17]. The advantage of such mathematic rigor is a precise formulation and understanding of preconditions and constraints, under which the approximation holds, as well as error bounds of the approximation. On the other hand, these works do not seem common reading, neither in signal processing nor in geoscience and remote sensing. The multidimensional stationary phase method has also appeared in book form, for example, [18] and [19], treating the subject in a comprehensive and thorough (excellent) way, but here, the price is paid in form of the length. This paper seeks a compromise; it substitutes an oversimplistic [13], [14] approach by an engineering mathematical approach, without, however, claiming the mathematic rigor in [15] and [16]. We start with the general multivariate (multidimensional) case, introducing vectorial time and frequency representations and replacing 1-D Taylor-series-based quadratic phase signals by multivariate signals with phase histories being given by more general quadratic forms. Summarizing the problem, we then look at integrals of the form
spectrum of an infinitely extended chirp and compare it with the analytical solution, which can be obtained. The second example considers the determination of the monostatic SAR point target reference spectrum, for which the results are commonly known, and the third example looks at an arbitrary chirp with some general quadratic form phase history and some other elliptical or hyperbolic envelope, to demonstrate the potential of the framework formulation. II. n-VARIATE MSP Recalling (2), we want to obtain the solution of integrals of the form ∞ A(x) · ejΦ(x) dx
I= −∞
=
···
A(x1 , . . . , xn ) · ejΦ(x1 ,...,xn ) dx1 , . . . , dxn (3)
n
where ejΦ(x) constitutes a scalar and quickly oscillating term with a large second-order derivative matrix (with all eigenvalues being sufficiently large). Furthermore, we assume that the assumptions concerning the window function A(x) and phase argument Φ(x), as discussed in the introduction, are met. Using the results in Appendix A, we suppose in the first step that we may approximate with sufficient exactness 1 x) + (x − x ˜)T · P0 (˜ x) · (x − x ˜) Φ(x) = Φ(˜ 2
A(x) · ejΦ(x) dx.
(2)
−∞
The dependence of the integral and the phasor φ(t, f ) on the parameter f has been omitted for simplicity. In the next section, the phase function Φ(x) will be approximated by a secondorder Taylor series approximation, expressed by a quadratic form. Such exponentials, containing quadratic forms, constitute multivariate Gaussians (cf. [20] for the scalar case). In the Appendix, we further show how the multivariate Dirac impulses can be interpreted as limiting cases of Gaussians, bringing a Dirac function into our integral. Conceptually, these Dirac functions, used in the integrand and being only nonzero, where their argument vanishes, give rise to line integration concepts, as discussed in [12]. The idea of expressing the implicit msp sifting property by Dirac impulses has already been used in [21] but here with the restriction to scalar signals. Having derived the formula framework, we apply this technique in the final chapter to three different applications. First, we calculate the n-variate
(4)
x ˜ is the point of stationary phase defined in Appendix A, given x) = (d/dx)[[(d/dx)Φ(x)]T ]x˜ by x ˜ : (dΦ(x)/dx)|x˜ = 0T . P0 (˜ is the second-order derivative of the multivariate phase history in the point of stationary phase, defined in Appendix A. We further assume that all eigenvalues of the matrix P0 are sufficiently large. Then, we have from (4) I∼ = I˜ =
∞ I=
2711
∞
1
A(x) · eΦ(˜x) · ej 2 (x−˜x) −∞ ∞
1
A(x) · ej 2 (x−˜x)
= eΦ(˜x) −∞
T
T
·P0 (˜ x)·(x−˜ x)
·P0 (˜ x)·(x−˜ x)
dx
dx .
(5)
I1
Looking at the integral I1 and then substituting x − x ˜ = u, where du = dx, we may write ∞
1
A(x) · ej 2 (x−˜x)
I1 = −∞ ∞
1
T
A(u + x ˜ ) · ej 2 u
=
·P0 (˜ x)·(x−˜ x)
T
·P0 (˜ x)·u
dx
du.
(6)
−∞
We now introduce an invertible orthogonal transformation of the form v = T · u => u = S · v where S −1 = T and T T = T −1 .
(7)
2712
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
Then, we have for the integration volume element
Substituting (15) into (5) and further using (3), the multivariate msp becomes
dv = |det(T )| · du
(8)
∞
−1
T
and due to T = T , we have det(T ) = 1. Substituting this into (6), we find that ∞ A(S · v + x ˜) · e
I1 =
j 12 (S·v)T ·P0 (˜ x)·S·v
−∞
dv
i=1
1
T
·S T ·P0 (˜ x)·S·v
1
T
·D(˜ x)·v
dv
−∞
∞
A(S · v + x ˜ ) · ej 2 v
dv.
(9)
−∞
Now, we require the matrix in the exponential to take on a diagonal form S T · P0 (˜ x) · S = D(˜ x) = diag (dii (˜ x)) .
(10)
i=1,n
It can be shown that the solution of (10) is a column vector matrix S, where the column vectors si , i = 1, . . . , n, obey the following equation: x) · si = dii (˜ x) · si , P0 (˜
i = 1, . . . , n
(11)
x) and si denotes saying that dii denotes the eigenvalues of P0 (˜ the eigenvectors corresponding to dii . Furthermore, due to the orthogonality of T and S, we also have: det(D(˜ x)) = x)). For the integral (9), we have then det(P0 (˜ ∞ I1 = −∞
1 T A(S · v + x ˜) · exp j v · D(˜ x) · v dv. 2
(12)
In Appendix B, we show that, provided that the eigenvalues of the diagonal matrix D are sufficiently large, the complex exponential with the quadratic form, given in (12), can be approximated by a multivariate Dirac impulse. Using this equivalence as specified in (57) in Appendix B, we obtain due to the diagonal form of D ∞
n
−∞ n n ∼ = (2π) 2 · (j) 2 ·
1 1 2
1
n
A(˜ x + S · v) · (2π) 2 · (j) 2 ·
I1 =
(16)
1
A(S · v + x ˜ ) · ej 2 v
=
−1
dii 2 (˜ x)
=det− 2 (P0 (˜ x)).
∞
!
n
n ∼ x) · ejΦ(˜x) · (j2π) 2 · = A(˜
−∞
=
A(x) · ejΦ(x) dx
I=
det D(˜ x)
1 2
det D(˜ x)
· δ(v)dv
· A(˜ x)
(13)
x ˜ : (dΦ(x)/dx)|x˜ = 0T is the multidimensional point of stationary phase, defined in (44), which, in the case of a Fourier transformation integral, depends on the frequency. x) denotes the eigenvalues of P0 (˜ x), where P0 (˜ x) = dii (˜ (d/dx)[[(d/dx)Φ(x)]T ]x˜ ) is the second-order derivative of the multivariate phase history in the point of stationary phase defined in (34). Two major items are worth noting in (16): First, we note that the product term of eigenvalues is equivalent to the x). Second, we note that rather than having determinant of P0 (˜ to deal with a convolution of Fourier transforms, we have a multiplication of the window function A(˜ x) with the phase history term, evaluated in the point of stationary phase. This essentially says that a multiplication with a strongly oscillating phase term carries over into the multiplication of the same parts in the frequency domain, indicating that multiplications do not necessarily always transform in convolutions in the frequency domain, if strongly oscillating terms are involved. This feature is implicitly introduced by the sifting property of the msp integral, being expressed by the Dirac impulse formalism in our derivation. It is essentially the deeper understanding of this feature which constitutes the implicit worth of the formal derivation. In the next part, we will look at application examples, further demonstrating the worth of this closed-form solution. III. E XAMPLES AND A PPLICATIONS A. Spectrum of the Infinitely Extended Chirp The first example is meant to provide a simple check of reasonableness with some infinitely extended chirp signal s(t) of elliptical geometry, for which its Fourier transform can be analytically determined. In order to ensure bounded total variation and, thus, existence of the involved Fourier integrals, we first assume s(t) = exp{jπ·tT ·Q·t} ⇒ S(f ) = F exp{jπ·tT ·Q·t} (17) where Q = Q1 + jQ2
where 1 2
D(˜ x) = diag (dii (˜ x)) → det (D(˜ x)) = i=1,n
n
1
dii (˜ x) 2 .
(14)
i=1
Hence, we obtain for the integral (13) n I1 ∼ = (j2π) 2 ·
n
i=1
−1
dii 2 (˜ x) · A(˜ x).
(15)
with Q1 being real, symmetric, and positive definite, and Q2 being real, symmetric, and nonnegative definite, with the following symmetries: QH = QT1 − jQT2 = Q1 − jQ2 = Q∗ . Conceptually, the matrix Q2 in the exponential ensures an exponential decay of the chirp envelope. Later on, we will obtain the spectrum of the infinite chirp as the limiting case, letting Q2 (more precisely its eigenvalues) go to zero.
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
Straightforward Solution: Rather than trying to evaluate the Fourier integral, we derive a differential equation describing the signal, which will then be transformed into the frequency domain, showing that the frequency spectrum basically obeys the same differential equation and thus will be again a chirp d s(t) = jπ·exp{jπ·tT ·Q·t}·tT ·Q = j2π·s(t)·tT ·Q. dt
Then, the corresponding phase at the stationary point is obtained as T φ(t˜) = − [2πf T − π t˜ Q1 ] · t˜ −1 = − [2πf T − πf T Q−1 1 Q1 ] · Q1 ·f = − πf T Q−1 1 f.
(18)
Transforming (18) into the frequency domain, we obtain, using (d/dt)s(t) = F −1 {j2πf T · S(f )} (19) j2πf T · S(f ) = F F j2π · s(t) · tT · Q .
2713
(26)
The eigenvalues of the matrix Q1 are given by 2 q11 + q22 q11 · q22 − q12 · 1± 1−4· (27) λ1,2 = 2 [q11 + q22 ]2 where
Now, using the dual relationship F {−j2π · s(t)} = (d/df )S(f ) for the spectrum, we obtain from (19) j2πf T · S(f ) = −
d S(f ) · Q df
d S(f ) = − j2π · S(f )·f T · Q−1 df
(20)
which is exactly the same differential equation as (18), except for the minus sign and the inverse matrix Q−1 . Hence, the wanted spectrum must be given by S(f ) = C · exp{−jπ · f T · Q−1 · f }
(21)
where the constant C can be determined by ∞ S(0) = C =
∞
n
Now, substituting the results of (25), (26), and (27) into (16), we obtain the final result for the spectrum −1 S(f ) = j2π · (2π)2 λ1 · λ2 · exp −π · f T · Q−1T Q2 · Q−1 1 1 f (28) · exp −jπ · f T · Q−1 1 ·f . Now, letting Q2 go to zero, we obtain the final result of (23) −1 S(f ) = j2π · (2π)2 λ1 · λ2 · exp −jπ · f T · Q−1 1 ·f 1 · exp −jπ · f T · Q−1 =j 1 · f . (29) 2 q11 · q22 − q12 C
−∞
B. Monostatic SAR Point Target Reference Spectrum
− 12
= j 2 · det(Q).
(22)
For n = 2, and letting Q2 go to zero, we obtain the spectrum of the infinite chirp S(f ) = j · det(Q1 ) · exp −jπ · f T · Q−1 1 ·f . − 12
(23)
We now determine the spectrum using the msp, given in (16) ∧ and t = x; A(x) = exp[−πxT · Q2 · x]. We have
φ(t) = − π 2f T · t − tT Q1 · t d φ(t) = − 2π[f T − tT Q1 ] dt T d d P0 (t) = φ(t) = 2π · Q1 . dt dt
(24)
The point of stationary phase t˜ is given by d ˜ T φ(t) = 0T ⇒ −2π[f T − t˜ Q1 ] = 0T dt −1 ˜ ⇒ t˜ = f T · Q−1 1 ⇒ t = Q1 ·f . T
d22 = 2πλ2
exp{jπ · tT · Q · t}dt
s(t)dt = −∞
d11 = 2πλ1
(25)
Our next application is the determination of the 2-D point target reference spectrum of a monostatic SAR sensor. This is an example, where the outcoming result is well known and well reported in the literature. Conventionally, the monostatic point target spectrum is calculated by first transforming the point target response into the range frequency domain and then subsequently into the azimuth frequency (Doppler) domain, using the 1-D msp. Furthermore, implicitly present in the classical understanding, the derivation must be done exactly in this sequence of 1-D transformations, since the other way round (first, an azimuth Fourier transformation, followed by the range Fourier transformation) would not work, since no stationary point could be identified. Our example will show, however, that, by addressing the problem in a 2-D way, we implicitly identify a 2-D stationary point. In doing so, the sequence of transform operations is not an issue anymore. The whole problem reduces to the solution of a system of equations, where we have two equations for the determination of a 2-D point of stationary phase. Considering the prerequisites of applying the msp (requiring a second-order derivative matrix with sufficiently large eigenvalues), there are conditions formulated for the 1-D case, such as in [22, pp. 272] or in the analysis given in [23] which, however, still need to be reformulated for the multivariate case. In some publications, like [24], this is reflected by requiring signals with sufficiently large time-bandwidth product, which
2714
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
seems to be necessary but not completely sufficient condition. We omit a deeper analysis here for simplicity and space. For the given example, the sensor moves along the azimuth axis, denoted by azimuth time τ , the fast time (range time) coordinate is denoted by t, the sensor’s azimuth velocity is v, and the point target position is described by the azimuth time instant of closest approach τ0 and by the (orthogonal) slant range vector R0 at that time, where R0 = R(τ, R0 , τ0 )|τ =τ0 = R02 + ν 2 (τ − τ02 )τ =τ0 . After quadrature down conversion to zero offset frequency, the point target response in the range azimuth time domain is then given by t − t0 (τ ) sT (t, τ ) = σ(R0 , τ0 ) · rect TR · exp jπkr (t − t0 (τ ))2 τ − τc · rect · exp (−j2πf0 · t0 (τ )) . (30) Taz kr is the chirp rate of the transmitted signal, and TR and Taz are the time extensions, specifying the time apertures of the signal, which are described by rectangular windows. τc is the azimuth time, when the point target is seen in the center of the (azimuth) antenna footprint (due to a squinted geometry), f0 is the carrier frequency, σ(R0 , τ0 ) describes the radar brightness of the point target, and the two-way propagation delay is given by 2 2 t0 (τ ) = R(τ, R0 , τ0 ) = R02 + v 2 (τ − τ02 ). (31) c c Implicitly, the point target response would be 4-D, since the signal clearly depends on the location of the point target, but this dependence is neglected for convenience. From (30), the Fourier spectrum of the signal normalized to its radar brightness is obtained by ST (f, fτ ) = ST (f ) ∞ ∞ t − t0 (τ ) τ − τc rect · rect = TR Taz −∞ −∞
∞
· exp [jφ(t, τ )] dtdτ
w(t) · exp [jφ(t)] dt
=
kr is the range chirp rate of the transmitted signal. In order to identify the point t˜ of stationary phase (16), the two partial derivatives of (33) must vanish d ! φ(t)|t˜ = − π · 2 [f − kr · [t − t0 (τ )]]|t˜ = 0 dt d d φ(t)|t˜ = − π · 2 fτ + [f0 + kr · [t − t0 (τ )]] · t0 (τ ) dτ dτ |t˜ !
= 0. From the first equation in (34), we obtain
f f ! τ)= ⇒ t˜ = t0 (˜ τ )+ . [f −kr · [t−t0 (τ )]]|t˜ = 0 ⇒ t˜−t0 (˜ kr kr (35) Substituting this result into the second part of (34), we get fτ + (f0 + f ) ·
fτ +(f0 +f ) ·
(36)
τ −τ0 ) ! τ) · c 2 v 2 · (˜ fτ · R(˜ · = 0 ⇒ τ˜ −τ0 = − 2 c R(˜ τ) 2v · (f0 +f ) (37)
which can be solved straightforward, yielding the following two results for the point of stationary phase and the slant range at this point: |f0 + f | R(˜ τ ) = R0 · (f0 + f )2 − fτ2 ·
c2 4v 2
R0 · sign(f0 + f ) · 2vc 2 τ˜ − τ0 = − fτ · . c2 (f0 + f )2 − fτ2 · 4v 2
(38)
Substituting these two results into (33), we readily obtain
φ(t˜) = − π · 2f · t˜ + 2fτ · τ˜ + 2f0 · t0 (˜ τ) −kr · [t − t0 (˜ τ )]2 R0 f2 − 2πfτ τ0 − 4π = −π· kr c 2 c · (f0 + f )2 − fτ2 · 2 4v
−∞
where
Equation (32) is in the form that we can readily apply (16). For convenience and shortness of the example, we will restrict our interest to the phasor and neglect the window function. The phasor in (32) is given by φ(t) = −π· 2f ·t+2fτ ·τ +2f0 ·t0 (τ )−kr ·[t−t0 (τ )]2 . (33)
d ! t0 (τ )|˜τ = 0. dτ
Differentiating (31) with respect to and substituting the result back into (36), we get
(32)
t = [t, τ ]T f = [f, fτ ]T w(t) = w(t, τ ) t − t0 (τ ) τ − τc = rect · rect TR Taz
(34)
(39)
which is the well-known result for the 2-D point target reference phasor. Summarizing, we obtain, applying (16) ∞ w(t) · ejΦ(t) dt
ST (f, f )τ = −∞
∼ = (j2π) ·
2
i=1
−1
dii 2 (t˜) ·w(t˜)
A(f,fτ )
·e
−jπ·
f2 kr
+2fτ τ0 +4·
R0 c
·
(f0 +f )2 −fτ2 ·
c2 4v 2
. (40)
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
2715
The first amplitude term (usually ignored in the literature) has not been evaluated in detail here. The window term w(˜t) with its stationary point dependence on the range and Doppler frequency provides the spectral limitation of the result, particularly that the analysis of rect((˜ τ − τc )/Taz ) gives the commonly known results for Doppler centroid frequency and Doppler bandwidth of the point target reference spectrum. C. Arbitrary Chirp Spectrum In our last application example, we look at a quickly oscillating signal, where the phase is again directly given by a quadratic form, and the rectangular envelope signal with its argument being given by some other quadratic form. Although appearing rather artificial, this example in 1-D describes the monostatic case in azimuth direction. The phase history of any point target is centered on the point of closest approach. If we approximate this phase history by a quadratic function (chirp approximation), we may think of the squinted antenna beam as cutting out a “windowed” (and due to the squint, shifted) part of the infinitely extended quickly oscillating azimuth chirp, where the slope of the window function (envelope) is given by the two-way antenna beam. It is readily observed from (16) and has been discussed previously that this envelope is reproduced in the frequency domain. Shifting the envelope in the time domain, cutting out a different “portion” of a chirp reproduces a spectrum with the envelope shifted in the frequency domain (introducing the azimuth Doppler centroid in SAR). For 2-D problems, we may think of more “fancy” signals where the center of the 2-D envelope is not aligned with the phase minimum (maximum) of the oscillating part, such as in azimuth time-variant bistatic SAR processing. Here, the transmitter and receiver move at different velocities on nonparallel tracks, each one contributing individually to the overall phase history. An example for this case is described in [25]. Further applications are in the design and analysis of aspheric lenses by means of Fourier optics. In such applications, one would be interested in the transfer function of some noncentric and/or noncircular cutout of some 2-D lens with arbitrary higher order phase response. In applying the method successfully to such signals, we demonstrate its further potential and compare the results with the numerical solution obtained by fast Fourier transform (FFT) techniques. Again, for simplicity, we look at a 2-D signal
s(t) = rect (t − t0 )T ·A · (t−t0 ) ·exp jπ · tT · Q · t
(41)
w(t−t0 )
where
2558.8354 −1687.8297 −1687.8297 3154.0552 114.25781 0 0.8 · tdr A= t0 = . 0 10.1562 0.5 · tdaz
Q=
The sampling time in range time direction is TR = 0.000256410, and in azimuth direction, the sampling time is Taz = 0.000384615. Q models an elliptical chirp with its phase
Fig. 1. Elliptical infinite chirp.
Fig. 2. Envelope w(t − t0 ).
minimum at zero, which is counterclockwise rotated by 40◦ . The matrix A in the second window models an elliptical support. The whole envelope window is shifted to the point t0 , where tdr = 0.262564 and tdaz = 0.393846 are the time extents of the support window. All the signals [elliptical chirp, envelope, and outcoming signal s(t)] are shown in Figs. 1–3. Conceptually, the shifted envelope sifts out a certain frequency support from the infinite chirp in Fig. 1, which leads to a spectrum with bandpass properties. The point of stationary phase is the same as that given in (25): ˜t = Q−1 · f . Directly applying (16), we obtain, using the simplifications already made in (28) S(f ) = j
1 det(Q) C
·w(Q−1 ·f − t0 ) · exp{−jπ·f T · Q−1 ·f }
2716
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
Fig. 3. Outcoming signal s(t). Fig. 5.
Modulus of FFT spectrum.
Fig. 6.
Modulus of msp spectrum.
Fig. 4. Envelope of msp spectrum.
where ⎞T
⎡⎛
⎤
⎟ ⎢⎜ ⎥ w(Q−1 ·f −t0 ) =rect⎣⎝f −Q·t0⎠ · Q−1TAQ−1 ·(f −Q · t0 )⎦ f0
˜ A
(42) with
0.00060400077 0.00032321896 −0.0040622155 ˜ A= −0.0020697351 205.11408 f0 = 266.57544 C = 0.00043760684 · j.
Q−1 =
0.00032321896 0.00049001633 −0.0020697351 −0.00094978919
In order to verify the result, we first calculate the FFT spectrum of s(t). First, we plot the envelope given in (42) in Fig. 4. The modulus of that spectrum is shown in Fig. 5. Finally, Fig. 6 shows the modulus of the msp spectrum, calculated by (42). Referring to the visual comparison, the results are in excellent agreement. In a next step, we calculate the interferogram of the FFT and msp spectra, shown in Fig. 7. No fringe patterns are visible, indicating no phase errors. The inverse transformation of the spectral interferogram gives the matched filter result of the given chirp, compressed in the frequency domain with the msp spectrum. This result is shown in Fig. 8. Comparing this result with the compression result of the chirp compressed with its FFT spectrum, which is shown in Fig. 9, we see a further excellent agreement.
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
Fig. 7.
FFT spectrum multiplied with conjugate msp spectrum.
2717
Fig. 9. Matched filter result of chirp compressed with FFT spectrum.
A PPENDIX A Quadratic Term Approximation We want to approximate the phase history by a second-order quadratic form around the stationary point, where the first-order derivative is zero 1 Φ(x) ∼ x) + (x − x ˜)T P (˜ x) · (x − x ˜). = Φ(˜ 2
(43)
P (˜ x) is a symmetric matrix, explained hereinafter, and x ˜ is the point of stationary phase being characterized by dΦ(x) dΦ(x) dΦ(x) dΦ(x) dΦ(x) = , , ,. . . , = 0T . (44) dx |˜x dx1 dx2 dx3 dxn |˜x Now, letting x − x ˜ = u, we obtain 1 x) + uT · P (˜ x) · u Φ(u) ∼ = Φ(˜ 2 Fig. 8.
Matched filter result of chirp compressed with msp spectrum.
IV. C ONCLUSION In this paper, we have presented the derivation and conceptual framework of the n-variate (multidimensional) msp. Rather than aiming at fundamentally new results, the focus is on the derivation, comprehensive formulation, and treatment of the problem, yielding a closed-form solution for the general case. The solution has been successfully tested in three examples of varying complexity, starting with a simple case, where an analytic solution is available, continuing with a monostatic SAR example, and ending up with some numerical problem, where the asymptotic solution was successfully verified with FFT results, which were numerically obtained. However, more than demonstrating the validity of the results, it was again the demonstration of the simple applicability of the derived framework which motivated the examples.
(45)
where Φ(u) is a scalar, depending on the vector u. For the second-order derivative, we have T d d d 1 P (˜ x) + P (˜ Φ(u) = x )T · u du du du 2 =
d P0 (˜ x) · u = P0 (˜ x) du
(46)
x) = (1/2)[p(˜ x) + P (˜ x)T ]. where P0 (˜ x), we obtain from (45) the Using the symmetric matrix P0 (˜ second-order phase history approximation as 1 x) + uT ·0 (˜ x) · u. Φ(u) ∼ = Φ(˜ 2
(47)
Resubstituting u, we may write 1 x) + [x − x ˜]T · P0 (˜ x) · [x − x ˜] Φ(x) ∼ = Φ(˜ 2
(48)
2718
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
where the matrix P0 and point of stationary phase are given by ' T '' d dΦ(x) '' d ' P0 (˜ x) = Φ(x) = 0T . (49) ' ' dx dx dx 'x=˜x x=˜ x
Conceptually, our assumption of a nonvanishing secondorder derivative in the point of stationary phase (formulated in the introduction) is expressed by requiring nonvanishing eigenvalues of P0 (˜ x), and as the next section will show, the larger x) are, the more exact the approximate these eigenvalues of P0 (˜ solution will be.
Both sides of (54) become equal if 1 x) = −P −1 => P = − P0 (˜ x)−1 = jP0 (˜ x)−1 . jP0 (˜ j
(55)
We note that the eigenvalues of P0 (˜ x) are equal in magnitude to the inverse of the eigenvalues of P . Hence, requiring that the eigenvalues of P be sufficiently small would force the x) to be sufficiently large. Normally, P0 (x) eigenvalues of P0 (˜ will not necessarily be of diagonal form, but by introducing a principal component transformation, we can assure the desired diagonal form. Hence, the inverse becomes very simple: x)−1 = diagi=1,n [d−1 P0 (˜ ii ]. For the determinant, we then obtain 1 1 n − 12
2 2 n n −1 det P = det jP0 (˜ x)−1 = (j) 2 · det P0 (˜ x) = (j) 2 · dii 2 .
A PPENDIX B Dirac Approximations by Gaussians
i=1
We further need a description for multivariate Dirac impulses in terms of multivariate Gaussian signals. This interpretation is very much inspired by [20] and extends the idea of Papoulis [21], who used this approach in the 1-D case. We start with a multivariate Gaussian with unit area −1 1 2 −1 n f (x) = (2π) det[P ] ·exp [x−x0 ]T ·P −1 [x−x0 ] (50) 2 where x0 = E{x} could be considered as some mean value if f (x) were interpreted as a probability distribution density. By Fourier transformation arguments, it can be shown that f (x) approximates a Dirac impulse in the sense that −1 1
(56) with dii denoting the eigenvalues of P0 (and D), which are nonzero by assumption. n gives the number of vector elements. Summarizing, we have the equivalence − 12 n n 1 T ∼ x) · x = (2π) 2 · (j) 2 · det P0 (˜ x) · δ(x) . exp j x · P0 (˜ 2 (57) ACKNOWLEDGMENT The author would like to thank the anonymous reviewers for their valuable comments and criticism.
2
lim f (x) = lim (2π)n det[P ]
P →0
R EFERENCES
P →0
1 T −1 · exp − [x − x0 ] · P [x − x0 ] 2 = δ(x − x0 )
(51)
where P → 0 means that all eigenvalues of P go to zero. For sufficiently small eigenvalues, we may then approximately write −1 1 2 n −1 T −1 2 [x − x0 ] · P [x − x0 ] (2π) det P · exp 2 ∼ = δ(x − x0 ). Conversely, we would have then 1 2 n 1 T exp − x · P −1 · x ∼ = (2π) 2 det P · δ(x). 2
(52)
(53)
Now, we intend to substitute the real exponential by a complex quadratic form, to express the complex integrand in (2) or (12) in terms of a Dirac impulse. Hence, we require 1 T 1 T −1 x) · x = exp − x · P · x . (54) exp j x · P0 (˜ 2 2
[1] J. Fortuny and A. J. Sieber, “Fast algorithm for a near-field synthetic aperture radar processor,” IEEE Trans. Antennas Propag., vol. 42, no. 10, pp. 1458–1460, Oct. 1994. [2] G. Franceschetti and R. Lanari, Synthetic Aperture Radar Processing. Boca Raton, FL: CRC Press, 1999. [3] O. Loffeld and A. Hein, “SAR processing by inverse scaled Fourier transformation,” in Proc. EUSAR, Königswinter, Germany, 1996, pp. 143–146. [4] J. Mittermayer, A. Moreira, and O. Loffeld, “Spotlight SAR data processing using the frequency scaling algorithm,” IEEE Trans. Geosci. Remote Sens., vol. 37, no. 5, pp. 2198–2214, Sep. 1999. [5] K. Natroshvili, O. Loffeld, H. Nies, A. M. Ortiz, and S. Knedlik, “Focusing of general bistatic SAR configuration data with 2-D inverse scaled FFT,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 10, pp. 2718–2727, Oct. 2006. [6] O. Loffeld, H. Nies, V. Peters, and S. Knedlik, “Models and useful relations for bistatic SAR processing,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 10, pp. 2031–2038, Oct. 2004. [7] S. A. Greenhalgh and L. Marescot, “Modeling and migration of 2-D georadar data: A stationary phase approach,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 9, pp. 2421–2429, Sep. 2006. [8] R. Mager and N. Bleistein, “An examination of the limited aperture problem of physical optics inverse scattering,” IEEE Trans. Antennas Propag., vol. AP-26, no. 5, pp. 695–699, Sep. 1978. [9] G. L. James, Geometrical Theory of Diffraction for Electromagnetic Waves, 3 Sub. London, U.K.: IET, 1986. [10] N. G. Van Kampen, “An asymptotic treatment of diffraction problems,” Physica, vol. 14, no. 9, pp. 575–589, Jan. 1949. [11] D. S. Jones, “Asymptotic expansion of multiple integrals and the method of stationary phase,” J. Math. Phys., vol. 37, pp. 1–28, 1956. [12] N. G. Van Kampen, “The method of stationary phase and the method of Fresnel zones,” Physica, vol. 24, no. 1–5, pp. 437–444, 1958. [13] B. York, “Stationary phase method,” in ECE 201C Antenna Theory, University of California in Santa Barbara, Electrical and Computer Engineering Department, 2010. [Online]. Available: http://my.ece.ucsb. edu/bobsclass/201C/Handouts/StationaryPhase.pdf
LOFFELD: ALTERNATIVE DERIVATION OF THE METHOD OF STATIONARY PHASE
[14] J. Lega, Method of stationary phase, vol. Math 583 B, 1999. [Online]. Available: http://math.arizona.edu/~lega/583/Spring99/lectnotes/AE4. html [15] N. Bleistein and R. A. Handelsman, “Multidimensional stationary phase an alternative derivation,” SIAM J. Math. Anal., vol. 6, no. 3, pp. 480–487, 1975. [16] R. Wong and J. McClure, “On a method of asymptotic evaluation of multiple integrals,” Math. Comput., vol. 37, no. 156, pp. 509–521, Oct. 1981. [17] G. Meinardus, “Remark on a lemma by R. Wong and J. P. McClure,” Math. Comput., vol. 45, pp. 197–198, Jul. 1985. [18] N. Bleistein, “Mathematical methods for wave phenomena,” Amsterdam, The Netherlands, Elsevier, 1984. [Online]. Available: http://knovel.com/ web/portal/browse/display?_EXT_KNOVEL_DISPLAY_bookid=2844& VerticalID=0 [19] R. Wong, Asymptotic Approximations of Integrals. Philadelphia, PA: SIAM, 2001. [20] M. J. Lighthill, An Introduction to Fourier Analysis and Generalised Functions. Cambridge, U.K.: Cambridge Univ. Press, 2003. [21] A. Papoulis, Systems and Transforms With Applications in Optics. New York: McGraw-Hill, 1968. [22] A. Papoulis, Signal Analysis. New York: McGraw-Hill, 1987. [23] G. Franceschetti and G. Schirinzi, “A SAR processor based on twodimensional FFT codes,” IEEE Trans. Aerosp. Electron. Syst., vol. 26, no. 2, pp. 356–366, Mar. 1990. [24] R. C. Thor, “A large time-bandwidth product pulse-compression technique,”IRE Trans. Mil. Electron., vol. MIL-6, no. 2, pp. 169–173, Apr. 1962. [25] R. Y. Wang, O. Loffeld, H. Nies, J. H. G. Ender, I. Walterscheid, and T. Espeter, “Processing the azimuth-variant bistatic SAR data by using monostatic imaging algorithm based on 2D principle of stationary phase,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 10, pp. 3504–3520, Oct. 2011.
2719
Otmar Loffeld (M’05–SM’06) received the Diploma degree in electrical engineering from the Technical University of Aachen, Aachen, Germany, in 1982, and the Eng. Dr. degree and the Habilitation in the field of digital signal processing and estimation theory in 1986 and 1989, respectively. Since 1991, he has been a Professor for digital signal processing and estimation theory with the University of Siegen. In 1999, he was a Principal Investigator (PI) on baseline estimation for the X-band part of the Shuttle Radar Topography Mission (SRTM). He lectures on general communication theory, digital signal processing, stochastic models, estimation theory, and synthetic aperture radar. Since 2005, he has been the Chairman of the Center for Sensorsystems (ZESS), University of Siegen. He is a PI for interferometric techniques in the German TerraSAR-X mission and a PI for a bistatic spaceborne airborne experiment, where TerraSAR-X serves as the bistatic illuminator, while the airborne “Phased Array Multifunctional Imaging Radar” system of the Fraunhofer Institute for High Frequency Physics and Radar Techniques (FHR) is used as a bistatic receiver. He is currently a PI in the “TerraSAR-X add-on for Digital Elevation Measurement mission” and in bistatic “Hitchhiker” experiments with stationary receiver and TerraSAR-X illuminator. He founded the International Postgraduate Program “Multi Sensorics” at the University of Siegen in 2002, and in 2008, he established the “North Rhine-Westphalia Research School on Multi Modal Sensor Systems for Environmental Exploration and Safety.” Prof. Loffeld is a member of the Information Technology Society (ITG/VDE) and a senior member of the IEEE Geoscience and Remote Sensing Society.
