Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Theory and application of the momentary fourier transform 1998

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_1998-0342.pdf [ 7.56MB ]
Metadata
JSON: 1.0065213.json
JSON-LD: 1.0065213+ld.json
RDF/XML (Pretty): 1.0065213.xml
RDF/JSON: 1.0065213+rdf.json
Turtle: 1.0065213+rdf-turtle.txt
N-Triples: 1.0065213+rdf-ntriples.txt
Citation
1.0065213.ris

Full Text

THEORY AND APPLICATION OF THE MOMENTARY FOURIER TRANSFORM by Sandor Albrecht M . S c . E . E . Technical University of Budapest, Hungary, 1993 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F A P P L I E D S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F E L E C T R I C A L A N D C O M P U T E R E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 1998 © Sandor Albrecht, 1998 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract The discrete Fourier transform (DFT) is a widely used tool in signal or image processing and its efficiency is important. There are applications where it is desirable to use relatively small, successive, overlapped DFTs to obtain the spectrum coefficients. The momentary Fourier transform (MFT) computes the D F T of a discrete-time sequence for every new sample in an efficient recursive form. In this thesis we give an alternate derivation of the M F T using the momentary matrix transform ( M M T ) . Recursive and non-recursive forms of the inverse M F T are also given, which can provide efficient frequency domain manipulation (e.g. filtering). Discussion on the properties and examples of the usage of the M F T is given, followed by a survey on its efficiency. In this work we investigate the applicability of the M F T to synthetic aperture radar (SAR) signal processing, and in particular show what advantages the M F T algorithm offers to the SPECtral ANalysis ( S P E C A N ) method and burst-mode data processing. In the S P E C A N algorithm, the received signals are multiplied in the time domain by a reference function, and overlapped short length DFTs are used to compress the data. The azimuth F M rate of the signal varies in each range cell, which leads to the issue of keeping the azimuth resolution and output sampling rate constant. After the introduction to S P E C A N , we show what advantages and disadvantages the M F T has compared to the F F T algorithms. i i When a S A R system is operated in burst-mode, its azimuth received signal has a segmented frequency-time energy in its Doppler history. It requires that EDFTs be located at specific points in the spectral domain to perform the azimuth signal compression. After the introduction of the burst-mode data properties, we show why the short EFFT (SIFFT) algorithm has the requirement of arbitrary-length, highly-overlapped EDFTs to process burst-mode data, in which case the I M F T is shown to have computational advantages. 111 Contents Abstract ii Contents iv List of Tables vii List of Figures viii Acknowledgements x 1 Introduction 1 1.1 Background 1 1.2 Thesis objectives and outline 4 2 Theory and Properties of the Momentary Fourier Transformation 6 * 2.1 Introduction 6 2.2 The Momentary Fourier Transformations derived form the Recursive Momentary Matrix Transformation 6 2.2.1 The Recursive Momentary Matrix Transformation 7 2.2.2 The diagonal form of the M M T 9 2.2.3 Inverse of the diagonalized M M T 11 2.2.4 Momentary Fourier Transform ,. 12 2.2.5 The non-recursive Inverse M F T 14 2.3 Properties of the M F T 16 2.3.1 Cosine windows using M F T ; 16 2.3.2 Implementation of the M F T algorithm , 19 2.3.2.1 Software Implementation of the M F T 19 iv 2.3.2.2 Hardware Implementation of the M F T 21 2.3.3 Example of M F T Usage 25 2.4 Computing Efficiency of M F T 27 2.4.1 Arithmetic of M F T 27 2.4.2 Comparison of M F T to F F T algorithms 28 2.4.3 Advantages and Uses of M F T 33 3 Overview of SAR Processing 35 3.1 Introduction 35 3.2 Ideal point-target model 35 3.3 S A R signal compression 39 4 Application of MFT to SPECAN SAR Processing Algorithm 42 4.1 Introduction 42 4.2 The S P E C A N Algorithm 42 4.3 Mul t i look processing in S P E C A N 48 4.4 The S P E C A N Algorithm Using the M F T 50 5 Application of MFT to Burst-mode SAR Data Processing 62 5.1 Introduction 62 5.2 Burst-mode S A R processing 62 5.3 Properties of fully exposed targets in burst-mode processing 64 5.4 Properties of partially exposed targets in burst-mode processing 67 5.5 The SLFFT Algorithm 69 5.5.1 Number of good output targets of a single EFFT 72 5.5.2 Real data simulation of burst-mode processing 74 v 5.6 Efficiency of SJFFT using the I M F T and the EFFT algorithms 77 5.6.1 Effect of varying S A R parameters and SNR/efficiency tradeoffs 78 5.6.2 Arithmetic of the SEFFT algorithm using the I M F T and the EFFT algorithms 80 6 Conclusions 89 6.1 Summary 89 6.2 Future work 91 Bibliography 93 vi List o f Tables Table 1 Memory requirement of the M F T 20 Table 2 Real operations in M F T for 7VC coefficients 28 Table 3 Resolution versus D F T length in S P E C A N for C-band satellite S A R 48 Table 4 Spaceborne and airborne S A R parameters for S P E C A N arithmetic calculation. 52 Table 5 Reduced and full M F T versus mixed-radix F F T 56 Table 7 Envisat swath parameters 78 Table 8 Min imum and maximum burst bandwidth of the seven swathes 79 Table 9 The length of the JTJFTs and the corresponding dSNR 84 List of Figures Figure 1 Windowing of the discrete-time function 7 Figure 2 Block diagram of the full M F T algorithm 22 Figure 3 Block diagram of the full M F T algorithm without the modulo-TV FIFO 23 Figure 4 Hardware structure of one M F T block 24 Figure 5 F S K signal analysis using M F T 25 Figure 6 Signal detection using M F T 26 Figure 7 Shift between D F T s when the M F T is more efficient 29 Figure 8 Arithmetic of M F T and Radix-2 F F T when qMvr =1 30 Figure 9 Arithmetic of M F T and Radix-2 F F T when qMrr = 4 30 Figure 10 Floating Point Operations of D F T algorithms 31 Figure 11 Fast Convolution with M F T and Radix-2 F F T 33 Figure 12 Synthetic aperture radar geometry 36 Figure 13 Deramping of multiple targets 43 Figure 14 Processing regions and the placement of successive D F T blocks in single look case 45 Figure 15 Division of the good output samples into looks and the location of the D F T operations in multi-look processing 49 Figure 16 Azimuth F M rate and the D F T length with varying range 51 Figure 17 Arithmetic of S P E C A N azimuth compression with different D F T algorithms 55 Figure 18 The arithmetic of S P E C A N when the M F T is more efficient 58 v i i i Figure 19 The arithmetic of the Range-Doppler algorithm and the S P E C A N algorithm using the M F T and F F T algorithms 58 Figure 20 The output sampling rate of the S P E C A N algorithm 60 Figure 21 Burst-mode operation in 2-beam ScanSAR case 63 Figure 22 Burst-mode processing of 16 fully exposed and evenly spaced targets in one range cell 64 Figure 23 Burst-mode processing of fully and partially exposed targets in one range cell 66 Figure 24 Effect of the circular convolution on the targets Doppler history 68 Figure 25 How minimum and maximum EDFT is placed to compress groups of targets from each burst 71 Figure 27 The Doppler history of real burst-mode data 76 Figure 28 Burst bandwidth and dSNR of IS1 swath 80 Figure 29 Arithmetic of the SIFFT algorithm when applied to the IS1 swath 86 Figure 30 Arithmetic of the SIFFT when it is applied to Envisat A P burst mode operation 88 ix Acknowledgements First of all, I would like to thank my Mother and Edina for providing consistent support and encouragement throughout my studies at U B C . Without their love and care this work would not have been completed. I would like to thank my supervisor, Dr. Ian Cumming, for his supervision and academic guidance, as well as providing the opportunity to continue my research of momentary Fourier transform in the field of synthetic aperture radar processing. I am grateful for the financial support provided for this research by the Natural Sciences and Engineering Council of Canada and by MacDonald Dettwiler and Associates. Finally, I would like to thank all of the members of the Radar Remote Sensing Group at U B C for providing a friendly and effective work environment. x Chapter 1 Introduction 1.1 Background Linear transformations, such as the discrete Fourier transform (DFT) are frequently used in digital signal processing, and their efficiency is very important. In applications where the D F T is applied to a signal, it is often desirable to use successive, possibly overlapping DFTs of smaller extent than the full length of the signal to obtain the spectrum coefficients. These transformations are normally off-line operations on blocks of data, requiring N samples of the signal before the transformation can be computed. The momentary Fourier transform (MFT) which is derived here is a method of computing the D F T of a sequence in incremental steps. It can be computed using an efficient recursive formula, and it is useful in cases where the detailed evolution of the spectra of a discrete series is wanted, and in cases where only a few Fourier coefficients are needed. The spectrum components of the M F T can be calculated independently and only one complex multiplication and two complex additions are needed to update each spectrum component. The inverse momentary Fourier transform (IMFT) is the dual of the M F T and 1 shares the same property, while the non-recursive form of the I M F T requires only additions to obtain a sample of the time sequence from its spectrum with N samples delay. The computational order of the M F T to update an TV-point D F T is N, a factor of log2N improvement over the radix-2 F F T algorithm i f all incremental results are needed. If only a sub-set of the transform domain components are needed, the computing load of the M F T can be further reduced, calculating only the coefficients of interest. The M F T does not rely upon on N being power of two to obtain its efficiency, in contrast to standard F F T algorithms. Uses of the incremental D F T were introduced by Papoulis in 1977 [1], and by Bitmead and Anderson in 1981 [5]. A detailed derivation of the momentary Fourier transform was given by Dudas in 1986 [6]. In 1991, L i l l y gives a similar derivation, using the term "moving Fourier transform", and uses the M F T for updating the model of a time-varying system [7]. In this thesis we further develop the theory of M F T , examine its applications and in particular, see what advantages it offers to synthetic aperture radar data processing. A synthetic aperture radar (SAR) is a powerful sensor in remote sensing which is capable of observing geophysical parameters of the Earth's surface, regardless of time of day and weather conditions [3]. S A R systems are extensively used for monitoring ocean surface patterns, sea-ice cover, agricultural features and for military applications such as in the 2 detection and tracking of moving targets. A S A R transmits radar signals from an airborne or spaceborne antenna which is perpendicular to the flight direction of the platform which is travels at a constant velocity. The back-scattered signal is collected by the antenna and stored in a raw format. Extensive signal processing is required to produce the output S A R image. The SPECtral ANalysis ( S P E C A N ) S A R processing algorithm was developed in 1979 by MacDonald Dettwiler and Associates, as a multi-look version of the deramp-FFT method of pulse compression. In S P E C A N , the received signals are multiplied in the time domain by a reference function, and overlapped short length DFTs are used to compress the data. In contrast, a precision processing algorithm such as the Range Doppler (RD) method requires both forward and inverse D F T operations, thus it is less computationally efficient. S P E C A N is an efficient algorithm for moderate to low resolution processing and generally implemented in quick look processors for viewing of magnitude detected imagery data. Burst-mode operation is used in S A R systems, to image wide swaths, to save power or save data link bandwidth. Several spaceborne remote sensing missions employ the ScanSAR mode in addition to other operational modes for radar imaging. Canada's Radarsat satellite, which was successfully launched in 1995, is a sophisticated Earth observation system developed to monitor environmental changes. The imaging platform supports various S A R operating modes, including a ScanSAR mode for the low- resolution (~100m) imaging of ground regions of width 500 km. 3 A n Advanced S A R ( A S A R ) system wi l l be flown on the Envisat satellite polar platform to be launched in 2000 by the European Space Agency. This system wi l l be able to operate in three burst-modes: alternating polarization mode (AP), wide swath mode (WS) and global monitoring mode (GM) . Alternating polarization mode provides high resolution in any swath with polarization changing from sub-aperture to sub-aperture within the synthetic aperture. This results in two images of the same scene in different polarizations combination with approximately 30 m resolution. In the wide swath mode the ScanSAR technique is used providing images of a wider swath (405 km) with medium resolution (150 m). 1.2 Thesis objectives and outline The objective of this research is to further develop the theory of the M F T , examine its properties and applications, and in particular, see what advantages it offers to S P E C A N processing and to the short IFFT (SIFFT) burst-mode processing algorithm. Chapter 2 presents the theory and properties of the momentary Fourier transform. Here, we introduce the recursive form of the momentary matrix transform ( M M T ) , and show when the M M T takes the form of the D F T or the IDFT, the resulting M F T and EvLFT have an efficient computational structure. The properties and computing efficiency of the M F T is also discussed in this chapter. 4 In Chapter 3, an overview of SAR processing is given, where the conventional compression method of the SAR signals is studied. This chapter gives a background knowledge for the research of the SPECAN algorithm and burst-mode data processing. The azimuth FM rate of the received signal varies in each range cell, which leads to the issue of keeping the azimuth resolution and output sampling rate constant. After the introduction of the SPECAN algorithm in Chapter 4, we show what advantages the MFT method offers vs. the FFT algorithms when they are applied to the SPECAN SAR processing algorithm. In Chapter 5, the ScanSAR operation mode will be introduced, and the received burst- mode data properties will be analyzed. After the effect of the varying SAR parameters and SNR/efficiency tradeoffs, a survey on the arithmetic of the SIFFT algorithm using IMFT and EFFT is given. Here, we show that the JJVEFT algorithm can improve the computational efficiency of the SIFFT algorithm in certain burst-mode data processing cases. Finally in Chapter 6, conclusions of the efficiency and applicability of the momentary Fourier transform to SAR processing will be drawn, and suggestions for possible future work will be given. 5 Chapter 2 Theory and Properties of the Momentary Fourier Transformation 2.1 Introduction In this Chapter, we give a derivation of the momentary Fourier transforms from the momentary matrix transform in Section 1.2. Section 1.3 gives a survey on the properties of the M F T , and in Section 1.4 a discussion on its computational efficiency is given. 2.2 The Momentary Fourier Transformations derived form the Recursive Momentary Matrix Transformation In this section, we introduce the matrix form of the momentary transform, and show that it has a recursive form. We also show that when the momentary matrix transform takes the form of the D F T or the inverse D F T , the resulting M F T has an efficient (recursive) computational structure. In the last part of the section, the inverse of the M F T is introduced, as well. 6 2.2.1 The Recursive Momentary Matrix Transformation Let xi be a sample of an arbitrary complex-valued sequence of one variable. The sequence wi l l be analyzed through an TV-point window, ending at the current sample i. In subsequent analyses, the window wi l l be advanced one sample at a time. A t sample i x, enters the window, while Xi.N leaves the window, as shown in Figure 1. Samples Figure 1 Windowing of the discrete-time function At samples i-1 and i, the windowed function can be represented by the following two column vectors: X, = v i - i v i - i X, (1) 7 Let T be an NxN non-singular matrix, which represent a linear transformation and has the inverse T" . The sequence of index vectors is transformed by T at each sample: • • • - y,-i =Tx1._„ y,. =Tx,., . . . (2) Let P be the NxN elementary cyclic permutation matrix. When the vector is multiplied by P, a one-element circular shift is performed, such that the index of each element is increased by one, and the first element becomes the last one: x, , where P= 0 1 . . 0 . 0 1 . 0 . . 0 1 . . . . 0 1 1 0 . . 0 (3) Using the result above, the x; vector can be expressed by the shifted xj.i vector and with an adjustment vector Ax\ made from the difference between the samples entering and leaving the window: Xi-(N-\) 0 + X i - l 0 Xi-N _Xi~Xi-N _ = P x w + A X i (4) 8 Substituting (4) into the transformation associated with the ith window in (2) and using the inverse transform XM = T"1 V M , the following relationships are obtained: y, =Tx, =T[Px,_ 1+Ax,.]=TPT- Iy,_ 1+TAx i (5) Equation (5) expresses the recursivity of the momentary matrix transforms (MMT), since calculation of the newly transformed index vector v* needs the previously transformed index vector and the difference between samples entering and leaving the window. 2.2.2 The diagonal form of the MMT The momentary matrix transform is particularly efficient and it can be calculated by components only i f the product of similarity matrix transform TFT"1 in (5) is diagonal. The P matrix has N distinct eigenvalues (Ao, ^N-I) which are the nth complex unit roots, Xk = w~k - ej2nkJN,k-0,1,2,...N-1. To each eigenvalue, N linearly independent eigenvector corresponds as follows: A 0 = w°=l <=> s0= w w w (N-l) 9 Xk=wk <=> st= w w -2k W (N-l)k ; - >• K-i= w ' ' <=> V i = w w 1 2(N-\; w •(N-l )(N-l) (6) If the eigenvectors are chosen to be the columns of the inverse of the T matrix, then TPT"1 is a diagonal matrix, with the eigenvalues of P along its diagonal: TPT 1 = S PS = S n S, s'1 s"1 3 N- 2 s N-l I I 0 -A, 0 0 0 A 2 0 0 0 . . A . (7) where S is the eigenvector matrix of P. The diagonalizing matrix S is not unique. A n eigenvector sk can be multiplied by a constant, and wi l l remain an eigenvector [2]. Therefore the columns of S can be multiplied by any nonzero constants and produce a new diagonalizing S. There is also no preferred order of the columns of S. The order of the eigenvectors in S and the eigenvalues in the diagonal matrix is automatically the same. Therefore, all T matrixes, which satisfy the above mentioned properties w i l l diagonalize the momentary matrix transform: 10 Xk 0 . . 0 0 A, 0 . 0 0 0 0 0 (8) where k, I, m e {0,1,...,N-l) and T N _i is the last column of the T matrix. 2.2.3 Inverse of the diagonalized MMT If y i is available at each sample and the columns of T are the eigenvectors of P , an efficient implementation of the inverse of the M M T can be obtained. The inverse MMT (IMMT) at time i: (9) 1 1 1 1 1 11 r y,.0 ^ XHNi) (10) The first row of T" 1 contains only ones (10), so the oldest element of X ; can be computed using adds only: AM i - ( N - l ) k = 0 (11) 11 from which the elements of the input sequence (JC , .^.;j.. .JC,) can be computed from the transform domain sequence y; with N-l sample delay. 2.2.4 Momentary Fourier Transform The matrix of the Discrete Fourier Transform (DFT) and the Inverse Discrete Fourier Transform (EDFT) have the properties described in Section 1.1.2, thus their columns are the eigenvectors of the matrix P: DFT = F = S w 2 w w 2iN-\) W 1 N-l W 2-(N-l) W (N-\HN-\) (12) IDFT= F 1 = S = N w 1 1 „-<N-\) W w-iN-l, w-l<N-l> w -Z(N-l) lHN-\HN-l) _1_ N 1 N - l w 2(N-l) W JN-lXN-l) W W 2 W „2(N-l) (13) 12 Using that w is the Nth complex unit root (i.e. w'k = wN'k), it can be seen that the columns of the IDFT matrix are the same as the D F T matrix, but they are in reverse order from the second column (13). Therefore, i f T performs the D F T (14) or the IDFT (15), diagonal forms of the M M T can be obtained: y,.= F P F - ' y ^ + F A x , = 1 0 0 wA 0 w'2 0 0 0 w x ^ F - ' P F x ^ + F - ' A y , . = 1 0 0 w1 0 w 2 0 0 0 . . 1 0 0 w(NA) 0 w w 0 0 0 N-l {N-l) Q 0 0 0 0 y,-. + 0 -(N-l) I w2 N-l w 0 " 0 . 0 -1 w w w -2 W -(N-l) ( X l - X i - N ) (14) (y,•-y,,A,)= w w 1 -(N-l) (N-2) W ( y i - y i - s ) (15) Equation (14) expresses the recursive equation of the momentary Fourier-transform (MFT). The TV-element vector y; contains the Fourier coefficients of the TV-point sequence xj ending at sample i. Note that each spectral component can be calculated independently, 13 y = w i.k k (yi-u  + xi-x, ) (16) which increases efficiency i f only a few frequency components need to be computed, as in the zoom transform. On the other hand, equation (15) is the dual pair of the M F T , the recursive inverse momentary Fourier-transform (IMFT), where the TV-element vector x; contains the TV- point time sequence and y-, contains TV Fourier coefficients ending at frequency bin /. Note that the each sample in Xi can also be obtained independently and the same twiddle factors, but in different order, can be used to calculate both the M F T and I M F T . Thus it has been showed that i f the D F T or the EDFT performs the momentary matrix transform of a sequence the elements of the transformed sequence can be computed recursively and independently using TV complex multiplies and N+l complex adds (computational savings are available i f the input sequence is real-valued). 2.2.5 The non-recursive Inverse MFT The non-recursive inverse momentary Fourier transform can be expressed using (11) and (13) as follows: x i-(N-l) i,k (17) 14 from which each sample of the input sequence (XJ) can be computed using adds only from the spectrum (y0 with N-l sample delay. In this way the M F T - non-recursive I M F T transform pair can provide an efficient frequency-domain manipulation method (e.g. filtering), especially i f many of the D F T coefficients are not calculated. If the elements of x; are real, taking advantage of the conjugate symmetry of the spectrum, the oldest element can be computed using only the real part of the spectrum components: t-(N-l) N k=0 (18) It has been showed in [6] that i f X ; is real, the Hilbert transform of Xj.(N-i) can be obtained to sum only^the imaginary part of the spectrum components: 1 N-l r -, H{x 1 = -Hm{y.k} i-(N-l) N k=0 (19) In this case the MFT/non-recursive I M F T pair can be useful for different signal processing applications where the in-phase (I) and quadrature component (Q) of the signal is needed (i.e. communications and radar systems). 15 2.3 Properties of the M F T In this section, some properties of the M F T are given. Section 1.3.1 shows how to implement cosine windows in the frequency domain using the M F T . In Section 1.3.2, discussion on the software and hardware implementation of the M F T algorithm is given, while section 1.3.3 gives an example of the use of the M F T . 2.3.1 Cosine windows using MFT The TV-point D F T treats the data sequence as i f it has a periodicity of N samples, xi = Xi+kN, for all integer k. In practice many signals do not have the above periodicity. If the boxcar window is applied to such a signal, the D F T wi l l treat it i f there were discontinuities at its edges. Ringing effects near the edges of filtered signals may occur as a result of these spurious discontinuities [4]. Such effects can be reduced by applying a more appropriate window. In addition to selecting a portion of the input sequence, the window modifies this portion to make it continuous at the edges when regarded as periodically repeated. Several types of window have been described in the literature [1], [4]. This section introduces how the Planning, Hamming and Blackman window can be implemented using the M F T . Given the discrete-time sequence xj, we wish to calculate the M F T of the windowed data xwJ = Wj-xi at time i, where wi is the window function. The ith element of the window may be expressed as follows: Hanning: w,= 0.5 1-cos ( 2ni { N J 16 Hamming: w,= 0.54 - 0.46 cos ( 2ni \ Blackman: w.= 0.42 - 0.5 cos ( 2ni ^ v N , + 0.08 cos r Ani ^ v N j The derivation for the Blackman window is given below: 0.42 - 0.5 cos =0.42 JC, - 0.25 0.04 exp exp ( Ani\ J V N ( 2ni { N 2ni \ N / f exp + 0.08 cos f Ani ^ v N y j + exp . Ani f 2ni ^ N x,- + N Taking the D F T of each part of (17) the spectrum of the windowed data at time i: y = 0.42y - 0.25 wi, k i , k y.M + yu-i + 0.04 (20) (21) (22) Therefore, the M F T of the windowed data can be obtained simply by maintaining the M F T of the non-windowed data and applying a weighted average in the spectrum. Similar results can be easily derived for the other cosine windows. These are: Hanning: y = 0.5 yik - 0.25 \yiM + y i t _, ] wt,k Hamming: y = 0.54 y - 0.23 Lv a + 1 +yik_l \ i,k wi,k (23) 17 Although only generalized cosine windows can be applied easily with the M F T , arbitrary windows can be approximated i f enough cosine terms are used. Note, that the memory requirement of the M F T algorithm gets larger, while its efficiency drops as the number of terms increases. The edge effect of the boxcar window can also be compensated i f the non-weighted moving average of the spectrum components is used. The moving average of a spectrum coefficient at times i, for L (L < N) consecutive samples can be expressed as: = 1 ' y mavg_i,k T y i,k L j=i-L+\ (24) It also has a recursive form where, the calculation of the averaged spectrum coefficient needs the previously calculated average and the difference between the spectrum coefficient entering and leaving the averaging window: y = y + mavg_i+l,k mavg_i,k ^ yt+\,k yi-L+\,k ^ (25) The above defined moving average with the long-term average (26) of the M F T coefficients can also be useful for statistical analysis of the input discrete sequence. 1 M yk = 7 7 s ^ ' M > > N k,avg M i=l (26) 18 2.3.2 Implementation of the M F T algorithm In this section, discussion on the software and hardware implementation of the M F T algorithm is given. Section 1.2.2.1 shows the computer coding of the M F T with its memory requirement, while the principle hardware structure of the M F T is given in Section 1.2.2.2. 2.3.2.1 Software Implementation of the M F T As it was shown earlier, the spectrum components in the M F T algorithm can be calculated independently from each other. Thus, the M F T can be built up from identical blocks, where a block refers to equation (10). The software implementation of one M F T block can be obtained using the trigonometric form of the equation: y.t =w~*(y,_ u+ *,.-*,•_„) i,k -j2nk . w ~ k - e N = cos (<&k }+- j sin ( ® k ) -jink where k N (27) Re /y } = cos(<Dj(Re/ ry J + Re{xrx.N J) - sin (Ok)(Im{y } + Im{xrXi_N}) i,k i-l,k Im{y } = cos (oj(lm/ r y J+lmfo-x^}) + sin(oJ(Re/> y" + R e 7 ) i,k i - l , k i-l,k (28) Equation (28) corresponds to the kth M F T block for complex xi, where is the kth spectrum component at sample i. Re{} means real part and Im{} means imaginary part of a complex number. 19 The M F T blocks can be organized in a for loop to calculate the needed D F T coefficients. The following pseudo-code segment illustrates the computer coding of the M F T algorithm, assuming the sine and cosine arrays (twiddle factors) have been precomputed: c a l c u l a t e ( x i - X i _ N ) ; for k = s t a r t to start+N c-l do MFTblock(k); endfor If the calculation of the spectrum coefficients is off-line, the difference of the entering and leaving samples of the window can be calculated for the whole data set and stored in a file or an array in the memory. If it is on-line, a modulo-A 7 array is needed to calculate xi-Xi-N. The index of the for loop in the pseudo-code indicates, that within the valid spectrum components only a smaller interval of the D F T coefficients can be computed. If the parameter start is zero and Nc = N , then all the spectrum coefficients are going to be calculated. Within the MFTblock procedure, the previously computed spectrum coefficient should be stored in an array for the computation of the recent one. The following table gives the memory requirement of the M F T algorithm: Array type Size Twiddle factors (sine and cosine arrays) 2-NVB bits Modulo-A^ FIFO for complex 2-N-B bits Spectrum coefficients at time /'-/ - y,./ L 2 N t B bits Spectrum coefficients at time i - y,.;^ 2-N c -B bits Table 1 Memory requirement of the M F T 20 In Table 1, parameter B is the number of bits used during the arithmetic operations, thus B = 32bits i f floating point arithmetic is used. B should be at least 24 bits for the fix point arithmetic, concerning the sensitivity of the M F T algorithm for the quantization error of the sine and cosine function. Note, the memory requirement of the M F T depends on the calculated spectrum coefficients. If the whole spectrum is computed, NB byte memory is needed for the computation. 2.3.2.2 Hardware Implementation of the M F T The trigonometric form of the M F T for one spectrum component (29) can be easily implemented in hardware. From the basic blocks of M F T a parallel hardware structure can be built for the computation of the D F T coefficients. Figure 2 illustrates the block diagram of the concurrent implementation of the M F T blocks of the full M F T algorithm. Note, the updating time of the fully concurrent implementation is equivalent to the propagation time of one M F T block, regardless of the number of the calculated spectrum coefficients. 21 Figure 2 Block diagram of the full M F T algorithm In Figure 2 the architecture of the full M F T contains a modulo-TV F I F O register to obtain Xi.N. If all the spectrum coefficients are computed, the leaving sample of the window at time i can be expressed using the I M F T algorithm: ^ N-l X i - N = T7 X y.,, N to "lk (29) Substituting (29) to (10), the recursive equation of one M F T coefficients becomes the following: 22 -k y -w yt-ik + xt (30) In (30), the data sample at time i and all the spectrum coefficients at time i-1 are needed to obtain yiik. The memory requirement of the M F T algorithm is reduced by 2-N-B bits, because there is no need to save the input data samples in a F IFO, while the arithmetic of the M F T increased by 2N-1 real operations due to the calculation of I M F T . The block diagram of the M F T corresponding to equation (30) is shown in Figure 3. x, - X,. MFT Block #0 M MFT Block #T Xi"(N-1) MFT Block #N-1 Yl.N.1 1/N Figure 3 Block diagram of the full M F T algorithm without the modulo-TV FIFO The detailed hardware structure of one M F T block for complex Xi is given in Figure 4. This implementation contains four multipliers ( M P Y ) and four adders ( A L U ) to obtain 23 the complex arithmetic of the M F T . The twiddle factors and the previously calculated spectrum coefficients are stored in registers. Because of the parallel computation of the real and the imaginary part of the M F T coefficients, the updating time of the spectrum coefficients is limited only by the propagation time of two multipliers and two adders. MFT Block #fc cos(k2»t/N) sin(k2n/N) MPY MPY MPY MPY A L U + + A L U Re{YM i k} + A L U + A L U + + Re{x,-x,.N} Re{Y,k) lm{YLk} lm{x,-x,.N} Figure 4 Hardware structure of one M F T block 24 2.3.3 Example of M F T Usage To illustrate the usage of the incremental form of the M F T , a frequency shift key (FSK) modulated sinusoidal signal of length 4N samples is used. Using an analysis window length 7V=100, and two frequencies of 5 cycles/window and 29 cycles/window, the magnitude of the evolving spectrum is shown in Figure 5, when the M F T is incremented by one sample at each analysis stage. FSK Signal 0 50 100 150 200 250 300 350 Time Time Dependent Spectrum of the FSK Signal Figure 5 F S K signal analysis using M F T The M F T begins with the initial conditions of yo = 0. This is equivalent to having N zeros precede the data vector. In Figure 5, note how the energy in the spectrum rises from zero to a maximum in the first /V samples. Also note how spectral leakage is observed in the 25 first N-l time samples, because the sinusoidal signal does not have an integer number of cycles/window over this time. At time N, there is an integer number of cycles/window, so all the energy in the spectrum lies in one bin. For the next A M samples, leakage occurs again as the window sliding towards to the next frequency component of the signal. The spectral energy of the 5th frequency bin decays to zero while the spectral energy of the 29th bin rises to its maximum. This spectrum energy 'swapping' between the two frequency bins is repeated as the window is moving through the two frequency components. In Figure 6, the same F S K signal is analyzed in the present of noise. The spectrum energy swapping between the two frequency bins is also noticeable, which shows how the M F T can be useful for signal detection in noise environment. FSK Signal in Noise Time Dependent Spectrum ot the FSK Signal Figure 6 Signal detection using M F T 26 2.4 C o m p u t i n g Efficiency of M F T In this section, a survey on the arithmetic of the M F T is given, followed by a discussion of its efficiency. The M F T is particularly efficient compare to F F T algorithms, when successive DFTs with high overlap ratio are to be computed or when only a few spectrum coefficients are needed. Examples of applications of the M F T to signal processing is also given, here. 2.4.1 Arithmetic of MFT The previously derived equation for one spectrum component of an N samples long M F T at time i: y., =w'k:(yM.*+*/ The twiddle factors (w"k) can be calculated only once and stored in an array before the M F T procedure. This computation is not included in the arithmetic of M F T . The difference of the sample entering and leaving the window - XJ-X;_N - can be pre- calculated at each time moment and used for the calculation of all spectrum coefficients. In this case, the spectrum of X ; can be updated from the spectrum of x\.\ using only N complex multiplies and N+l complex adds, i f x j contains complex-valued data. If Xj is real, N/2 complex multiplies and N/2+1 real adds are needed to obtain the N/2 new spectrum coefficients. Table 2 gives a summary of the number of real operations for these 27 cases when only Nc coefficients are calculated (N c <N for complex data and Nc <N/2 for real data). Input data Real Multiplies Real Adds Real Operations Real 4 N C 3NC+1 7N c +2 Complex 4 N C 4N c+2 8N c+2 Table 2 Real operations in M F T for Nc coefficients Note, the number of operations in each case can be reduced with one, i f the D C component is calculated, because in that case the twiddle factor equals to 1 (wk = 1 when k=0). 2.4.2 Comparison of MFT to FFT algorithms Consider the case where TV point DFTs are used to analyze an M-point complex-valued data record. If the window is shifted by q samples between each D F T application, where M-N 1 <q <N, then +1 DFTs are needed to spectrum analyze the record, in the case q of F F T . If the M F T is applied, M M F T s are needed, because the spectrum coefficients have to be calculated in each time samples, irrespectively of the value of q. Then, when radix-2 FFTs are used: OPS — f M - N ^ + 1 [5Mog 2 ( iV)] (31) 28 real operations, while in the case of M F T : OPSMFr=M[8Nc+2] (32) real operations are needed to analyze the whole record. From (31) and (32) the number of shift between DFTs when the M F T is more efficient than the radix-2 F F T can be expressed: (M -N)[5N\og2(N)] QMFT < M(SNc-l)-5N\og2(N) (33) A s we can see from (33), qMrr is function of the length of the data record (M), the size of the window (AO and the calculated M F T spectrum coefficients (Nc). In Figure 7, the shift between DFTs when the M F T is more efficient is shown as a function of the window length, with two values of M and Nc: Number of shift between DFTs when MFT is more efficient than FFT Number of shift between DFTs when MFT is more efficient than FFT 1 . . , . . , ! « 1 Total sample* an alyzed = 5000 i MFT i r - » i> i i | \ - i I:::*- .0 o FullMFT 0 ! j< |18 5 "16 t a 14 6 Total samples analyzed: 0••! • 100 200 300 400 500 600 Window size [sample] 700 800 900 1000 100 200 300 400 500 600 700 Window size [sample] 800 900 1000 (a) (b) Figure 7 Shift between DFTs when the M F T is more efficient 29 The full M F T is more efficient compared to the radix-2 F F T , i f the shift between D F T s is very small (qMFT ^ 5), while for the reduced M F T (Nc = N/4), the M F T is more efficient even for larger values of shift. Note, i f the data record is longer (Figure 7 (b)), the values of qMFT are larger for all window sizes. The computational load for small amount of shifts is illustrated in Figure 8 and 9: Arithmetic of MFT and Radix-2 FFT - Shift Between DFTs = 1 sample(s) Arithmetic of MFT and Radix-2 FFT - Shift Between DFTs = 1 sample(s) — i 1 1 1 1 1 1 1 1 .Total samplesanalyzed.?. 5QO0, 400 500 600 700 Window size [sample] Total samples ana|yzrf=M000 400 500 600 700 Window size [sample] (a) (b) Figure 8 Arithmetic of M F T and Radix-2 F F T when qMFT = 1 Arithmetic of MFT and Radix-2 FFT - Shift Between DFTs = 4 sample(s) Arithmetic of MFT and Radix-2 FFT - Shift Between DFTs = 4 sample(s) nples an 1 1 1 1 1 1 ilyzcd = flNM) 1 1 Radix-2 FFT / — \ ; Full MFT 1/4NAU FT J.4 - -0-- _ . • 1 1 1 1 1 101) 200 3(H) 41)0 500 600 700 800 900 1000 Window size [sample] 100 200 300 400 500 600 700 800 900 1000 Window size [sample] (a) (b) Figure 9 Arithmetic of M F T and Radix-2 F F T when qMFT - 4 30 The arithmetic of M F T is linear with the number of the computed spectrum coefficients (7VC) and the length of the data record (M). For a given record size (e.g. Figure 8 (a) and Figure 9 (a)) the M F T arithmetic remains the same, with varying shifts, while the F F T arithmetic drops down considerably as the value of shift gets larger. In Figure 10 the floating point operation per D F T is shown for radix-2 F F T , mixed-radix FFT, M F T and the direct D F T , when the consecutive windows are overlapped by N/2 samples (i.e. qMFr = N/2). The arithmetic of the mixed-radix F F T was estimated using the Matlab's fft and flops functions. The radix-2 F F T is very efficient i f the D F T length is power of two, while the M F T is more efficient when TV is a non-composite number (e.g. prime). Floating Point Operations per DFT Window size [samples] Figure 10 Floating Point Operations of D F T algorithms When there is no overlap between the two consecutive DFTs (i.e. qMFT = AO during the spectrum analysis, the M F T has to be applied N consecutive times to obtain the next valid 31 D F T . In this case, the arithmetic of the M F T becomes the same as the direct D F T in Figure 10, while the arithmetic of the other algorithms remains the same. As an example of efficiency of the M F T , consider an TV-point frequency domain filter with a fast convolution method applied to the complex-valued data record. The filter coefficients are pre-calculated and stored and the filter is applied with a radix-2 F F T (or M F T ) , an array multiply, then an inverse F F T (or IMFT) . When radix-2 F F T is used to obtain the convolution real operations are needed. In the case of the M F T the number of real operations needed are: (34) COPSMFT =M[10NC +2] (35) The comparison of efficiencies is shown in Figure 11: 32 Arithmetic of Frequency Domain Convolution Window size [samples] Figure 11 Fast Convolution with M F T and Radix-2 F F T If the M F T is used to compute all the spectrum coefficients, then the radix-2 F F T is more efficient for most of the D F T lengths. But i f only a subset of the spectrum coefficients need to be computed (i.e. sub-band filtering), the M F T / I M F T transformation pair can be more efficient for many values of N. 2.4.3 Advantages and Uses of MFT The computational order of the M F T to recursively calculate the coefficients of an Ap- point D F T is N, a factor of log2N improvement over the F F T . If only a sub-set of the spectrum components are needed, the computing load of the M F T can be further reduced, calculating only the frequency coefficients of interest. The M F T does not rely upon on N being power of two to obtain its efficiency, in contrast to standard F F T algorithms. In this way, the M F T can provide more efficient computation of the D F T when any or all of the following conditions apply: 33 • DFTs are highly overlapped • only a few Fourier coefficients are needed • a specific, non-composite D F T length is needed. Concerning the above properties of the M F T , we can say that it can be useful in different applications of signal processing such as: • on-line computations in real-time spectral analysis • on-line signal identification and detection • speech processing • radar and sonar processing. 34 Chapter 3 Overview of SAR Processing 3.1 Introduction In this chapter the basic geometry of the synthetic aperture radar (SAR) system, the mathematical form of the ideal received signal of a point target and the traditional S A R compression algorithm, the range-Doppler algorithm are introduced. 3.2 Ideal point-target model Assume, the airplane or satellite carrying the S A R antenna travels across the surface of the earth at a constant velocity (V r ) while transmitting microwave pulses at a given pulse repetition frequency (PRF) to the ground with a squint angle © as it is shown in Figure 12. The direction of travel of the S A R antenna is called the azimuth direction while the direction of travel of the transmitted pules is referred to as the range direction ( azimuth and range directions are perpendicular to each other). 35 The transmitted pulses travel at the speed of light (c = 3x10 m/s) which is much faster than the velocity of the antenna. In this way, the antenna can be treated as stationary from the time when it sends out one pulse and receives the reflections from the ground targets. Then the antenna moves to one position to the next azimuth position, sending out another pulse and receiving the back scatter again. Because of the large disparity in time duration of the two directions, azimuth is referred to as the "slow time" axis while range as the "fast time " axis. Figure 12 Synthetic aperture radar geometry 36 Let r\ represent the slow time variable in azimuth direction. Then at each r\ there is a fast time variable t corresponds to the signal in the slant range (R) direction. The transmitted pulse is a chirp signal, exp(y'7TKr t2), and the ideal received S A R signal from a point target can be written as a two dimensional signal [3]: s (t, r\) = P(t) A(r]) exd 2R(V) A (36) In (36) P(t) is the pulse envelope in range direction, A(r\) is the azimuth antenna pattern, Kr is the range F M rate, and X is the radar wave length. The received signal s(t,r]) can be separated to the range signal sft,rj) and the azimuth signal sa(rj) as follows: s(t,rj) = sr(t,ri)xsa(ri) sr(t,rj) = P(t) exp] jnKr sa(T]) = A(77)exp ( t 2R(M v c ,2 \ f 4R(r])^ (37) After the chirp travels through the slant range R(r\) and back, the received signal sr has a time delay 2R(r])/c. In this case, for the same target, but in a different azimuth position, the time delay of the received chirp is different, causing range cell migration in the data memory. 37 From the geometry of Figure 12, the closest slant range of the target (R0) is at azimuth time r\ = no. When the target is at some arbitrary position with respect to the antenna, the slant range R can be expressed as: R(n) = JR?>+V?(n-noy (38) Because Ro» Vfn-rio), equation (37) can be approximated by a parabola, expanding the equation in a second order Taylor series around r|o: R(n) = R0 + 2Rn (39) Combining equation (37) and (39) together, the received azimuth signal can be written as: sa(n) = A(r\) exp = A(r\) exp 77T 4 V exp r 2V2 2^ • i»7t ( l - 1 o ) K0 A ]7l- 4Rn X exp(-jn Ka(n-r}0)2) (40) where the Ka is the azimuth F M rate. Note, that the value of Ro changes for each range cell, therefore the azimuth F M rate changes also, and this must be taken into account when processing the data from different range locations. Also note, that sjr]) has a constant phase -4nRo/X proportional to Ro. This constant phase must be preserved or recovered after the azimuth compression. It is needed for further S A R processing applications, such as S A R interferometry (InSAR). 38 The time variable 77 in the azimuth signal is valid within the exposure period of the target, which is determined by the antenna pattern A(r]). When the antenna length is L , the foot print width of the antenna beam at the target is XRo/L, so the target exposure time is _ A R0 e~LVr (41) Let r\c represent the azimuth time when the beam center crosses the target. Then, r\c = Rotan(©/Vr), and the valid interval of t] for sa(r]) is: T T n c ~ * V < (42) 3.3 SAR signal compression The received S A R data in both range and azimuth can be modeled as the convolution of a linear F M chirp and the ground reflectivity. Using the form of the ideal received signal in equation (36) a matched filter can be derived for each dimension and a pulse compression can be performed on the received data. The pulse compression rearranges the energy received from each ground targets into a single focused pulse. The location of the maximum energy of the pulse corresponds to the location of the target in range and azimuth, while the strength of the pulse represents the reflectivity of the target. The Range/Doppler (RD) algorithm is a traditional, highly accurate and efficient method for compressing S A R data. It consist of the following major stages [3]: 39 • range F F T , • range matched filter multiplication, • range EFFT, • azimuth F F T , • range cell migration correction ( R C M C ) , • azimuth matched filter multiplication and • azimuth EFFT. In the R D algorithm, the range and the azimuth signal are compressed separately using different matched filters. In either case, the matched filtering can be implemented via time domain convolution or frequency domain multiplication. In the rest of this section, the pulse compression is introduced through the frequency domain azimuth compression. The azimuth signal in (40) can be simplified without loss of generality by ignoring the phase term exp(-JTV4RO/A) and the antenna pattern A(n): sa(n) = exp(-jn Ka(t]-rj0)2) The spectrum of the signal in equation (43): Sa(f) = rect dc exp K. 2f% (43) Where Fdc is the Doppler centroid frequency and Fdc = -Ka(nc-r]o). The Doppler bandwidth of the azimuth signal BW=TeKa. The matched filter in the azimuth frequency (Doppler) domain is the complex conjugate ofSa(f): 40 M(D=s:(f\o=0 (44) The frequency domain compressed signal is the product of the spectrum of the matched filter and the spectrum of the azimuth signal: Sc(f) = M(f)-S(f) = red exp(-;'27z:/770) (45) The compressed signal in the time domain is the inverse Fourier transform of Sc(f): Sc(TI) = F{SJf)} = KA Te exP(j2nFdc(r] -Vo)) s i n c ( ^ Te(r)-r),)) (46) In (48), the compressed peak is at the point of the target's closest approach (r\o). This compressed peak location can be changed to other position, such as the target's starting time or Doppler centroid location, by changing the format of the matched filter. In the following two chapters, survey on the applicability of the Momentary Fourier Transform to S A R signal processing algorithms is given. Chapter 4 shows how the M F T can be applied to the S P E C A N S A R processing algorithm, while in Chapter 5 it is shown what advantages the M F T offers when it is used for burst-mode data processing. 41 C h a p t e r 4 A p p l i c a t i o n o f M F T to S P E C A N S A R Process ing A l g o r i t h m 4.1 Introduction As it was mentioned in the previous chapter, S A R signal compression in range and azimuth can be accomplished by cross-correlation in the time domain using time domain convolution, or in the frequency domain, using the fast-convolution variant Range- Doppler method. Alternately, advantage can be made of the linear F M structure of signals by replacing the cross-correlation operation with a frequency filtering operation performed by a D F T . This method is called SPECtral ANalysis ( S P E C A N ) [11]. In this chapter after the theory of the S P E C A N algorithm discussion on the application of M F T to azimuth S P E C A N processing is given. 4.2 The SPECAN Algorithm The S P E C A N algorithm consists two major computational steps: • Deramping 42 • D F T extraction. Deramping is the operation of multiplying a linear F M signal with a complex conjugate reference signal with the same F M rate, but opposite F M slope. The deramping of a signal containing multiple targets is shown in Figure 13. 1 2 3 4 5 6 7 7 8 9 10 11 12 13 14 15 16 17 (a) Frequency - time history of a set equally spaced point reflectors (b) Frequency - time history of the reference function 15 14 13 12 F„ Hz Time F„ Hz Time F. Hz Time (c) Frequency • time history of the product function when the signal in (a) is multiplied by the reference function (b) Figure 13 Deramping of multiple targets 43 After deramping, the targets whose zero Doppler frequency is located within one reference function cycle wi l l have frequencies ranging from -FJ2 to FJ2, where Fa is the azimuth sampling frequency. These targets w i l l constitute one processing region. The formation of parallelogram shaped regions from the deramped targets is shown in Figure 13 (c). Note that each incremental step in the time direction in Figure 13 (a) results an incremental step in the frequency direction in (c), and that the frequency continuity is reset by Fa H z (i.e. / = 0 a n d / = Fa are connected). Consecutive processing regions are separated by lines of slope of the azimuth F M rate (Ka) and constitute a parallelogram shaped area. Figure 14 shows one processing region in more detail. The parameters defined on the figure are calculated as follows: 2V2 • Azimuth F M rate: Kn = r— [Hz/s] AR0 F2 • One cycle of the reference function: M - —3- [samples] • The processed region of the total Doppler spectrum: Mp=(\.-/3)M, where parameter j3 denotes the guard band. • Velocity of a sub-satellite point: Vr [m/s] • Wavelength: A [m] • Closest slant range: Ro [m] • Sampling frequency of the azimuth signal: Fa [Hz] 44 Each deramped target in the parallelogram has a unique frequency ranging from -FJ2 to +FJ2. It is this unique frequency which defines the position of each target with respect to other targets in the same parallelogram. N = input DFT samples •< • Time : G good DFT output samples : discarded DFT output samples Figure 14 Processing regions and the placement of successive D F T blocks in single look case The parts of the energy of any target which originate near the ends of each trajectory (i.e. near the sloped lines) has poor S N R because they arise from the low gain part of the azimuth beam profile and because of the relatively high presence of aliased energy there. 45 The guard band is represented by dashed lines in Figure 14. The total width of the guard band in the time direction is (3 times the distance between the sloped lines. This means that only (i-/3) fraction of the total available Doppler spectrum of each target trajectory is processed. The next step in the processing is to separate the targets into different energy cells with regard to their position in the time domain. This is done by performing short length DFTs across the deramped data. The placement of the DFTs is also shown in Figure 13, where the D F T length is N samples. The location of the first D F T block is arbitrary, but for the sake of illustration it is placed so that the last sample corresponds to the bottom right corner of the processing region. The first D F T rectangle is divided into two parts by a horizontal line where the left-hand side of the rectangle touches the right side of the guard band. The upper section of the rectangle contains invalid output samples which must be discarded because the targets in that region are not fully exposed, while output points corresponding to the lower section of the rectangle are kept as the good D F T output points. If N = ccM, then the unused portion of the available time axis is (l-a-(3) M samples long and the height of the valid part of the rectangle is G = (1-a-p) N (28) D F T output samples. G is the number of good points retained from the D F T operation. 46 The next D F T must be applied so that the frequency of the lowest frequency cell corresponding to the good output region is exactly one cell higher than the highest frequency cell of the current D F T . This placement is shown in Figure 14, as well . The gap between successive DFTs in the input time domain: g = (l-2cc-P)M (29) The D F T length is governed by the desired azimuth resolution pa, by the relation: N=0X9VrFa _ 0 . 8 9 F a A / g 0 (30) where cr is the weighting factor used in the application of D F T such that o~N is the effective number of samples used in the D F T input. From (30) it is seen that the D F T length (AO is directly proportional to the range (Ro) and inversely proportional to the resolution (pa), while the variables Vr, Fa and Ka are defined by the S A R system. The azimuth F M rate strongly depends upon the range while Vr depends weekly upon R for satellite systems and constant for airborne systems. The azimuth resolution can be expressed from (30): _ 0 . 8 9 y r F a _0.89 FaXR0 Pa~ oNKa ~ oN2Vr (31) Table 3 gives available azimuth resolutions for various D F T lengths, with the following C-band S A R satellite parameters: Ka = 2100 Hz/s, Fa = 1650 Hz, Vr = 6800 m/s, a = 0.68. 47 D F T length - N Resolution - p a 256 23 m 512 11 m Table 3 Resolution versus D F T length in S P E C A N for C-band satellite S A R 4.3 M u l t i look processing in S P E C A N The processing scheme of Figure 14 is a single look processing because the azimuth compression operation produces only one output point for each azimuth location. A multi look processing scheme is introduced as follows. The large gaps between DFTs in Figure 14 shows that much of the input data is not used, which is indicative of the excess azimuth bandwidth available when this length of D F T is used. This excess bandwidth can be used to generate multiple looks. Extra looks can be generated by dividing the G good output points from a single D F T into Ni equal groups, and assigning each group to a different look number. Ni is the number of looks. The next D F T is then placed so that the beginning of its good output points are contiguous with the end of the first look of the first D F T . Such a division into looks and placement of the second D F T is illustrated in Figure 15 for the four-look case. In this case the target is fully exposed in four consecutive DFTs , and the fifth D F T would have the same location as the second D F T in Figure 14. 48 Using the good points, obtained from each D F T output (28), the number of good output points per look is L=INT N, (32) Time Figure 15 Divis ion of the good output samples into looks and the location of the D F T operations in multi-look processing Note that G/Ni is not in general an integer so that usually LNi < G and that a small fraction of the good output points are not used. L must be an integer so that the looks fit together evenly at look summation time for each successive D F T . 49 The shift between successive DFTs is: LF R0 A NKa 2 pa Fa{l-P) 0.89 (33) In general q is not an integer, and the nearest integer must be chosen in order to decide where to place the second D F T . For complete efficient data utilization q ~ o~N, so the number of looks that are normally taken can be expressed as follows: I V | =_S_(l.^)-I = A P a _ ( i . / 5 ) - l » ^ a . ' oNKa 0.89V ' a Vr (34) Note that Ni is linearly proportional to pa and does not depend on RQ. 4.4 The SPECAN Algorithm Using the MFT The azimuth F M rate of the received S A R signal is inversely proportional to range, so it changes as the range varies in each cell. In order to keep the resolution and output sampling rate constant across the swath, there is a need to choose different D F T lengths which vary with one or a few samples at a time. The affect of the varying range on the F M rate and the required D F T length for spaceborne and airborne case are shown in Figure 16. 50 Azimuth F M Rate and DFT Length in SPECAN / i ~ / A dmuth FM Rate—» r ' / / - •* / / . j \ : , *- / : >w DFT length i I I r - ' t / ~ / 310 j§ 302 300 830 835 840 845 850 855 860 865 870 Range [Km] (a) Spaceborne S A R processing Azimuth F M Rate and DFT Length in SPECAN ^-DFT lengt • * * • • Azimuth FM Rate^ I o I 300 £ 12 14 Range [Km] (b) Airborne S A R processing Figure 16 Azimuth F M rate and the D F T length with varying range The radar parameters for the two cases are given, in Table 4. Note that in the spaceborne case there is a need for only 16 different D F T sizes to keep the resolution constant through the whole swath, while in the airborne case a wide range of window length is needed. 51 Radar parameter Spaceborne S A R Airborne S A R Units Velocity (V'r) 7000 [m/s] Wavelength (A) 0.057 0.057 [m] Weighting parameter (<T) 0.68 JlllllIMI Guard hand ((I) 0.15 0.15 - Slant range (R) 830 - 870 Lkm| Sampling frequency (F,,) 1680 300 [Hz| Number of looks (A-'/) HBil Is I - Azimuth resolution (pa) 25 4 [m] Table 4 Spaceborne and airborne S A R parameters for S P E C A N arithmetic calculation The Radix-2 F F T algorithm can be only used when the D F T length is a power of two. In other cases of window length a mixed-radix F F T algorithm is used to achieve efficiency only for highly composite N. It means for each different D F T a different F F T algorithm should be implemented within the S P E C A N processor, which makes the D S P architecture rather complex when many lengths of F F T are needed. In contrast to F F T algorithms, the structure and the efficiency of the M F T does not rely upon on the size of the D F T . The same simple algorithm can be used to calculate all of the necessary D F T s during the azimuth compression. The number of real operations of M F T to process an M samples long region, when all of the spectrum coefficients are computed (Nc = AO is given in equation (32). 52 During the D F T extraction only a portion of the spectrum coefficients - good output samples (G) - are used at the same time to obtain the compressed data. Although the amount of these spectrum components remains the same through the processing, but their position changes with the DFTs . So, the conventional form of the reduced-MFT algorithm cannot be used in this case. The sub-band of the calculated spectrum coefficients has to change its positions in the frequency domain in phase with position of the good output samples. The arithmetic of the required reduced-MFT algorithm is introduced below. For the first M F T (DFT 1 in Figure 13) all of the good output samples are computed for the first time so the arithmetic of the first M F T is: MFTOPSDFT, =N [8 G + 2 ] (35) For the next M F T (DFT 2 in Figure 13) the position of the sub-band of the good output points shifted towards to the higher frequencies with L samples (the good output points in a look). Thus, there are L new frequency components to calculate beside the (Ni-L) continuously calculated ones. The number of real operations needed for the L new coefficients of the second D F T are: MFTOPSDFT 2 NEW =N [ 8 L + 2 ] (36) The arithmetic of the previously computed spectrum coefficients of the second D F T is: MFTOPSDFT20LD =q [8(/V, - l ) L + 2 ] 53 (37) The shifting of the interval of the good output points in the frequency domain continues during the whole azimuth processing, which means that the 'new' and the 'old' spectrum coefficients have to be computed previously introduced equations R M - N ^ + 1 times through the whole region. Using the MFTOPS = MFTOPS D F T L + R M - N ^ + 1 (MFTOPSDFT2NEW + MFTOPS D F T 1 0 L D ) = N[SG + 2] + R M - N ^ + 1 N[SL + 2]+ (M-N +1 4 [ 8 ( ^ - 1 ) ^ + 2 ] (38) real operations are needed to process the whole processing region with the reduced M F T . Although equation (38) looks rather complex, the implementation of the above described reduced M F T algorithm is the same as the full M F T algorithm, except the timing and synchronization of the sub-band of the spectrum coefficients. Figure 17 shows the number of operations of the S P E C A N azimuth compression for spaceborne and airborne cases. 54 D F T Operations in S P E C A N D F T Operations in S P E C A N 250 300 350 Window size [samples] (b) Airborne S P E C A N processing Figure 17 Arithmetic of S P E C A N azimuth compression with different D F T algorithms The D F T algorithms used to obtain the D F T output samples are the direct D F T algorithm, the full and reduced M F T , the mixed-radix and radix-2 F F T . Although the radix-2 F F T algorithm is the most efficient, it can be used only once during the process, when the D F T length is 256 samples (Figure 17 (b)). As the window length gets larger ( K a smaller) 55 the size of the processing region also gets larger, thus more operation is needed to process one range cell. The envelope of the plots in Figure 17 (a) and (b) represents the tendency of the growing arithmetic of S P E C A N . Note that the arithmetic of the full and reduced M F T is more uniform in contrast with the arithmetic of the mixed-radix F F T algorithms. Also note, that the arithmetic of the mixed-radix algorithm is equal to the direct D F T i f the window length is a non-composite number (i.e. prime). In Table 5 the ratio between the maximum and minimum of real operation of the reduced-MFT and mixed-radix F F T are given, followed by the ratio between the average operation of reduced and fu l l -MFT and mixed-radix F F T of the full swath. Ratio of flops Spaceborne S A R Airborne S A R NU v o l icd ikvd Ml '1 M i n . of reduced M F T 11.15 Max. of mixed radix F F T M i n . of mixed radix F F T 18.00 141.25 reduced M F T per swath mixed radi \ F F T per swath 1.74 full M F T per swath mixed radix F F T per swath 1.42 2.17 Table 5 Reduced and full M F T versus mixed-radix F F T 56 From Table 5 we can see that in the spaceborne case there is less than 10 % difference between the minimum and maximum of the reduced-MFT flops, while the maximum of the real operations of the mixed-radix F F T more than 18 times bigger than the minimum. In the airborne case, the azimuth F M rate changes more dramatically, thus the difference between the minimum and maximum flops are much larger. When the M F T is applied to S P E C A N there is a factor of 11 difference between the maximum and minimum arithmetic, while in the case of the F F T the maximum of flops is 141 times bigger than the minimum. This large difference in the processing arithmetic makes the timing of the data flow in S P E C A N rather difficult when F F T algorithms are implemented. Note, the results in Table 5 are strongly depend on the azimuth resolution and on the width of the swath. In the S P E C A N algorithm, the resolution is inversely related to the D F T length, thus larger DFTs are needed to obtain finer resolution. As the transformation length is getting wider, the interval of the good output points wi l l shorten, therefore more D F T blocks with higher overlap ratio wi l l be needed to cover the processing region (Figure 14). In this case, the S P E C A N algorithm requires more computation to obtain the azimuth compression. Figure 18 and 19 illustrate the real operations per input samples of the traditional Range-Doppler (RD) algorithm and the S P E C A N algorithm with the M F T and radix-2 F F T algorithms for a spaceborne S A R system. The used system parameters for the analysis are the same as in Table 4, except the number of looks, Ni = 1 in this case. 57 Figure 18 The arithmetic of S P E C A N when the M F T is more efficient Operations/Input samples vs. DFT length in SPECAN and RD Operations/Input samples vs. DFT length in SPECAN and RD Full M F T Reduced M F T Radix-2 FFT Range-Doppler S 0.03 o ! » 0-025; 1 1 1 1 1 1 i :... - vC-SRECAN-wiuvRadiX" 2. FFT- - Range—Do ppicr 801) 900 1000 Window size [samples] (a) [) 900 1000 Window size [samples] (b) Figure 19 The arithmetic of the Range-Doppler algorithm and the S P E C A N algorithm using the M F T and F F T algorithms As it was shown in Section 2.4.2, the M F T algorithm gets more efficient as the Radix-2 F F T when the shift between successive DFTs are not bigger than five samples. Small amount of shift between the DFTs and the corresponding resolution are shown, in Figure 58 18 (a). The shift between DFTs is smaller than six samples only for large window lengths (1181, 1182, 1183, 1184, 1185 samples), close to the maximum (Nmax = (l-fi)M = 1186). Thus the M F T algorithms are more efficient than the Radix-2 F F T algorithm only for very fine resolutions. In Figure 18 (b) the real operations per input samples of the S P E C A N algorithm with the fu l l -MFT, reduced-MFT and radix-2 F F T are shown. Figure 19 (a) illustrates the arithmetic of the R D algorithm and the S P E C A N with different D F T implementations. A s the D T F length gets larger, thus the azimuth resolution gets finer (i.e. N > 700 samples, pa < 10.85 m) the R D algorithm gets more efficient, even if the Radix-2 F F T is applied in S P E C A N (Figure 19 (b)). Thus, when fine resolution is needed, the R D algorithm is used for S A R signal compression. From Figure 18 (b) and 19 (b) it can be seen that the F F T algorithms are more efficient than the full- or reduced-MFT for all transformation lengths usually used in S P E C A N , thus the M F T does not improve significantly the computational efficiency of the S P E C A N algorithm. Beside the complexity and computational efficiency, another important issue in the S P E C A N algorithm to keep the output sampling rate constant. In other words, targets which are T seconds apart in azimuth input time must appear T seconds apart in the output data. It was shown in [11] that the azimuth output sample rate is ^ = ^ - H z (40) 59 The output sampling rate strongly depends on the azimuth F M rate, so when Ka changes trough the swath, there is also a need for slowly varying D F T length to keep the output sampling rate constant. Figure 20 shows Fou, as the function of range, when M F T and radix-2 F F T is applied in the S P E C A N algorithm for spaceborne and airborne systems. Note when the M F T algorithm is applied, the output sampling rate is more uniform for both cases. During the application of the radix-2 F F T only two transformation lengths - 256 and 512 - can be used for the airborne case, while only one, N = 512 can be used for the spaceborne case. This is the reason of the large migration of the output sampling rate when the radix-2 F F T is applied to the S P E C A N . Output Sampling Rate Output Sampling Rate : ! Radixr2 FFT: N = 512 MFT: £ 6 5 a <2 ,s6 0 o. E i 55 Br S 5 0 Radix-2FF'T: N -.256and 512 . 835 8411 845 850 855 860 865 870 Range [Km] 12 14 Range [Km] (a) Spaceborne case (b) Airborne case Figure 20 The output sampling rate of the S P E C A N algorithm In this chapter, the applicability of the M F T to S P E C A N S A R processing algorithm has been investigated. Although, the M F T does not improve the computational efficiency of 60 the S P E C A N algorithm, it has the following advantages over the radix-2 and mixed-radix F F T when they are applied to S P E C A N : • The M F T has consistent computing load as the D F T length changes. • It is easier to implement the M F T algorithm for variable window length. The architecture of a S P E C A N processor using M F T is less complex, because the same M F T algorithm can be used for the different window lengths. • The full radar resolution can be achieved, because the full Doppler spectrum of the targets can be used for the compression by using high-overlapped DFTs . • The output sampling rate of S P E C A N is more uniform. 61 Chapter 5 Application of MFT to Burst-mode SAR Data Processing 5.1 Introduction Burst-mode operation is used in S A R systems, such as R A D A R S A T , to image wide swaths, to save power or to save data link bandwidth. In this operational mode the received data has segmented frequency-time energy in its Doppler history. There are several algorithms for burst-mode data processing: one of them is the Short JLFFT (SIFFT) algorithm which was proposed by Dr. Frank Wong at M D A [17]. In this chapter, after the introduction of the burst-mode data and the SLFFT algorithm properties, a survey on the efficiency of the SLFFT and the applicability of the M F T to the SLFFT algorithm is given. 5.2 Burst-mode SAR processing Burst-mode is commonly used in S A R systems in ScanSAR mode, where the beam is switched between two or more swaths to maximize the imaged swath width. A 2-beam ScanSAR mode is illustrated Figure 21. 62 Figure 21 Burst-mode operation in 2-beam ScanSAR case In this operation mode, the radar beam scans through one sub-swath for a certain time interval, and switches to the next one. After scanning through the second sub-swath, the radar switches back to the first one to start the next burst cycle. The burst cycle has to be short enough to make sure each target is fully exposed in at least one burst. The data from one sub-swath have to be processed separately from other sub-swaths, because the radar beam covers a different ground area in different sub-swaths. In the azimuth direction, the data is segmented into discrete bursts (shaded area in Figure 21), while in range the signal of any sub-swath is continuous, thus the data is not acquired in discrete range bursts. 63 16 fully exposed targets N S . Continuous signal ^ ^ ^ ^ ^ Burst signal o IS 3 a. a. m — *. 5 "> O O) Target 1 Target 2 Target 3 Target 4 Target 5 Target 6 Target 7 Target 8 Target 9 Target 10 Target 11 Target 12 Target 13 Target 14 Target 15 Target 16 Frequency-time diagram of 16 targets slope = KJHz/s] burst 1 burst 2 l _ burst 3 Target 1 Target 16 burst 4 Frequency Doppler history of 16 targets Figure 22 Burst-mode processing of 16 fully exposed and evenly spaced targets in one range cell 5.3 Properties of fully exposed targets in burst-mode processing A typical 2-beam burst-mode data collection pattern is shown in Figure 22. Data from 16 evenly spaced, fully exposed targets in one range cell are shown, where the burst length is 20% of the aperture length. Dashed lines show the azimuth exposure time of each 64 target i f the S A R were operating in continuous mode, and the solid parts of each line show that part of each target actually exposed in burst mode. Note that the part of the target's exposure captured in burst mode varies with each target, which is illustrated in the frequency-time diagram of Figure 22. Each successive target is received at a lower Doppler frequency within a given burst, but is later captured at a higher frequency in the next burst, as long as it stays within the beam. The pulse repetition frequency (PRF or Fa) and the aperture length in azimuth time (Ta) are connected through Ka, Fa=TaKa [Hz] (47) Ta in the time-domain consist of 5 burst lengths, so Fa in the frequency-domain also consists of 5 burst bandwidths (Figure 22). The Doppler history of the 16 targets is also shown in Figure 22, where it is seen that it takes up to 2 bursts to cover all of them. The plot shows the distribution of target spectral energy that would appear i f an azimuth D F T were taken over the full 4 bursts, thus the D F T length was 4 bursts plus 4 gaps long in this case. Note that some targets appear in 2 full bursts (e.g. Target 6), some appear in 3 full bursts (e.g. Target 11), while others appear in two full and one partial burst. In this case, the average number of target exposures or bursts per aperture is 2.5. The number of bursts/aperture in ScanSAR systems is typically between 1.5 and 3. 65 If single-look complex processing is to be done, then there is a choice of which bursts to use for each target. Normally, the target exposures closest to the Doppler centroid (Fdc) wi l l be selected, as shown by the heavier lines in the lower part of Figure 22. However, other bursts may be chosen, when the data is processed for InSAR purposes. r - £ a E < First partially exposed target Frequency-time diagram of fully and partially exposed targets Last partially exposed target a o co 3 a. a. in j? <u O ut First fully exposed target • Last fully exposed target H burst 1 exposed target burst 2 From targets previous to synthetic aperture From targets after synthetic aperture Last partially exposed target First partially exposed target slope = K, [Hz/s] burst 3 First fully exposed target Last fully exposed target burst 4 First partially exposed target burst 1 burst 2 burst 3 Frequency Doppler history of fully and partially exposed targets Figure 23 Burst-mode processing of fully and partially exposed targets in one range cell 66 5.4 Properties of partially exposed targets in burst-mode processing In the previous section the properties of fully exposed targets of a synthetic aperture were introduced when they are processed in burst-mode. These are the valid targets of continuous-mode processing: they have complete frequency-time history and get compressed using the full Doppler aperture (Fa) matched-filter. Beside the fully captured targets, there are also partially exposed targets both at the leading edge and trailing edge of the azimuth D F T . These targets are incomplete, so they are discarded during continuous-mode processing. The frequency-time diagram and Doppler history of the partially exposed targets are shown in Figure 23, where the light gray region corresponds to signals from targets that begin previous to the start of the D F T and the dark gray region corresponds to signals from targets which end after the end of the D F T . Note that most of the partially exposed targets are also completely captured by one or two bursts, so they can be fully compressed using one burst width of their spectrum. In this way more targets with lower resolution can be fully compressed in the same synthetic aperture as in the continuous-mode for a given D F T length. In Figure 23, both partially exposed regions in the Doppler history have a triangle shape in frequency-time space and they complete each other to a rectangular shape. Originally, the position of the light gray region is before the region of the fully exposed targets, corresponding to the position of targets to the synthetic aperture in azimuth time domain. The change in the position of the light gray region is caused by the wrap-around properties of the D F T . 67 First partially Frequency Figure 24 Effect of the circular convolution on the targets Doppler history Figure 24 shows the original parallelogram region of the Doppler history and how the light gray region gets 'wrapped-around', when frequency-domain circular convolution is used instead of time-domain linear convolution. Note, targets in the parallelogram region get compressed into different output cells, while when the top triangle gets wrapped- around different targets can be compressed into the same cell depending on which Doppler sub-band is used. This property of the wrapped-around target Doppler history has to be taken under consideration, when selecting good targets during burst-mode signal compression using DFTs . 68 5.5 T h e S I F F T A l g o r i t h m Most S A R processing algorithms are based on the fast convolution principle where a frequency-domain matched filter is used in the azimuth or Doppler frequency domain. When this method is applied to burst-mode data, the inter-burst gaps are filled with zeros and all the bursts are compressed at once using a full length matched filter followed by an IFFT. However, the compressed targets are then left with a burst-induced modulation. The SIFFT algorithm differs from the conventional fast convolution algorithm in that short, overlapped IFFTs are taken after the matched filter multiply in the Doppler domain [ 2 2 ] . So, when one burst of a target is fully captured by the IFFT, little or no energy from adjacent bursts of the same target is present in the same IFFT. In this way, each IFFT compresses a group of targets without interference (modulation) from other bursts, and an accurate impulse response is obtained. The IFFT is acting like a band-pass filter to extract target energy from the segmented form of the targets' spectra. The filter is time varying in the sense that each successive IFFT is applied to a different frequency band. To capture a target fully, the length of the IFFT must be at least as long as the length of the bandwidth of one burst. The bandwidth of one burst is N K BWHz = [Hz] F a BWbin = N b K " ^ F F T [frequency bins] (48) 69 where Nb is the burst length and NFFT is the length of the azimuth F F T in time samples. Then, the minimum LFFT length is A W „ = BWbin = [samples] (49) The EFFT cannot be longer than the bandwidth of one burst plus one gap, so that a fully- exposed target is not contaminated by a partial exposure of the same target at a different frequency. In the 2-beam ScanSAR case, the length of the gap is often equal to the burst length, so the maximum length EFFT is N =7BW = 2 N b K " N f f t J V / F F T max ^ u y v bin ^ j-,2 (50) In practice, these length limits must be modified a little because of the spreading of target energy in the frequency domain, i.e. a guard band is used when locating the IFFTs. Note, NIFFT max and N IFFT min are proportional to Ka and NFFT- The effect of this property is discussed in Section 5.6, where the arithmetic of the SEFFT is given. Usually, NIFFT max is smaller than NFFT, so less than the whole Doppler spectrum is used for azimuth compression, which means that the output resolution of the SEFFT algorithm is smaller than the maximum available by a factor of NIFFT/NFFT. Locations of NIFFT max and N IFFT min to extract targets with the highest energy using the closest burst spectrum to Fdc (thick lines) are shown in Figure 25. Targets, which are after or before the synthetic aperture and partially exposed, are noted as A 1 - A 2 4 and B1-B24, respectively. 70 F igure 25 H o w m i n i m u m and m a x i m u m L D F T is p laced to compress groups o f targets f rom each burst It can be seen that IFFT Max I captures the complete energy of a single burst spectra of Targets 1, 6-11, 16-A5 B10-B15 and B20-24. F o r these targets, IFFT Max I does not 7 1 extract any energy from other bursts' spectra, so their impulse response is not corrupted by modulation. Similarly, IFFT Max 2 captures the complete energy of a single burst spectra of Targets 1-6, 11-16, A 5 - A 1 0 and B15-20 and between the two IFFTs, all the targets are correctly compressed. In order to form a continuous output image, the results of successive IFFTs are stitched together. If only bursts with the highest energy are used to compress targets, then each output target gets to a different output cell. Note, not all the partially exposed targets can be compressed with the best S N R and Targets A16-24 cannot be compressed at all , because none of the bursts covers them completely. Also note, targets before the full synthetic aperture (B10-B24) get compressed at last and their position in the output array is unique and correct. 5.5.1 Number of good output targets of a single IFFT It was shown in the previous section that each IFFT compress different groups of targets correctly. We refer to these correct results as "good" output targets of one IFFT (GIFFT)- The General form of GIFFT w i l l be derived as below. The number of groups of good output targets of an IFFT in the 2-beam burst-mode processing case is N _  l y FFT N group (51) Note, Ngr0Up is equal to the number of complete bursts in a D F T . In our example in Figure 25, the D F T consists of 4 burst and 4 gaps, so there are 4 groups of good output targets. 72 Also note, that Target 1 and Targets B20-B24 are in the same group or burst ( but Target B20-B24 is wrapped-around by the circular convolution). The shift q between two consecutive targets in the Doppler frequency domain is shown on the top in Figure 25 and K N2 q = — | — [ f r e q u e n c y bins] FA ^ IFFT (52) Note, q is proportional to Ka and N FFT, SO the shift varies with range and the length of the D F T . Note, also that although targets are more than one samples apart in the Doppler domain, they placed one cell apart in the output space. The number of good output targets in a group depends on the length of the IFFT and shift between targets: N - BW n — IFFT ""bin , 1 kgroup ~ " +  1 (53) Note, when NIFFT = NIFFTMIN, then G g r 0 u p = i as it is shown in Figure 25. Then, the number of good output targets from an EFFT can be expressed as: ( T , „ ™ — (j • N IFFT group group N FFT 2Nh (NIFFT BWBIN)FA NIFFT ^ KANFFT (54) 73 5.5.2 Real data simulation of burst-mode processing A real data experiment was done to verify the above described properties of the SIFFT algorithm. A burst-mode processor was created by combining the M D A / U B C dtSAR processor and our Matlab azimuth processing program. In the experiment, the raw data ingestion, range compression and R C M C are done by the dtSAR processor. The R C M C - e d data, along with the required processing parameters, such as Fa, Ka and the Doppler centroid frequency are read into our Matlab program for azimuth processing. Before we can apply the burst-mode algorithm to the d tSAR output data, we have to generate the burst mode signal and correct the antenna pattern to avoid scalloping in the output. The antenna pattern correction is done by summing the azimuth F F T of each range cell, and polyfitting the summed data. The burst mode data is emulated from S A R continuous-mode data by windowing the signals in the azimuth direction in the Matlab program. The parameters of the ERS-1 data are shown in Table 6. The single look detected S A R image generated with the SIFFT algorithm is shown in Figure 26. Processing parameter Value Unit Sampling lrequv.'iK> i / i 1679.90 [Hz] Doppler Centroid (F^ c) 447.01 [Hz] Burst length (N,,) Samples Range cells 1024 Samples Azimuth samples (TV/v /) 2048 Samples Table 6 ERS-1 parameters for real data simulation 74 o O o o o o o o un o i n o un o « CO CO CSI (M i— r- < injiunzv Figure 26 S L C product of burst-mode simulation Although, the synthetic aperture length varies with range, the burst length is approximately l / 5 t h of the aperture through the whole processing region. During the azimuth compression, NIFFT MIN was used for bursts close to the Doppler centroid, so the burst-mode data is extracted with the best SNR. The resolution of the output image is l / 5 t h of the original continuous S L C , because maximum only l / 5 t h of the targets' Doppler spectrum can be used for signal compression. Figure 27 The Doppler history of real burst-mode data The Doppler history of one range cell of the real data was generated using the LMFT algorithm is shown in Figure 27. LMFTs, with size of NIFFT MIN were taken, at each frequency bin, and the result of the transforms were stored in a matrix. Figure 27 shows 76 the values of the matrix, where the brightness of the image corresponds to the magnitude of the output result. The bright points represents the good output answers, while the darker parts show the gaps in the frequency domain. Note, the energy leakage of the good output points in the figure. When the I M F T is applied to process a data record, it needs N samples to fully contain data. So, in the output space there is an ./V sample long ramp before the first valid result. Similarly, after the last valid result there is an N samples long ramp as the I M F T is sliding off from the data. These properties of the frequency analysis cause the leakage in the output space. Ideally Figure 27 would look like the Doppler history in Figure 23, thus it would have the same intervals of gaps and targets. 5.6 Efficiency of SIFFT using the IMFT and the IFFT algorithms To show the efficiency of the rMFT vs. the IFFT algorithm we use the parameters of the alternating polarization (AP) mode of E S A Envisat satellite. Alternating polarization mode provides high resolution products (approximately 30m) in any of the seven swath located over a range of incidence angles spanning from 15 to 45 degrees with polarization changing from sub-aperture to sub-aperture within the synthetic aperture. Effectively, a 2-beam case ScanSAR technique is used but without varying the sub-swath. The velocity of the satellite Vr = 6700 m/s, the radar wavelength X = 0.0567 m and the azimuth F F T is 2048 and 4096 samples long during the efficiency evaluation. Other parameters of the seven swath are given in Table 7. 77 Swathes Sampling frequency F a [Hz] Burst/Gap length in [samples] Range [km] IS 1 •f "';'::!':E"""I:.":E! h !̂h!-L:JP4pNî ^̂ ^̂ ^̂y|̂ l̂i:!i!!i!i:H!::: 825 - 864 IS 2 1645 196 843 - 891 IS 3 2096 887-934 IS 4 1680 218 929 - 990 IS 5 H H 9 N H I <>S3 - 1032 IS 6 1698 238 1027- 1087 IS 7 IOMI Table 7 Envisat swath parameters 5.6.1 Effect of varying SAR parameters and SNR/efficiency tradeoffs As it was shown previously, the azimuth F M rate of the received S A R signal is inversely proportional to range, so it changes as the range varies in each range cell. While the length of the bursts is constant in the azimuth time domain, the bandwidth of the bursts (49) varies with range because it is proportional to Ka. Table 8 shows the maximum and minimum value of the burst bandwidth of each swath in FIz and in frequency bins. A t close range the bandwidth has its maximum, while at far range has its minimum as it is shown in Figure 28 for swath IS1. In a swath, the minimum IFFT length (ND?FT min) should be at least as long as the maximum BWbin to compress all the targets in each range cell correctly. Note, the lower and upper limit length of the IFFTs related to the burst bandwidth so they are different for each swath as well . 78 BW H z [HZ] BWbin [frequency bins] Swathes N F F T = = 2048 N F F T -= 4096 M L N M A X M L N M A X M L N M A X IS 1 212 222 259 271 ^517 542 IS 2 212 221 2o4 279 527 557 • - I S 3 208S/:-': 219 203 214 4Mo 428 IS 4 208 221 253 :~o 506 539 IS 5 204 214 201 211 402 422 IS 6 204 216 246 261 493 521 IS 7 210 198 397 l l o Table 8 Min imum and maximum burst bandwidth of the seven swathes The S N R of the SLFFT algorithm depends on the ratio of the LFFT and the burst bandwidth [22]. The S N R is maximum when the LFFT length is equal to the burst bandwidth (NIFFT = NIFFT min), while the S N R is 3 dB lower when the LFFT window is one burst plus one equal-sized gap long (NIFFT = NIFFT max)- So, in order to keep the S N R loss below a certain value across the swath, there is a need to choose specific azimuth LFFT lengths. The change in S N R for a given LFFT length is: dSNR = 10 log 1 0 J V IFFT KBWbinj (55) If we choose NIFFT as small as possible at near range and have it stay the same through the whole swath, the S N R wi l l decrease, as the burst bandwidth decreases with range. 79 Figure 28 shows how dSNR changes with range, in the case of IS1 swath and NFFT = 2048. The decrease in S N R is zero at near range and rises to about 0.2 dB at far range. Note, although BWbin in (54) depends on NFFT, the gradient and maximum of dSNR is the same for different azimuth F F T lengths. Figure 28 Burst bandwidth and dSNR of IS1 swath 5.6.2 Arithmetic of the SIFFT algorithm using the IMFT and the IFFT algorithms During the efficiency evaluation of the SLFFT algorithm, the ful l -LMFT, the reduced- I M F T and the mixed-radix LFFT algorithms are considered. A formula of the arithmetic of each algorithm is introduced below. It can be shown that the number of LFFTs to compress bursts with the highest energy is 80 NIFFT = 2 N N b G FFT group 2NlFFTNbq NFFT(NlFFT-BWbin+q) (56) We can see from (56) that NIFFT inversely proportional to Gsr0uP, thus i f more targets are in a group less LFFTs are needed for data extraction. Than, the number of operations needed to compress all the targets using the LFFT algorithm is OPIFFT = NIFFT • OPN ,FFT (57) where OPNiFFT is the number of operations needed for one N sample long mixed-radix LFFT. Note i f N is power of 2 than OPN,FFT = 5N log2(N). As we saw in section 1.4.2, in case of the ful l-LMFT algorithm OPIMFT = M(8NIFFT + 2) real operations are needed to analyze an M-point complex data record. In the case of 2- beam burst processing M = 3BWbin, so the arithmetic of the ful l -LMFT algorithm is OPIMFTFULL=3BWBIN(%N IMFT + 2) (58) During the LDFT extraction in the SLFFT algorithm, only a portion of the output target space - 'good' output targets (GIFFT) - is compressed correctly. Although the amount of the 'good' targets remains the same through the processing of a range cell, but their position changes with the position of the LDFTs. So, a simple reduced-LMFT algorithm cannot be used for the target extraction. The position of the computed Doppler - spectrum 81 coefficients has to change in phase with position of the not-compressed good output samples, and the Doppler-frequency coefficients of the already compressed targets do not have to be computed during the rest of the processing. The arithmetic of the required reduced-MFT algorithm is given below. At first, the reduced - I M F T algorithm has to be applied NMFT times to give the first valid compression result. This requires OPIMFTReducedl=N IMFT (SNIMFT+2) (59) real operation. When the first I M F T is done, GIFFT number of targets are compressed correctly, so in the next I M F T NReduced 1 = NMFT - GIFFT number of targets have to be compressed. Than, the I M F T window is shifted q times, a sample at a time, t i l l it fully covers the next target in the group. Now, Ggr0Up number of targets can be extracted correctly, so during the next q shift Nneduced - GSROUP targets need to be extracted. The deduction of the spectrum coefficients from the previously reduced whole output target cell repeats (3BWbin - NiFFrVq times trough the processing region, until all the targets get compressed. It can be shown that the arithmetic of this procedure is: 3BW -N OPIMFTReduced2 =2[3BWbin-NIMFT][4(NlMFT-GIFFT)-2GIFFT ( * ^ L - l ) + l] q (60) Using the previously introduced equations 82 0P1MFT, Reduced = OPIMFT, Reduced 1 + OPIMFT, Reduced 2 = NIMFT(SN, IMFT + 2) + IMFT IMFT - l ) + l] (61) real operations are needed to compress all targets using the reduced-IMFT algorithm. Note that the arithmetic of the reduced-IMFT depends on GIFFT and G&ROUP, thus i f more targets get extracted by an I M F T than less computation is needed. Note, the implementation of the reduced-IMFT is the same as the fu l l - IMFT algorithm, except the timing of the position of the output targets need to be compressed. During the efficiency evaluation, we choose the IDFT lengths on the principles that we want: • maximum S N R at near range, • minimum sampling rate at near range, and • the sampling rate constant with range, and we set the forward azimuth F F T to be 2048 and 4096 samples long. First, we make the IDFT as small as possible at near range (i.e. NIDFT = Max BWbin), and have it stay the same with range, even though the burst bandwidth decreases with increasing range. Thus there is only one IDFT window length to choose from (w = 1) for the IFFT or the I M F T algorithms. 83 Second, we consider the case where the IDFT is allowed to be up to 4 samples longer than the minimum (i.e. Max BWun ^ NIDFT < Max BWbm + 4). This allows some flexibility in choosing a favorable IFFT length from 5 different window sizes, at the expense of a small decrease in S N R . Table 9 shows the IDFT lengths and the maximum of the corresponding dSNR for the seven swathes of the Envisat. Swathes N F F T = 2048 N F F T = 4096 w = 1 w = 5 w = 1 w = 5 NIDFT Max dSNR NIDFT Max dSNR NIDFT M a x dSNR , NIDFT Max dSNR IS 1 271 0.20 275 0.26 542 0.20 546 0.24 IS 2 279 0.24 280 0.26 557 0.24 560 0.26 is ̂  - 214 (»23 216 0.27 42S 0.23 • 432 0.27 IS 4 270 0.28 270 0.28 539 0.27 540 0.28 IS 5 211 0.21 215 0.29 422 0.21 425 0.24 IS 6 261 0.26 264 0.30 521 0.24 525 0.27 208 : : 210 0.26 416 0.20 420 0.24 Table 9 The length of the IDFTs and the corresponding dSNR The burst bandwidth is direct proportional to NFFT (48), so it changes with the same ratio as NFFT does. Note that i f the value of BWbin is bigger, than it is easier to find a highly composite number in its neighborhood, thus it is easier to pick a suitable window length for the IFFT algorithm. Also note in Table 9, that the change in dSNR when a more 84 suitable window is used for the LDFTs is less than 0.1 dB, thus the S N R decrease is practically negligible. Figure 29 shows the number of real operations of the SLFFT algorithm when it is used to compress the data of IS1 swath. The LDFT algorithms used to obtain the azimuth compression are the mixed-radix LFFT and the full- and reduced-LMFT. In Figure 29 (a), we see that the I M F T algorithms are more efficient trough the whole swath when only the minimum window length can be used. The LFFT window begins to cover the bandwidth of more than one target as BWbm decreasing with range, so there is a change in the value of Ggr0up and GIFFT when the LFFT window starts to fully cover two or more targets in the Doppler domain. The I M F T arithmetic decreasing slowly with rage and the arithmetic of the mixed-radix LFFT drops down to its half when Ggr0Up doubles from 1 to 2. The LFFT arithmetic is constant on both side of the down-step. Note that there is more than a factor of 10 difference between the arithmetic of I M F T and LFFT, because the window length is prime (271), so only the direct LDFT can be used to obtain the LFFT results. 85 IS1 Swath Number of IDFT windows to choose from = 1 IS1 Swath Number of IDFT windows to choose from = 5 B25 £.20 JlS IFFT full-IMFT reduced-IMFT 3.5 2.5 IFFT full-IMFT reduced-IMFT 825 830 835 840 845 850 855 860 Range [km] 825 830 835 840 845 850 855 860 Range [km] a) b) Figure 29 Arithmetic of the SLFFT algorithm when applied to the IS1 swath In Figure 29 (b), we can see that the mixed-radix LFFT arithmetic dramatically drops down when the LFFT window length can be a composite number (275). The arithmetic of the I M F T algorithms did not change significantly and they are efficient only in a part of the swath, where only one target can be fully compressed in each group. The down-step of the LFFT arithmetic happens in closer range, because the LFFT window is larger, so it starts to fully cover two targets earlier. Note that the reduced-IMFT arithmetic also drops down when Ggr0Up doubles, but this change is not significant compare to the change in the LFFT arithmetic. Note, the arithmetic of the I M F T and LFFT algorithms of the other swaths is some what similar to the arithmetic of swath LSI. In figure 30, the average millions of operations (MOPS) are shown for all the Envisat swathes when the azimuth F F T 2048 and 4096 samples long. The trend of the arithmetic of the full- and reduced-LMFT is similar in all cases, while the LFFT arithmetic is quite 86 variable, depending upon the composition of the window length. When NFFT = 2048, the I M F T is more efficient for most of the swathes even i f there is an option to chose a suitable window length for the LFFT algorithm. When NFFT = 4096, the I M F T is almost always better than the LFFT i f only the smallest LDFT window can be used (Figure 30 (c)). When there is a possibility to chose a favorable LFFT length, the LFFT is more efficient for all the swathes except one. Note, there is a higher possibility to find a high composite number in the neighborhood of the smallest window size value and more group of good targets (51) can be extracted, when the azimuth F F T larger. So, the arithmetic of the LFFT algorithm decreases more dramatically, than in the case of smaller NFFT length. Number of IDFT windows to choose from = 1, N = 2048 Number of IDFT windows to choose from = 5, N = 2048 Envisat AP mode burst length [samples] Envisat A P mode burst length [samples] a) b) 87 Number of I D F T windows to choose from = 1, N = 4096 Number of I D F T windows to choose from = 5, N = 4096 Envisat AP mode burst length [samples] Envisat AP mode burst length [samples] c) d) Figure 30 Arithmetic of the SIFFT when it is applied to Envisat A P burst mode operation From the above given arithmetic survey we can see that the I M F T algorithm can improve the computational efficiency of the SIFFT algorithm when, • the azimuth F F T is relatively small, • the maximum near range S N R is required and • the IDFT window length is a non-composite number. Beside its efficiency, the I M F T has the following advantages when it is applied to the SIFFT algorithm: • the I M F T has more consistent computing load as the burst bandwidth changes with range and • it is easier to implement the I M F T algorithm for different burst and NFFT lengths, because the same I M F T algorithm can be used for the different I D F T window lengths. 88 Chapter 6 Conclusions 6.1 Summary The objective of this work has been to further develop the theory of the momentary Fourier transform, to examine its arithmetic and efficiency and to show what advantages it offers when it is applied to the S P E C A N S A R algorithm and the SLFFT burst-mode data processing algorithm. The momentary matrix transform was introduced and it was shown when it takes the form of the D F T or the LDFT, the resulting M F T / L M F T have an efficient recursive computational structure. The spectrum coefficients of the M F T / L M F T can be calculated independently and only one complex multiplication and two complex additions are needed to update each spectrum component. This is a factor of log2N improvement over the radix-2 F F T algorithm if all incremental D F T results are needed. The efficiency of the M F T / L M F T do not rely upon the transform length being a power of two, in contrast to 89 standard FFT algorithms. The MFT/IMFT transformation pair can provide more efficient computation of the DFT when: • DFTs are highly overlapped • only a few Fourier coefficients are needed • a specific, non-composite DFT length is needed, and they can be useful in different applications of signal processing such as: • on-line computations in real-time spectral analysis • on-line signal identification and detection • speech processing • radar and sonar processing • narrow-band filtering. After the introduction of the SPECAN SAR processing algorithm, the applicability of the MFT to SPECAN has been investigated. Although, the MFT does not improve the computational efficiency of the SPECAN algorithm, it has the following advantages over the radix-2 and mixed-radix FFT when they are applied to SPECAN: • The MFT has consistent computing load as the DFT length changes. • It is easier to implement the MFT algorithm for variable window length. The architecture of a SPECAN processor using MFT is less complex, because the same MFT algorithm can be used for the different window lengths. • The full radar resolution can be achieved, because the full Doppler spectrum of the targets can be used for the compression by using high-overlapped DFTs. 90 • The output sampling rate of S P E C A N is more uniform. In burst-mode S A R processing, the time-varying spectral properties of the azimuth received data requires that highly-overlapped inverse DFTs used at specific points in the spectral domain to obtain the azimuth compression. It was shown that the I M F T can be more efficient than the LFFT when it is applied to the SLFFT burst-mode data processing algorithm, when: • the azimuth F F T is relatively small, • the maximum near range S N R is required and • the LDFT window length is a non-composite number. 6.2 Future work The research of this thesis project raised the following topics for further study: • Further develop the theory of the momentary matrix transform, and see what other discrete transforms, which are used in signal processing can be efficiently implemented using the M M T . • Investigate i f there is a closed recursive formula of the M M T for shifts greater than one, using higher order permutational matrixes. • Further examine the applicability of M F T to other S A R processing techniques, such as interferometry S A R (InSAR) processing. The M F T can improve the computational efficiency of the local frequency estimation of interferograms, because frequency estimates are needed at every sample point of the interferogram. In addition, there is a possibility of obtaining better resolution estimates of the local frequency using two 91 passes of M F T , in order to provide more accurate estimates of unwrapped interferogram phase. Compare the accuracy, signal-to-noise ratio and computational efficiency of the S P E C A N and the SIFFT algorithm, when they are applied to burst-mode data processing. Investigate how efficiently and with what parameters the SIFFT algorithm could be used to process continuous-mode S A R data. 92 Bibliography [1] A . Papoulis, Signal Analysis. M c G r a w - H i l l , 1977 [2] G . Strang, Linear Algebra and Its Applications. Saunders College Publishing, Third Edition, 1988 [3] J. Curlander and R. McDonough, Synthetic Aperture Radar: System and Signal Processing. Wiley, New York, 1991. [4] J. G . Proakis and D . G . Manolakis, Digital Signal Processing. Prentice Ha l l , Third Edition, 1996 [5] R. R. Bitmead and B . D . O. Anderson, "Adaptive frequency sampling filters," IEEE Trans. On Circuits and Systems, vol. CAS-28 , pp. 524-534, June 1981. [6] J. Dudas, The Momentary Fourier Transform. Ph.D. thesis, Technical University of Budapest, 1986. [7] H . L i l l y , "Efficient DFT-based model reduction for continuous systems", IEEE Trans, on Automatic Control, vol.36, pp. 1188-1193, Oct. 1991. [8] B . G . Sherlock and D . M . Monro, "Moving discrete Fourier transform", IEE Proceedings-F, vol. 139, pp. 279-282, Aug . 1992. [9] S. Albrecht, I. Cumming and J. Dudas, "The momentary Fourier transformation derived from recursive matrix transformations", in Proceedings of the 13th International Conference on Digital Signal Processing, DSP'97, (Santorini, Greece) vol. 1, pp. 337-340, July, 1997. [10] I. Cumming and J. R. Bennett, "Digital processing of S E A S A T S A R data", IEEE 1979 International Conference on Acoustic, Speech and Signal Processing, (Washington, D .C . , U S A ) , Apr i l 2-4, 1979. 93 [11] I. Cumming and J. L i m , "The design of a digital breadboard processor for E S A remote sensing satellite synthetic aperture radar", technical report, MacDonald Dettwiler, Richmond, BC, Canada, July 1981. Final report for E S A Contract No. 3998/79/NL7HP(SC). [12] M . Sack, M . Ito, I. Cumming, "Application of Efficient Linear F M Matched Filtering Algorithms to S A R Processing", IEEE Proc-F, vol.132, no. 1, pp. 45-57, 1985. [13] T. Ngo and C. M . Vigneron, "Project Report: U B C S Q L P v. 1.7.", technical report Radar Remote Sensing Group, Electrical and Computer Engineering, University of British Columbia, 1995. [14] H . Hobooti, "Radiometric Correction in R a n g e - S P E C A N S A R Processing", M.A.Sc. thesis, Electrical and Computer Engineering, University of British Columbia, Apr i l 1995. [15] Sandor Albrecht, "The Momentary Fourier Transform and Its Application to the S P E C A N S A R Processing Algorithm", technical report SA-97-01, Radar Remote Sensing Group, Electrical and Computer Engineering, University of British Columbia, 1997. [16] S. Albrecht and I. Cumming, " Application of Momentary Fourier Transform to the S P E C A N S A R Processing Algorithm", in Proceedings of the IX. European Signal Processing Conference, EUSIPCO'98, (Rhodes, Greece) September 8-11, 1998. [17] F. Wong, "Processing Envisat A P mode with Range Doppler algorithm", technical report, MacDonald Dettwiler, Richmond, BC, Canada, January 1996. [18] I. Cumming Y . Guo and F. Wong, "Analysis and Precision Processing of Radarsat ScanSAR Data", In Geomatics in the Era of Radarsat, GER'97, (Ottawa, Canada), M a y 25-30, 1997. [19] F. Wong, D . Stevens and I. Cumming, "Phase-Preserving Processing of ScanSAR Data with Modified Range Doppler Algorithm", in Proceedings of the International 94 Geoscience and Remote Sensing Symposium, IGARSS'97, (Singapure), pp.725-727, August 3-8, 1997. [20] I. Cumming Y . Guo and F. Wong, " A Comparison of Phase-Preserving Algorithms for Burst-mode ScanSAR Data Processing", in Proceedings of the International Geoscience and Remote Sensing Symposium, IGARSS'97, (Singapure), pp.731-733, August 3-8, 1997. [21] Y . Guo, "Precision Processing of Burst-Mode ScanSAR Data", M.A.Sc. thesis, Electrical and Computer Engineering, University of British Columbia, October 1997. [22] I. Cumming Y . Guo and F. Wong, "Modifying the R D Algorithm for Burst-mode S A R Processing", in Proceedings of the European Conference on Synthetic Aperture Radar, EUSAR'98, (Friedrichshafen, Germany), pp.477-480, M a y 25-27, 1998. [23] S. Albrecht and I. Cumming, " The Application of Momentary Fourier Transform to SEFFT S A R Processing ", in Proceedings of the IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis, (Pittsburgh, U S A ) October 7-9, 1998. 95

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
Japan 5 0
France 5 0
China 4 32
Canada 3 0
United States 2 0
City Views Downloads
Unknown 6 0
Tokyo 5 0
Beijing 4 1
Waterloo 3 0
Ashburn 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items