9

STOLT INTERPOLATION PROCESSING ON SAR IMAGES

This chapter discusses an alternative approach for processing radar images from the wavenumber (or spatial frequency) domain. Section 9.1 reviews the general background of wavenumber theory. Section 9.2 describes the direct interpolation based on unevenly spaced samples. Section 9.3 reviews the algorithm of Stolt interpolation. Computer simulation for six-target broadside SAR system based on the Stolt interpolation technique is covered in Section 9.4. Computer simulation for six-target squint SAR system is described in Section 9.5. Section 9.6 shows the results of a satellite image file processed by both range–Doppler and Stolt interpolation algorithms. Comparison between range–Doppler and Stolt interpolation algorithms is discussed in Section 9.7.

9.1 WAVENUMBER DOMAIN PROCESSING OF SAR DATA

Several algorithms are used to process SAR image data in wavenumber domains. Examples are direct interpolation from unevenly spaced data, Stolt interpolation (or mapping), time-domain correlation and backprojection, and range stacking.

Consider a group of N targets randomly distributed on the ground with each target located at (x, y) = (x_n, y_n), n = 1, 2,..., N. Assume that the radar is located at (x, y) = (0, u), and moves at speed V along the y axis. Let p(t) be the transmitting signal; then the total received signal s(t, u) can be expressed as

where s(t, u) is a function of both t and u, and p(t) is the transmitting signal. The Fourier transform of s(t,u) with respect to time variable t becomes

where P(ω) is the Fourier transform of p(t).

By taking the spatial Fourier transform of S(ω, u) with respect to spatial variable u, as discussed in Section 7.2.1.3 (of Chapter 7), one obtains

Here ω_D is the Doppler angular frequency with respect to the radar position u and k_u = ω_D/V is the spatial wavenumber in the Doppler frequency domain.

Letting k_x be defined as

and k_y = k_u, one can then rewrite Eq. (9.3) as follows:

Consider an ideal target function in the spatial domain defined as

Here δ(x – x_n, y – y_n) corresponds to an ideal impulse-like target located at (x, y) = (x_n, y_n). By applying the 2D spatial Fourier transform on f₀(x, y), one obtains

Here k_x and k_y are the wavenumber in the spatial frequency domain, and correspond to x and y in the spatial domain, respectively.

Since we are dealing with a moving spatial variable u, which represents a moving radar position, the corresponding spatial frequency will be considered as Doppler frequency in the following discussion.

Comparing Eqs. (9.5) and (9.7), one can see that if the signal S(k, k_y) can be mapped, or transformed into a spatial frequency k_x as shown in Eq. (9.4), then the target function can be obtained by taking the 2D inverse spatial Fourier transform on S(k_x, k_y). However, the measured signal S(k, k_y) is evenly spaced in the time–frequency (ω) and spatial wavenumber (k_y = ω_D/V) domains, but the transformation from k to k_x, as shown in Eq. (9.4), causes the data to become unevenly spaced in the k_x domain. The inverse 2D spatial transform on S(k_x, k_y) cannot work correctly with a 2D data array, where the data are evenly spaced in the k_y domain but unevenly spaced in the k_x domain. This issue must be resolved to correctly reconstruct the radar image.

Figure 9.1 illustrates the data relationship before and after the square root transformation defined in Eq. (9.4). The black dots in Fig. 9.1a represent the data S(k, k_y) distributed in the (k, k_y) or (ω, ω_D) domain. The white circles in Fig. 9.1b represent the same but transformed data S(k_x, k_y) distributed in the (k_x, k_y) domain. As can be seen, the evenly spaced black dots become unevenly spaced white circles.

In order to use the inverse 2D spatial transform to obtain the target function, the unevenly spaced data in the k_x domain must be interpolated to become evenly spaced data. Figure 9.2 displays the ideal relationship of the data distribution before and after interpolation. The white circles shown in Fig. 9.2 are the unevenly spaced data S(k_x, k_y) after transformation. The black dots are the evenly spaced data S′(k_x, k_y) after interpolation.

The conversion of unevenly spaced data to evenly spaced data in the (k_x, k_y) domain, based on the received data in the (k, k_y) domain, is addressed next.

FIGURE 9.1 Data distribution before (a) and after (b) transformation.

FIGURE 9.2 Data distribution before () and after (•) interpolation.

9.2 DIRECT INTERPOLATION FROM UNEVENLY SPACED SAMPLES

As discussed in Section 1.7.5 (of Chapter 1), the unevenly spaced data can be interpolated to obtain the evenly spaced data. Let S_c′(k_x, k_y) and S_c(k_x, k_y) be the evenly spaced and unevenly spaced 2D data in the (k_x, k_y) domain, respectively. Let S_c(k, k_y) be the linearly measured samples in the (k, k_y) domain, with the unevenly spaced data S_c(k_x, k_y) as the direct transformation, through Eq. (9.4), from S_c(k, k_y).

In order to obtain the 2D spatial domain (x, y) signal from the inverse 2D transform on signals in the (k_x, k_y) domain, one method is to interpolate the unevenly spaced samples of S_c(k_x, k_y) to obtain the evenly spaced data S′_c(k_x, k_y), as illustrated in Fig. 9.2. Let the received signal be an array of (M × N) data. With respect to the mth row of data array {S_c(k_xmn, k_ym)}, n = 1, 2,..., N, the direct interpolated S_c′(k_xm, k_ym) can be expressed as

where

Here J_m(k) is the Jacobian of transformation from k to k_x, and index n must satisfy |k_xm – k_xmn| ≤ N_sΔk_x, with N_s as half the number of sinc lobes used for interpolation. The parameters Δk and Δk_x are described below.

Let f_s be the sampling frequency, with 2X₀ as the swath of target area; then Δk and Δk_x can be expressed as follows:

Let the carrier frequency and the bandwidth of the transmitting signal be f_c and 2f₀; then

The corresponding k_x is therefore band-limited as

For broadside SAR, the minimum value of Doppler frequency is zero and 2(k_c – k₀) ω_Dmax/V; therefore, the range of k_x becomes

The number of samples in the new k_x domain is then

The interpolation function h(k_x) is similar to the one used in range cell migration correction (RCMC). In practical applications, only a finite number of resolution and sidelobes are adopted for interpolation filtering. In the RCMC case, the interpolation filter is chosen with the number of fractional sample shifts equal to 16, and the number of sidelobes equal to 7; therefore the variable n is chosen to be n = 1, 2,..., 8. The direct implementation of Eq. (9.8) is time-consuming. An alternative wavenumber domain processing technique, namely, the Stolt interpolation (or transformation), is discussed next.

9.3 STOLT INTERPOLATION PROCESSING OF SAR DATA

To simplify the discussion, only the targets with ideal reflection coefficients, that is, σ_n = 1 for n = 1, 2,..., N, are considered here. Let s₀(t, u) be the received signal from a reference target located at (x, y) = (X₀, 0). Then the 2D Fourier transform of s₀(t, u) can be expressed as follows:

The function S₀(ω, ω_D) serves as a 2D reference function for matched filtering on both the ω and ω_D domains. Here P*(ω) is considered as the range matched filter, while

is the azimuth matched filter. The 2D compressed function S_c(ω, ω_D) can be represented as

or

where k = ω/c, and k_y = k_u = ω_D/V.

The 2D compressed signal S_c(k, k_u) shown in Eq. (9.11) serves as the basis for wavenumber domain image processing. Since the azimuth reference function was chosen at (X₀, 0), only the target located at (X₀, 0) can be correctly focused. Other targets away from (X₀, 0) will be unfocused and require more accurate compression or correction.

Let f_b, f_c and 2f₀ be the baseband frequency, carrier frequency, and bandwidth of transmitting signals, which satisfy the relationship ω_c – ω₀ ≤ ω = ω_c + ω_b ≤ ω_c + ω₀, or, in terms of wavenumber representation k_c – k₀ ≤ k = k_c + k_b ≤ k_c +k₀. The square root term in Eq. (9.11) can then be approximated as

where

Since k_c k_u and k_c > k_b, Eq. (9.12) can be represented as follows:

Equation (9.11) therefore becomes

Letting k_x = k_c + k_xb, where k_xb is the baseband of k_x, and from Eq. (9.13), one obtains

Therefore

The function can be approximated as

and

Equation (9.14b) then becomes

This equation illustrates the transformation from k to k_x in the baseband situation. The Stolt interpolation simplifies the nonlinear transformation of k_x = k_c + k_xb = into the form shown in Eq. (9.14c). The baseband component k_xb consists of three items: (1) a constant k_c that is independent of k or k_y; (2) a shift that varies for different k_y in the k domain, namely, and (3) a stretch of k_b in the k domain with stretch factor In other words, to perform the transformation of received baseband data from (k_b, k_y) domain to (k_xb, k_y) domain, the received baseband data must be adjusted as shown in Eq. (9.14c).

By assuming that |P(k)| = 1, and substituting Eq. (9.14c) in Eq. (9.11), one obtains

where and

The term

in Eq. (9.15) represents a stretched and new spatial frequency in spatial domain signal s_c(x, ω_D). The amount of stretched frequency is constant for a given value of k_y. It can be implemented by varying the sampling frequency f_s. This item can be neglected if only the magnitude of the radar image reconstruction is concerned.

Given s(x) and S(k) as the Fourier transform pair, the well-known modulation and shift property shows that

S(k – k₀) s (x)exp(jk₀x).

The term k₀ in Eq. (9.15) states that in the wavenumber (k_b) domain, S_c(k_b, k_y) has a k_y-dependent shift of This implies that the corresponding spatial domain function s_c(x, k_u), which is in the range–Doppler domain, must be multiplied by a factor

to correct or compensate the constant frequency shift:

Since k_c is a constant wavenumber corresponding to the carrier frequency of the radar-transmitting signal, it will not cause any difference in the reconstructed image quality. Therefore, Eq. (9.16a) can be simplified as follows:

Equation (9.16) is considered as differential azimuth compression (DAC) because it not only performs the spatial shift of Δx_n on s_c(x, k_y) but also provides a shift quantity of that depends on the azimuth value of k_y. This technique, originally developed by Stolt for seismic data processing, is called Stolt interpolation (or approximate Stolt interpolation in other textbooks).

Alternatively, the Stolt interpolation can be derived in a simple way as

FIGURE 9.3 Block diagram of Stolt interpolation algorithm.

Equation (9.15) can then be expressed as, assuming |P(k)| = 1.

where k′ = 2k and Again, the parameter k′ can be considered as a stretch of the k-axis; while k₀ is a constant phase and depends on both k_u and k. Since k = k_c + k_b and k_c k_b; therefore This is the same results as derived by Stolt interpolation.

Figure 9.3 shows the block diagram representation of the Stolt interpolation processing of SAR data.

The Stolt interpolation technique will be used to process some synthesized image data, employing the same radar and signal parameters as listed in Section 8.2. The broadside SAR with six targets is described first, followed by the squint SAR system.

9.3.1 System Model of Broadside SAR with Six Targets

A 2D system model of a broadside SAR with six targets is depicted in Fig. 9.4, where the six targets are located at (R_0a, y₁), (R_0b, y₂), (R_0c, y₂), (R_0d, y₂), (R_0e, y₂), and (R_0a, y₃), respectively. The vertical axis represents the azimuth lines (or cross-range samples), and the horizontal axis shows the range that is perpendicular to the flight path. Lsa, Lsb, Lsc, Lsd, and Lse are the synthetic aperture lengths corresponding to ranges R_0a, R_0b, R_0c, R_0d, and R_0e, respectively. The total synthetic aperture length that covers the six-target region is L_tot = Lsa + (y₃ – y₁).

FIGURE 9.4 System model of a six-target broadside SAR.

The received echo signal array based on Fig. 9.4 is displayed in Fig. 9.5. The horizontal axis is changed from time samples, as discussed in Chapter 8, to slant range samples. The slant range sample is shown with index n, where n = r/R_s with r and R_s as the slant range and slant range sample spacing, respectively. Normally, the curved line will occur at the beginning and end of the received signal array due to slant range differences at various radar positions. The differences are negligible for broadside SAR and are replaced with straight lines instead. The slant range is therefore considered approximately equal to range in the broadside case.

FIGURE 9.5 Received signal array based on Fig. 9.4.

The symbols <r1>, <r2>...<r6> are the received signal arrays corresponding to targets 1, 2,..., 6, respectively. Naza, Nazb, Nazc, Nazd, and Naze are the number of samples within the associated synthetic aperture length, which corresponds to ranges R_0a, R_0b, R_0c, R_0d, and R_0e, respectively. The parameters D₁, D₂,..., D₅ are defined as follows:

D₁ = integer of {[Lsa – (y₃ – y₁)]/A_s}

D₂ = integer of {[(y₂ – y₁) – (Lsb2 – Lsa2)]/A_s}

D₃ = integer of {[(y₂ – y₁) – (Lsc2 – Lsa2)]/A_s}

D₄ = integer of {[(y₂ – y₁) – (Lsd2 – Lsa2)]/A_s}

D₅ = integer of {[(y₂ – y₁) – (Lse2 – Lsa2)]/A_s}

Here A_s is the sample spacing along the azimuth axis (or m axis), and Lsx2 equals half of Lsx with x equal to a, b, c, d, and e, respectively. The size of total array samples along the vertical axis (or m axis) is 1 to 2Naza – D₁. N_Rb, N_Rc, N_Rd, and N_Re are the sample differences between R_0a and R_0i for i = b, c, d, e. They are defined as follows:

N_Rb = integer of ((R_0b – R_0a)/R_s)

N_Rc = integer of ((R_0c – R_0a)/R_s)

N_Rd = integer of ((R_0d – R_0a)/R_s)

N_Re = integer of ((R_0e – R_0a)/R_s)

The array size along the horizontal axis (or n axis) is 1 to N_Re + N_r. The sizes of signal arrays for <r1> to <r6> are listed below:

9.3.2 Synthesis of Broadside SAR Data Array

The following target location parameters are used to generate the received data array <ri>, i = 1, 2,..., 6 for the broadside SAR system.

y₁ = 0 m

y₂ = 100 m

y₃ = 150 m

R_0a = 7500 m

R_0b = 7650 m

R_0c = 8000 m

R_0d = 8350 m

R_0e = 8500 m

In terms of the given parameters of target locations and radar signals, the following data can be computed:

Lsa = 225 m

Lsb = 229.5 m

Lsc = 240 m

Lsd = 250.5 m

Lse = 255 m

Naza = Lsa/A_s = 563

Nazb = Lsb/A_s = 573

Nazc = Lsc/A_s = 601

Nazd = Lsd/A_s = 627

Naze = Lse/A_s = 637

D₁ = integer of ((Lsa – (y₃ – y₁)) /A_s) = 187

D₂ = integer of ((y₂ − y₁) − (Lsb2 − Lsa2))/A_s = 244

D₃ = integer of ((y₂ − y₁) − (Lsc2 − Lsa2))/A_s = 231

D₄ = integer of ((y₂ − y₁) − (Lsd2 − Lsa2))/A_s = 218

D₅ = integer of ((y₂ − y₁) − (Lse2 − Lsa2))/A_s = 212

Number of azimuth lines M = 2 × Naza – D₁ = 939

Number of range samples N = N_r + N_Re = 281

The waveforms of the real part of the synthesized signal arrays <r1>, <r2>,..., <r6> are shown in Fig. 9.6. From left to right, Figs. 9.6a and 9.6b are <r1> and <r2>, Figs. 9.6c and 9.6d are <r3> and <r4>, and Figs. 9.6e and 9.6f are <r5> and <r6>. The vertical axis represents the number of azimuth lines, where 1 out of every 30 azimuth lines is sequentially displayed throughout the total synthetic aperture length. The horizontal axis represents the range samples, with sample spacing of 10 m each and corresponding to a 30-MHz sampling rate.

Let {N_ai: N_af, N_ri: N_rf} be the size of a signal array, where N_ai is the starting sample number and N_af is the ending sample number along the azimuth axis (or y axis), while N_ri is the starting sample number and N_rf is the ending sample number along the range axis (or x axis).

The sizes of signal arrays for <r1> to <r6> can be computed as follows:

FIGURE 9.6 Waveforms of the real part of individual echo signal based on Fig. 9.4.

The overall received signal array s(t, u), which is the combination (or sum) of six signal arrays <r1>, <r2>,..., <r6>, is depicted in Fig. 9.7. The real part is in Fig. 9.7a; the imaginary part, in Fig. 9.7b. The size of the overall received signal array is {1:939, 1:281}. Accordingly, the FFT size used to process this time-domain signal array will be chosen as N_y = 2048 and N_x = 512.

9.3.3 Simulation Results

The data array s(t, u), represented by Eq. (9.1) and shown in Fig. 9.7, is first 2D Fourier-transformed to become S(ω, ω_D) as shown in Eq. (9.3). The range reference function P*(ω), which is the complex conjugate of the transmitting function P(ω), is then multiplied by the data array S(ω, ω_D) to become the range-compressed data array S_1c(ω, ω_D). The range-compressed data are then inverse-Fourier-transformed to become s_1c(t, ω_D). Figure 9.8 displays the range-compressed data array s_1c(t, ω_D) in the range–Doppler frequency domain. The origin of the Doppler frequency spectrum is at bin 1024. The pulse repetition frequency (or sampling frequency) along the azimuth axis is f_PRF, and equals 500 Hz in this example.

FIGURE 9.7 Waveforms of received signal based on Fig. 9.4.

Figure 9.9 is the top view of s_1c(t, ω_D) in the range–Doppler frequency domain. Only 1 out of every 50 rows along the Doppler frequency bins is displayed. As can be seen, the range-compressed signal has five columns of spectrum array along the range axis. Each spectrum array extends uniformly (about 80%) along the Doppler frequency axis. The first column array is located at range sample 181, which corresponds to the number of samples within the LFM pulse duration. The second column array appears at range sample 196 and is 15 range samples away from the first column array. Given the range sample spacing R_s = 10 m, the 15 range samples equal 150 m. The third column array appears at range sample 231, and is 50 range samples away from the first column array, which corresponds to a range distance of 500 m. The fourth column array appears at range sample 236, has a range sample number difference of 85 from the first column array, and corresponds to a range distance of 850 m from the first column array. The fifth column array appears at range sample 281 and is 1000 m away from the first column array with a range sample number difference of 100. Notice that the Doppler frequency spectrum appears similar for column arrays 2–5, but column array 1 seems different from the remaining four column arrays. This is because two targets contribute to the Doppler frequency spectrum on column array 1, and their phases cancel each other at some frequency bins.

FIGURE 9.8 3D view of s_1c(t, ω_D) in range–Doppler frequency domain.

FIGURE 9.9 2D view of s_1c(t, ω_D) in range–Doppler frequency domain.

Since the simulated targets are located between x = 7500 m and x = 8500 m, the center of the target locations is chosen at X₀ = R_0c = 8000 m. The ideal target location used by the reference function (or matched filter) for azimuth compression is then chosen at (8000, 0). The corresponding azimuth reference function can be synthesized in two ways.

1. The 1D azimuth reference function in s (slow time) domain, given the f_Dc (= 0 for broadside case) and frequency changing rate β, can be designed as

   

   where Δs = 1/f_PRF is the pulse repetition interval and – Nazc2 ≤ m ≤ Nazc2. The sample number Nazc corresponds to the reference point located at X₀ = R_0c = 8000 m. The reference function H_az(ω_D) can be obtained by taking FFT on h_az(s), with s = mΔs.

2. The 2D azimuth reference function can be obtained by substituting X₀ = 8000 in Eq. (9.10) and discarding the range matched filter function P*(ω):

   

   Notice that the frequency ω is in passband; that is, ω_c – ω₀ ≤ ω ≤ ω_c + ω₀ with ω_c as the carrier frequency and 2ω₀ as the bandwidth of the transmitting signal. The frequency ω_D corresponds to the true Doppler frequency.

Since no range migration correction is needed in broadside case, the 2D azimuth reference function H_az(ω, ω_D) will be approximated by 1D function H_az(ω_D), while both 1D and 2D reference functions will be applied to the squint SAR case later.

The azimuth compression is performed by multiplying every column of the range-compressed data array s_1c(t, ω_D) by the reference function H_az(ω_D) to become the 2D compressed signal s_c(t, ω_D). Since only one azimuth reference function is used to perform the azimuth compression on the whole data array, the result is roughly compressed along the azimuth direction. Figure 9.10 is a 3D view of the target function after taking the inverse FFT on the roughly compressed signal s_c(t, ω_D).

Figures 9.11 and 9.12 display two side views of Fig. 9.10, one from the range direction into the targets and the other from the azimuth direction into the targets, respectively.

Figure 9.11 shows that the targets appear as pulses located at five different ranges. The center one is located at range sample 231, and corresponds to the reference point located at x = X₀ = 8000 m. The large magnitude of this narrow pulse proves that it is accurately compressed in both range and azimuth directions. The other four pulses, appearing at range samples 181, 196, 266, and 281, have smaller magnitudes and are the targets roughly compressed in the azimuth direction.

FIGURE 9.10 3D view of the roughly compressed six-target function.

Figure 9.12 shows three pulses located at different positions along the azimuth axis. The center pulse is at azimuth sample 832, and corresponds to the reference target for azimuth compression. This position corresponds to the reference function located at y = 0. The other two pulses appear at azimuth samples 582 and 958, and correspond to targets located at y = −100 m [(582 – 832) = −250 samples away] and y = 50.4 m [(958 − 832 = 126 samples away], respectively.

FIGURE 9.11 Side view, from the range direction, of Fig. 9.10.

FIGURE 9.12 Side view, from the azimuth direction, of Fig. 9.10.

The roughly compressed data array s_c(t, ω_D) will now be processed by the Stolt interpolation, which improves the image quality through the differential azimuth compression (DAC) technique. The DAC technique is based on the implementation of Eq. (9.16b), which requires computation of two items. The first item is the range difference Δx_n = x_n – X₀, where x_n is the range value along the range axis and X₀ = 8000 m corresponds to range sample 231 as shown Fig. 9.11. For n = 1, 2,..., N_x along the range axis with N_x = 512 in this example, the range difference can be computed as

Δx_n = (n − 231) R_s.

The second item is where k_c = ω_c/c and k_y = (2πmf_PRF/VN_y). Here m = 1, 2,..., N_y and N_y = 2048 in this example. For a given value of m along the Doppler frequency axis, one can compute k_y first, then multiply every range sample of s_c(t, ω_D) by

After the Stolt interpolation process, the final target locations can be obtained by taking the IFFT with respect to ω_D. The more accurately compressed and reconstructed target function is obtained and shown in Fig. 9.13.

By comparing the two plots in Figs. 9.10 and 9.13, the difference is obvious. The new plot with DAC improves greatly in both magnitude and sidelobe compression of the target pulses.

FIGURE 9.13 3D view of refocused six-target function.

Figures 9.14 and 9.15 display two side views of Fig. 9.13. The first view is from the range direction into the targets, and the second is from the azimuth direction into the targets. The magnitudes of the compressed pulses, located at the two sides of the center one, reach about 130 in Figs. 9.14 and 9.15. In Figs. 9.11 and 9.12, their magnitudes are less than 40.

FIGURE 9.14 Side view, from the range direction, of Fig. 9.13.

FIGURE 9.15 Side view, from the azimuth direction, of Fig. 9.13.

The principles of synthesizing radar image data and the radar signal parameters used in broadside SAR will be applied to the squint SAR system. In addition, a nonzero squint angle θ_q and different specification of target locations will be incorporated into the squint SAR system. These are described next.

9.3.4 System Model of Squint SAR with Six Targets

The 2D system model of a squint SAR with six targets and a squint angle θ_q is shown in Fig. 9.16. Here the vertical axis represents azimuth samples (or lines), and the horizontal axis represents range. The six targets are located at (R_0a, y₁), (R_0b, y₂), (R_0c, y₃), (R_0d, y₄), (R_0e, y₅), and (R_0a, y₆), respectively. The synthetic aperture length labeled as Lsa, Lsb, Lsc, Lsd, and Lse corresponds to targets located at ranges R_0a, R_0b, R_0c, R_0d, and R_0e, respectively.

The radar signal array based on the system model described above is displayed in Fig. 9.17. The horizontal axis is changed from the range to slant range sample index n with n = r/R_s, where r is the slant range. The vertical axis is changed from the azimuth position to azimuth lines (or cross-range samples). The curved lines, which appear at the beginning and end of the radar pulse duration, are caused by the slant range difference during the time period when the target is illuminated by the radar beam. The slant range sample difference at two ends of the curved line is represented as Δ_x with x = a, b,..., e. The symbol Δ_a is defined as,

FIGURE 9.16 System model of a 6-target Squint SAR.

Other symbols, Δ_b, Δ_c, Δ_d, and Δ_e, can be computed accordingly.

The symbols <r1>, <r2>,...,<r6> are the received signal array corresponding to targets 1, 2,..., 6, respectively. Naza, Nazb, Nazc, Nazd, and Naze are the number of samples within the associated synthetic aperture length. The parameters D₁, D_ae, D_be, D_ce, and D_de are defined as follows:

FIGURE 9.17 Received signal array derived from Fig. 9.16.

D₁ = integer of {[Lsa − (y₆ − y₁)]/A_s}

D_ae = integer of [(Lse2 − Lsa2)/A_s]

D_be = integer of [(Lse2 − Lsb2)/A_s]

D_ce = integer of [(Lse2 − Lsc2)/A_s]

D_de = integer of [(Lse2 − Lsd2)/A_s]

The parameters N_Rb, N_Rc, N_Rd, and N_Re are defined in the same form as in the broadside SAR case.

9.3.5 Synthesis of Squint SAR Data Array

The following are the specifications of target locations used to generate the received data array:

R_0a = 7500 m

R_0b = 7650 m

R_0c = 8000 m

R_0d = 8350 m

R_0e = 8500 m

θ_q = 6°

y₁ = R_0a × tanθ_q = 788.3 m

y₂ = R_0b × tanθ_q = 804 m

y₃ = R_0c × tanθ_q = 840.8 m

y₄ = R_0d × tanθ_q = 877.6 m

y₅ = R_0e × tanθ_q = 893.4 m

y₆ = y₁ + 150 = 938.3 m

Given the target locations and the radar signal parameters, the following data can be obtained:

Lsa = R_0a × [tan(θ_q + 0.5 θ_H) − tan(θ_q − 0.5 θ_H)] = 227.5 m

Lsb = 232 m

Lsc = 242.6 m

Lsd = 253.3 m

Lse = 257.8 m

Naza = Lsa/A_s = 569

Nazb = 581

Nazc = 607

Nazd = 633

Naze = 644

D₁ = 193

D_ae = 37

D_be = 32

D_ce = 18

D_de = 5

N_Rb = 15

N_Rc = 50

N_Rd = 85

N_Re = 100

Δ_a = 2

Δ_b = 2

Δ_c = 3

Δ_d = 3

Δ_e = 3

Number of azimuth lines M = 2 × Naza – D₁ + D_ae = 982

Number of slant range samples N = c × T_p + (R_0e – R_0a)/R_s = 281

The real parts of the synthesized signal waveform corresponding to targets 1–6 are shown in Fig. 9.18. The vertical axis shows the number of azimuth lines, where 1 out of every 30 azimuth lines is sequentially displayed throughout the total synthetic aperture length. The horizontal axis shows the slant range samples, with a sample spacing of 10 m and corresponding to a 30 MHz sampling rate.

The sizes of signal arrays for <r1> to <r6> are shown below.

FIGURE 9.18 Waveforms of the real part of individual echo signal from Fig. 9.16.

The overall synthesized signal array s(t, u), which is the combination (or sum) of the six signal arrays <r1> to <r6>, is displayed in Fig. 9.19. The real part is shown in Fig. 9.19a; the imaginary part, in Fig. 9.19b. The size of the overall synthesized signal array is {1:982, 1:283}. Accordingly, the FFT size used to process this signal array remains the same as in the broadside SAR case; that is, N_y = 2048 and N_x = 512.

9.3.6 Simulation Results

The data array s(t, u), represented by Eq. (9.1) and shown in Fig. 9.19, is first 2D Fourier-transformed to become S(ω, ω_D). The range reference function P*(ω) is multiplied by the received baseband signal S(ω, ω_D) to become the range-compressed signal S_1c(ω, ω_D). The result is then inverse-Fourier-transformed to become s_1c(t, ω_D). Figure 9.20 is a 3D view of the range-compressed function s_1c(t, ω_D) in the range–Doppler frequency domain with frequency origin at bin 1. Comparing the spectrum with that in the broadside case, one can see that (1) the magnitude of the spectrum varies significantly along the Doppler frequency axis, and (2), the spectrum splits into two portions. As discussed previously, this is because the squint angle θ_q causes significant range migration, and the observed spectrum is the fallback version of the true spectrum of s_1c(t, ω_D).

Figure 9.21 is a 2D view of s_1c(t, ω_D) in the range–Doppler frequency domain. Only 1 out of every 30 rows is displayed. There are five columns of spectrum array in the display, and each array corresponds to one of the targets located at ranges R_0a, R_0b,..., R_0e (there are two targets contribute to the spectrum at location R_0a). In addition, every column array splits in two parts and occupies about 80% of the whole Doppler frequency range.

FIGURE 9.19 Waveforms of received signal from Fig. 9.16.

FIGURE 9.20 3D view of s_1c(t, ω_D) in range–Doppler frequency domain.

FIGURE 9.21 2D view of s_1c(t, ω_D) in range–Doppler frequency domain.

The first column array of the Doppler spectrum corresponds to targets 1 and 2 and tilts from range samples 181 to 185. The number 181 corresponds to the number of samples in the LFM pulse duration. The second column array of the Doppler spectrum tilts from range samples 196 to 200, which corresponds to a target that is 150 m away from the first column array. The third column array tilts from range samples 231 to 235, and corresponds to the target located 500 m away from the first column array. The fourth column array tilts from range samples 266 to 270, and corresponds to the target that is 850 m away from the first column array. The fifth column array tilts from range samples 281 to 285, and corresponds to a target that is 1000 m away from the first column array. The Doppler frequency spectrum looks similar for column arrays 2 to column 5, but column array 1 seems different from the other four column arrays. This is because two targets contribute to the Doppler spectrum on column array 1, and their phases cancel each other out at some frequency bins.

The matched filter for azimuth compression is chosen with a target located at the center of the target area, namely, (x, y) = (R_0c, 0) = (8000, 0). For image quality comparison purpose, both the 1D reference function H_az(ω_D) and 2D reference function H_az(ω,ω_D) will be designed as the azimuth matched filter to reconstruct the targets. The 1D azimuth filter H_az(ω_D) will first be used to reconstruct targets, and the 2D azimuth filter H_az(ω, ω_D) will be used later.

FIGURE 9.22 Synthesized 1D azimuth reference function for squint SAR system.

The azimuth matched filter H_az(ω_D), used in Chapter 8 for squint SAR processing based on the range–Doppler algorithm, is adopted here to serve as the 1D azimuth reference function. The frequency spectrum of H_az(ω_D) is shown again in Fig. 9.22 for convenience.

Azimuth compression is performed by multiplying every column of s_1c(t, ω_D) by the reference function H_az(ω_D) to become the 2D compressed signal s_c(t, ω_D). Figure 9.23 is a 3D view of the reconstructed target function after taking the inverse spatial DFT on the roughly compressed signal s_c(t, ω_D). The roughly reconstructed signal array has five targets located at the same y value along the azimuth direction as expected from the slant range viewpoint. The relative positions between the targets are not reconstructed correctly and will be corrected through the Stolt interpolation.

Two side views of Fig. 9.23, one from the range direction into the targets and the other from the azimuth direction into the targets, are displayed in Figs. 9.24 and 9.25. Figure 9.24 shows that the targets appear as pulses located at five different ranges. The center pulse is located at range sample 233, and corresponds to the azimuth reference function with range sample chosen to be x = X₀ = 8000. It appears to be accurately compressed in both range and azimuth directions. The other four pulses appear at range samples 183, 198, 233, and 283 and have smaller magnitudes. These four pulses are roughly compressed in the azimuth direction.

Figure 9.25 shows that the targets appear as two pulses located along the azimuth axis. The left pulse is at azimuth sample 626 and corresponds to the azimuth reference function located at y = 0. The right pulse appears at azimuth sample 1002. This pulse corresponds to 376 samples, or 150.4 m, away from the reference target.

FIGURE 9.23 3D view of roughly compressed target function.

FIGURE 9.24 Side view, from the range direction, of Fig. 9.23.

FIGURE 9.25 Side view, from the azimuth direction, of Fig. 9.23.

The Stolt interpolation algorithm requires computation of two parameters:

1. The range difference Δx_n = x_n − X₀. Given that X₀ corresponds to a range sample 233, the range difference between any range sample n and the reference target can be computed as Δx_n = (n − 233) R_s.
2. The ratio Here, k_c = (ω_c/c) and k_y is more complicated than in the case of the broadside SAR. As shown in Fig. 9.20, the Doppler frequency spectrum of the squint SAR system has two characteristics: (1) the true or original Doppler spectrum is a passband signal having a lower frequency edge f_DL and upper frequency edge f_DU; and (2) because of the natural digitization at sampling frequency f_PRF, the Doppler frequency spectrum folds from passband to baseband and splits into two parts to become the observable baseband frequency spectrum. The high-frequency part appears at the lower side of the baseband spectrum. Accordingly, the observed Doppler baseband frequencies f_DL and f_DU need to be adjusted to obtain their true values in the passband. Given squint angle θ_q, wavelength λ, and radar speed V, as shown in Eq. (6.16), the true values of f_DL and f_DU can be obtained. Let N_y be the FFT size used to transform the Doppler frequency spectrum, and N_yL and N_yU be the Doppler frequency bin numbers corresponding to the true (or passband) f_DL and f_DU. The corresponding true value of k_y can be computed as k_y = m Δk_y, where Δk_y = (2π f_PRF/VN_y) and m = N_yL to N_yU. Differential azimuth compression is then implemented by multiplying every range sample of s_c(t, ω_D) by

FIGURE 9.26 3D view of refocused target function.

After the Stolt interpolation, the final target function can be obtained by taking the inverse FFT with respect to the Doppler frequency. The result is shown in Figure 9.26.

By comparing Fig. 9.26 with Fig. 9.23, it can be seen that the new target function with differential azimuth compression improves greatly in magnitude with the target pulses. It also displays the correct target positions along the azimuth axis.

Two side views of Fig. 9.26, one from the range direction into the targets and the other from the azimuth direction into the targets, are depicted in Figs. 9.27 and 9.28, respectively.

Figure 9.27 shows that there are five pulse-like signals along the range axis, with the center signal located at range sample 233 and corresponding to the azimuth reference function. The magnitude of this center signal remains the same as the one in Fig. 9.24, while the magnitudes of the other four signals are improved from less than 30 in Fig. 9.24 to over 40 in Fig. 9.27.

Figure 9.28 shows that there are six pulses located at different azimuth lines (or y axes), which differ from the two pulses shown in Fig. 9.25. This demonstrates that the differential azimuth compression accurately compressed the targets and reconstructed the target locations correctly. The six pulses, corresponding to targets 1–6, are located at azimuth samples 495, 535, 626, 718, 758, and 871 respectively.

By setting target 1 at the position of y = 0, the sample number differences between target 1 and targets 2–6 are 40, 131, 223, 263, and 376, respectively. Given the azimuth sample spacing A_s = 0.4 m, the positions of y₂ to y₆ (corresponding to targets 2–6) can be computed and compared against the specifications listed below:

FIGURE 9.27 Side view, from the range direction, of Fig. 9.26.

FIGURE 9.28 Side view, from the azimuth direction, of Fig. 9.26.

Simulated	Specification
y1 = 0	0
y2 = 16 m	15.7 m
y3 = 52.4 m	52.5 m
y4 = 89.2 m	89.3 m
y5 = 105.2 m	105.1 m
y6 = 150.4 m	150 m

The second azimuth reference function H_az(ω, ω_D), based on Eq. (9.10) with X₀ = 8000 m, is now applied to process the same image data. By taking the inverse DFT on the 2D azimuth reference function H_az(ω, ω_D) with respect to the frequency domain, the range–Doppler frequency-domain view of H_az(t, ω_D) is shown in Fig. 9.29 with bin 1 as the Doppler frequency origin. As can be seen, the Doppler frequency spectrum of H_az(t, ω_D) is folded from passband to baseband and divided in two portions, the upper one representing the low-frequency band and the lower one, the high-frequency band. The thin-wall-like Doppler frequency spectrum is not equal in magnitude and aligns around range sample 440.

A different view of Fig. 9.29 is shown in Fig. 9.30, where 1 out of 50 rows along the Doppler frequency axis is displayed. As can be seen, the column-like Doppler frequency spectrum has two portions along the Doppler frequency axis and tilts to the left, opposite to the squint targets spectrum shown in Fig. 9.21. The column spectrum is located at range sample 442 and tilts toward the left at range sample 438 along the Doppler frequency direction.

FIGURE 9.29 3D view of H_az(t, ω_D).

FIGURE 9.30 2D view of H_az(t, ω_D).

The 2D azimuth reference filter H_az(ω, ω_D) is then multiplied by the 2D range compressed data array S_1c(ω, ω_D). The result becomes the roughly compressed data array. By using the same process of differential azimuth compression as in the 1D azimuth matched filtering case, the final target function can be obtained.

Figure 9.31 shows the reconstructed targets based on the 2D azimuth reference function. Compared with the results from the 1D azimuth matched filter, the 2D azimuth filter appears to have a better pulse compression ratio, and therefore better image quality.

Two side views of Fig. 9.31, one from the range direction into the targets and the other from the azimuth direction into the targets, are depicted in Figs. 9.32 and 9.33, respectively.

By comparing Figs. 9.32 and 9.33 with Figs. 9.27 and 9.28, one can see that the targets appear sharper in Figs. 9.32 and 9.33 than those in Figs. 9.27 and 9.28. The magnitudes of five pulses are near or above 70 in Figs. 9.32 and 9.33, while in Figs. 9.27 and 9.28, the pulse magnitudes are around 43 only.

From Figs. 9.32 and 9.33, one can see that the reference target, located at range sample 161 and azimuth sample 377, is different from that using a 1D azimuth matched filter. The range sample 161 is the result of the convolution of 2D reference function H_az(t, ω_D) and the 2D range-compressed data array s_1c(t, ω_D). The H_az(t, ω_D) peaks around range sample 442, while in the 2D range-compressed data array the reference target is located at range sample 231. The convolution of these two functions will have a peak at Mod[(442 + 231), 512] = 161. Similarly, azimuth sample 377 can be computed accordingly. The relative positions among other targets are quite accurate. The third (center and reference) target, located at range sample 161, is 50 range samples (50 × 10 = 500 m) away from target 1 (located at range 111) and targets 5 (at range 211). The distance between the targets 1 and 2 is 15 samples (15 × 10 = 150 m), which is also true for targets 4 and 5. Similarly, the locations of six targets along the azimuth direction are 244, 284, 377, 469, 509, and 620, respectively. With target 1 be assigned as y₁ = 0 and with azimuth sample spacing A_s = 0.4 m, the simulated results and original specifications are listed below:

FIGURE 9.31 3D view of a reconstructed target function.

FIGURE 9.32 Side view, from the range direction, of Fig. 9.31.

FIGURE 9.33 Side view, from the azimuth direction, of Fig. 9.31.

Simulated	Specification
y1 = 0	0
y2 = 16 m	15.7 m
y3 = 53.2 m	52.5 m
y4 = 90 m	89.3 m
y5 = 106 m	105.1 m
y6 = 150.4 m	150 m

9.4 RECONSTRUCTION OF SATELLITE RADAR IMAGE DATA

A satellite-based (RADARSAT-1) raw radar image file, generated by the Canada Space Agency and processed and distributed by MDA Geospatial Services Inc., will be used to reconstruct the radar image. The key parameters of this radar image data are listed below:

Antenna size (L × W): 15 × 1.5 m

Closest range between target and radar R₀ = 989,340 m

Radar image data array size:

Number of azimuth lines = 1536

Number of range samples = 2048

Effective radar velocity V_r = 7062 m/s

Radar signal:

f_c = 5.3 GHz

θ_H = λ/L = 0.0377 (radian)

θ_q = − 1.6°

Pulse duration time T_p = 41.74μs

Range FM rate α = 0.72135 MHz/μs

Pulse bandwidth: 30.111 MHz

Data format: baseband complex data, 4 bits each

A/D sampling frequency f_s = 32.317 MHz

Pulse repetition frequency f_PRF = 1257 Hz

On the basis of this information, the following data can be computed:

The complex baseband data, which appear in an in-phase–quadrature-phase format with 4-bit accuracy each, is first decoded and modified with automatic gain control provided with the data set. The real and imaginary part of the received signal, s_b(m, n) with m = 1, 2,...,1536 and n = 1, 2,...,2048, is shown in Fig. 9.34: Fig. 9.34a shows the real part and Fig. 9.34b the imaginary part, of the baseband signal. Only 1 out of 50 azimuth lines is shown in the plot.

FIGURE 9.34 Waveforms of the real and imaginary parts of a received satellite baseband signal. (With permission from MDA Geospatial Services.)

A range matched filter is required to perform range compression on the received raw data. Given the sampling frequency f_s, the chirp rate α, and the radar beamwidth, the range matched filter can be designed as

h_r(n) = exp [jπα (nΔt)²],

where Δt = 1/f_s is the sampling interval, α = 0.72135 MHz/μs, and −674 ≤ n ≤ 674.

Range compression is then performed by first transforming both the matched filter h_r(n) and every line of the range signal s_b(m, n) into the frequency domain using a 4096-point FFT. The product of the transformed signal and matched filter is then inverse-Fourier-transformed to become a time-domain range-compressed signal s_1c(m, n). The range-compressed signal s_1c(m, n), which appears in the range-azimuth (or time–spatial) domain, is shown in Fig. 9.35 for n = 1349–3396. Here the plot is shown in image format, with a gray level used to represent the magnitude of the signal (the white color corresponds to a large magnitude and the dark color corresponds to a small magnitude).

Since the sample length of the range matched filter is 1349, the range-compressed signal s_1c(m, n), which is obtained through 4096-point IFFT, will have range samples 1–4096. Because of the edge effect of circular convolution of IFFT, only range samples 1349–3396, corresponding to n′ = 1,2,...,2048, are chosen for further image compression.

The range-compressed signal array s_1c(m, n′), for n′ = 1,2,...,2048, is then transformed into range–Doppler frequency-domain signal S_1c(m, n). This is done by taking a 2048-point FFT on every column of the range–compressed signal. Since the sample length of the azimuth pulse is Ls/A_s = 664, and the number of azimuth lines of raw data is 1536, M_FFT = 2048 is chosen as the azimuth FFT length. Figure 9.36 displays an image of the range-compressed signal S_1c(m, n) in the range–Doppler frequency domain.

FIGURE 9.35 Image of a received satellite signal after range compression.

FIGURE 9.36 Image of a range-compressed signal in range–Doppler frequency domain.

The range-compressed data array S_1c(m, n) in the range–Doppler frequency domain serves as the basis for image reconstruction. This data array will first be processed in the wavenumber domain based on Stolt interpolation, followed by the range–Doppler algorithm.

The Stolt interpolation starts with a reference function serving as the azimuth matched filter to perform the bulk azimuth compression, followed by differential azimuth compression. The 2D azimuth matched filter is designed as follows:

Here X₀ = 1,000,000 m is the reference range, ω = 2π[f_c + (nf_s/2048)] and ω_D = 2πf_PRF(m + 2048(M − 1))/2048. The row and column variables m, n = 1,2,...,2048, and M = −5. The variable ω refers to the passband frequency of the time-domain signal, and ω_D is the true Doppler frequency along the azimuth direction.

The range-compressed signal S_1c(m,n), in range–Doppler frequency (t,ω_D) domain, is then transformed into (ω,ω_D) domain by taking a 2048-point FFT on every row of S_1c(m,n) to become S₂(m,n). The 2D matched filter H_az(ω,ω_D) is then applied on S₂(m,n) to become the roughly compressed signal S_2c(m,n) in the (ω,ω_D) domain. A roughly reconstructed image u(m,n) can then be obtained by taking 2D IFFT on the bulk compressed signal S_2c(m,n) with respect to ω and ω_D, and the result is shown in Fig. 9.37. As can be seen, the roughly reconstructed image provides quite detail and accurate map of the ground. Since the azimuth matched filter has sample length 664, and the data file along the azimuth direction is 1536, the edge effect of circular convolution will cause distortions on both ends of the image. Only azimuth lines ranging from 332 to 1867 will provide a correct image.

FIGURE 9.37 Radar image after bulk compression.

The bulk compressed signal S_2c(m, n) in the (ω,ω_D) domain is then inverse Fourier transformed into (t,ω_D) domain as S(m,n). It is then further processed by the differential azimuth correction, which requires computation of Δx_n and dk(m) and are defined as

where

k_fc = (2πf_c/c), M = −5, and m = N_DL,..., N_DU, n = 1,2,...,2048.

The reference sample N_rf = 1148 is chosen and corresponds to the slant range X₀ = 1,000,000 m. On the basis of the estimated Doppler frequency centroid f_Dc = − 6968 Hz and Doppler bandwidth = 940 Hz, the true upper and lower frequency band edges can be computed as f_DU = −7338 Hz and f_DL = −6398 Hz. The corresponding folded or observed baseband frequencies are f_Dc′ = −683 Hz, f_DU′ = − 1153 Hz and f_DL′ = −213 Hz. Given the Doppler sampling frequency f_PRF = 1257 Hz and N_FFT = 2048, differential azimuth compression is operated on S(m, n) from N_DL = 347 (corresponding to f_DL′) to N_DU = 1879 (corresponding to f_DU′), and for the whole range samples from 1 to 2048. That is, for every Doppler frequency column, only Doppler frequency bin numbers from 347 to 1879 will be adjusted with a phase factor of exp [− jdk(m)Δx_n]. In equation form,

U(m,n) = S(m,n) exp[− jdk(m)Δx_n]

After the differential azimuth compression, the reconstructed image is obtained by taking IFFT on (m, n). The result is shown in Fig. 9.38.

The image quality after differential azimuth compression is about the same as compared with the bulk compressed image shown in Fig. 9.37. This implies that the Stolt interpolation with bulk compression provides good quality of images for radar with steady moving speed V and stable squint angle θ_q, and for ground targets that cause smooth variation in the Doppler frequency centroid f_DC.

The range–Doppler algorithm will now be applied to the same image file. It operates on data in the range–Doppler frequency domain, namely, S_1c(m, n) as shown in Fig. 9.36. The sample length of azimuth matched filter equals N_az = R₀θ_H/A_s = 674, where R₀ = 1000,000 m, θ_H = 0.00377 radian, and A_s = V/f_PRF = 5.62 m are used. Given the slant range sample spacing R_s = 9.28 m, for a swath of 2048 range samples, the maximum slant range difference is 9.28 × 1148 = 10654 m, which is less than 1.1% of R₀. Therefore, the same azimuth matched filter h_az(m) can be used throughout the whole image frame without significant errors. The 1D azimuth matched filter is designed as follows:

FIGURE 9.38 Radar image after differential azimuth compression.

h_az(m) = exp [− j2πf_DCmΔs + jπβ(mΔs)²].

Here Δs = 1/f_PRF is the pulse repetition interval, β = 1780 Hz/s, f_Dc = −6968 Hz, and −332 ≤ m ≤ 332. The Doppler frequency-domain filter H_az(m), for m = 1,2..., 2048, is obtained by taking a 2048-point DFT on h_az(m).

The range migration of the squint SAR system is caused by the slant range difference at the edges of the 3-dB radar beamwidth. The amount of range migration can be computed, from Eq. (8.8b) as follows:

The integer part of ΔN_RK will be used for the range sample shift, while the fractional part will be used for interpolation on the range samples. In this example, the 16-set 8-tap sinc filter is adopted again as the interpolation filter. The maximum amount of range migration is 22 samples in this example.

After computing the range migration amount, every row of the range-compressed signal S_1c(m, n) along the Doppler frequency axis is then filtered with the corresponding sinc interpolation filter. The sinc-filtered output is then further adjusted by two elements. The first one is 4-sample delays caused by the 8-tap sinc filter. The second one is the range sample shift, which varies along the Doppler frequency axis. Every column of the adjusted output of the sinc filter is then multiplied by the azimuth reference filter H_az(m) to become U(m, n). The inverse DFT is then applied on U(m, n), and the result is shown in Fig. 9.39.

FIGURE 9.39 Radar image processed by range–Doppler algorithm.

By comparing Figs. 9.38 (or 9.37) and 9.39, it is difficult to tell which one has better image quality. Both the range–Doppler and Stolt interpolation algorithms generate good quality of radar images.

For comparison purposes, a different image data file was used and processed by the two algorithms described above. Figure 9.40 displays the image using the Stolt interpolation technique (without DAC), while Fig. 9.41 presents the image using the range–Doppler technique. Again, both algorithms appear to have similar image quality.

FIGURE 9.40 Radar image processed by Stolt interpolation technique.

FIGURE 9.41 Radar image processed by range–Doppler algorithm.

9.5 COMPARISON BETWEEN RANGE–DOPPLER AND STOLT INTERPOLATION ON SAR DATA PROCESSING

The difference between the range–Doppler and Stolt interpolation algorithms on processing SAR data can be summarized as follows: (1) the range–Doppler algorithm requires different range-dependent reference functions for azimuth compression, while the Stolt interpolation uses only one reference function for azimuth compression; (2) the range–Doppler algorithm uses a synthesized 1D azimuth reference function, while the Stolt interpolation utilizes a 2D azimuth reference function; (3) the range–Doppler algorithm corrects the range cell migration problem based on interpolation filtering of range samples for all Doppler frequency bins in the range–Doppler domain, while the Stolt interpolation algorithm corrects the problem through the multiplication of an azimuth and range-dependent phase factor, with respect to a reference point, on Doppler frequency samples for every range column in the range–Doppler domain.

The major computation requirements for range–Doppler and Stolt interpolation algorithms are based on the following general assumptions:

1. The basic computation requirement of a radix-2 FFT (or IFFT) are

   One complex multiplication (four real multiplications and two real additions)

   Two complex additions (four real additions)

2. The sample length of the range matched filter is N_fr.
3. The range matched filter is predesigned in the frequency domain and stored in memory as a table ROM (read-only memory)
4. The radar raw image data array is M_i × N_i.
5. The FFT (or IFFT) length used along the range or x axis is N ≥ N_i +N_fr − 1, and N is chosen to be a power of 2.
6. The reconstructed radar image data array is M_i × N_i.

In addition, the following assumptions are made for the range–Doppler algorithm:

• The sample length of the azimuth 1D matched filter is M_fa.
• The FFT (or IFFT) length used along the y axis is M ≥ M_i + M_fa − 1, and M is a power of 2.
• The interpolation filter coefficients are real numbers, and the filter length is N_fi.
• The group of interpolation filter is predesigned in the range domain and stored in memory as a tabular ROM.

The major computation requirements for range–Doppler processing algorithms are

1. Range compression
   1. FFT along x axis for N range samples
      1. Complex multiplications: 0.5M_iN log₂N
      2. Complex additions: M_iN log₂N
   2. Range matched filtering on M_i × N data array (azimuth–Doppler domain)
      1. Complex multiplications: M_iN
   3. IFFT along x axis for N range samples
      1. Complex multiplications: 0.5M_iN log₂N
      2. Complex additions: M_iN log₂N

      (After the IFFT operation, the range-compressed output data are in a range–azimuth domain, and the size of the data array is rescaled to M × N_i. The azimuth sample length is zero-padded to render M ≥ M_i +M_fa − 1 sample length for the following range cell migration correction and azimuth compression.)

2. Range cell migration correction
   1. FFT along y axis for M azimuth samples
      1. Complex multiplications: 0.5N_iM log₂M
      2. Complex additions: N_iM log₂M
   2. Fractional interpolation and sample shift
      1. Complex multiplications:* 0.5MN_i N_fi.
      2. Complex additions:* MN_i(N_fi –1.5)

         (* Since the interpolation filter coefficients are real, only two real multiplications are needed instead of 4 multiplications and two additions.)

3. Azimuth compression
   1. Azimuth matched filtering on M × N_i data array
      1. Complex multiplications: MN_i
   2. IFFT along y axis for M azimuth samples
      1. Complex multiplications: 0.5N_iM log₂M
      2. Complex additions: N_iM log₂M

      The data array is rescaled to M_i × N_i.

In summary, the total computation requirements for the range–Doppler algorithm are

Complex multiplications: M_iN (log₂N + 1) + N_iM(log₂M + 1 + 0.5N_fi)

Complex additions: 2M_iN log₂N +N_iM (2log₂M + N_fi − 1.5)

For the Stolt interpolation algorithm, two additional assumptions are made:

1. The 2D matched filter is in the (ω, ω_D) domain, and its array size is M × N.
2. The FFT (IFFT) length used along the y axis is M ≥ M_i + M_fa − 1, and M is a power of 2. Here M_fa is the azimuth sample length within the synthetic aperture length.

The major computation requirements for the Stolt interpolation algorithm are

1. Range compression
   1. FFT along x axis for N range samples
      1. Complex multiplications: 0.5M_iN log₂N
      2. Complex additions: M_iN log₂N
   2. Range matched filtering on M_i × N data array
      1. Complex multiplications: M_i × N

   (The range-compressed output data are in the azimuth–frequency domain with array size M_i × N. The data array is rescaled to M × N by zero-padding the azimuth sample to render M ≥ M_i + M_fa − 1 with M equal to a power of 2.)

2. Rough azimuth compression
   1. FFT along y axis for M azimuth samples
      1. Complex multiplications: 0.5NM log₂M
      2. Complex additions: NM log₂M
   2. 2D Matched filtering for M × N data array
      1. Complex multiplications: MN
   3. IFFT along x axis for N range samples
      1. Complex multiplications: 0.5MN log₂N
      2. Complex additions: MN log₂N

   (The 2D bulk-compressed output data are in the range–Doppler domain with array size M × N. The data array is rescaled to M × N_i, which corresponds to the effective values of 1D convolution between the 1D range matched filter and 2D data.)

3. Differential azimuth compression
   1. Complex multiplications: MN_i
   2. IFFT along y axis for M azimuth samples
      1. Complex multiplications: 0.5N_iM log₂M
      2. Complex additions: N_iM log₂M

   (The 2D differential azimuth-compressed output data are in a range–azimuth domain with an M × N_i array. The data array is rescaled to M_i × N_i, which corresponds to the effective values of 2D convolution between the 2D matched filter and 2D data.)

In summary, the total computation requirements for the Stolt interpolation algorithm are

For comparison purposes, letting the input image array be 1536 × 2048 (M_i = 1536 and N_i = 2048), the intermediate array M × N be 2048 × 4096, and the interpolation filter be an 8-tap filter, one obtains the following data:

For the range–Doppler algorithm:

Number of complex multiplications: 1.42 × 10⁸

Number of complex additions: 2.58 × 10⁸

For the Stolt interpolation without DAC:

Number of complex multiplications: 1.42 × 10⁸

Number of complex additions: 2.56 × 10⁸

For the Stolt interpolation algorithm:

Number of complex multiplications: 1.68 × 10⁸.

Number of complex additions: 3.0 × 10⁸.

Since the range–Doppler algorithm and the Stolt interpolation algorithm process the M_i × N_i image data array in the same spatiotemporal and frequency–Doppler frequency domains, both algorithms require roughly the same memory capacity for SAR processing. Taking the double-buffering requirement for the FFT/IFFT process, the range cell migration correction, and the differential azimuth correction, the total memory required for SAR processing can be estimated at about 2.5MN complex samples. For a data array of M × N = 2048 × 4096, the memory requirement is 2 × 10⁷ complex samples.