Advanced Intermittent Clutter Filtering for Radar Wind Profiler: Signal Separation through a Gabor Frame Expansion and its

A new signal processing method is presented for the suppression of intermittent clutter echoes in radar wind profilers. This clutter type is a significant problem during the seasonal bird migration and often results in large discrepancies between profiler wind measurements and independent reference data. The technique presented makes use of a discrete Gabor frame expansion of the coherently averaged time series data in combination with a statistical filtering approach to exploit the different signal characteristics between signal and clutter. The rationale of this algorithm is outlined and the mathematical methods used are presented in due detail. A first test using data obtained with an operational 482 MHz wind profiler indicates that the method outperforms the previously used clutter suppression algorithm.


Introduction
Radar wind profilers (RWP) were developed from MST-Radars (Van Zandt, 2000) and have meanwhile become standard instruments for measuring wind velocities in the atmosphere. Overviews of the technical and scientific aspects of RWP including its signal processing have been provided, among others, by Gage (1990); Röttger and Larsen (1990); Doviak and Zrnic (1993) and Muschinski (2004). Especially the routine application by weather services and the assimilation of the data in Numerical Weather Prediction Models is an indicator for the degree of maturation that this technology has achieved, see e.g. Monna and Chadwick (1998); Bouttier (2001); Benjamin et al. (2004b); St-James and Laroche (2005); Ishihara et al. (2006). However, it is a matter of fact Correspondence to: V. Lehmann (volker.lehmann@dwd.de) that sometimes large and unacceptable differences are observed between the profiler data and independent reference measurements. In many cases these differences are clearly attributable to either clutter echoes or Radio Frequency interference. Spurious signals are often easily discernible in the Doppler spectrum by human experts, but not always adequately handled by the automatic processing. For that reason, research on improvements in wind profiler signal processing has remained a very active field over the last decade.
In this paper we deal with so-called intermittent clutter and propose a new filtering algorithm for the detection and suppression of these clutter signals in the profiler raw data. Of particular importance are echoes caused by migrating birds in spring and fall. It is well known that birds are effective targets for a wide range of radars from X-band to UHF (Vaughn, 1985;Bruderer, 1997a). In fact, most of the knowledge about migrating birds come from radar observations. That concerns in particular their flight behavior under the influence of environmental factors (Bruderer, 1997b). Radar ornithology is meanwhile a mature field and it is no surprise, that birds are also detected by the sensitive radar systems used for wind profiling. The susceptibility of wind profiler radar systems to bird echoes depends primarily on wavelength and antenna characteristics. It mostly affects L-band and UHF-systems, that is Boundary Layer profilers and Tropospheric profilers, as discussed in Wilczak et al. (1995). Intermittent clutter is an issue for the standard Doppler-beam swinging radars as well as for spaced antenna and imaging radar systems, where new mitigation techniques like adaptive beamforming have recently been proposed (Cheong et al., 2006;Chen et al., 2007). We mention in passing that other remote sensing instruments used in Meteorology are also affected by migrating birds (Mastrantonio et al., 1999;Zhang et al., 2005;. Intermittent clutter echoes caused by aircraft were already mentioned by Hogg et al. (1983), and a few years later it Published by Copernicus Publications on behalf of the European Geosciences Union. became obvious that especially echoes from migrating birds can be a serious issue in wind profiling (Ecklund et al., 1990;Barth et al., 1994). If present, such spurious signals can cause a significant deterioration of the quality of the derived winds. To give an example, the investigation of low-level jets using RWP data is hampered by bird migration clutter (Stensrud, 1996). This makes it necessary to either use extensive quality control procedures to identify and skip contaminated data (Daniel et al., 1999;Song et al., 2005) or to limit the studies to periods where bird migration is negligible (Anderson and Arritt, 2001). Many other investigations using RWP data have mentioned the bird contamination problem, e.g. Ralph et al. (1998); Locatelli et al. (1998); Parker and Johnson (2000); Lundquist (2003); Nielsen-Gammon et al. (2007). While the need for an extensive manual data quality control and cleaning might be acceptable for research activities, it is surely not feasible in any operational setting. Nevertheless it is mandatory to avoid the assimilation of bird contaminated profiler wind data, as this can have significant effects on the quality of the forecasts (Semple, 2005). Due to the nature of the problem, a bird migration check at the operational center itself is not the best approach (Benjamin et al., 2004a). While current state-of-the art profilers nowadays run more or less sophisticated algorithms on site to reduce bird contamination (Merritt, 1995;Jordan et al., 1997;Ishihara et al., 2006), practical experience supports the statement that the problem has not been fully resolved.
The problem of bird contamination has been well-known for more than a decade (Wilczak et al., 1995;Engelbart et al., 1998) and it still is a research topic in RWP signal processing. The first successful attempt to reduce bird contamination was made by Merritt (1995), who suggested a selective averaging method of the individual Doppler spectra based on a statistical criterion. The same method can also be applied off-line to averaged spectra, when data with higher resolution are not available (Pekour and Coulter, 1999). Weber (2005) used neural networks for a classification of contaminated single spectra, followed by a selective averaging. Other proposals have concentrated on modified peak detection in the Doppler spectrum to address spurious flier returns, among other clutter types (Griesser and Richner, 1998;Cornman et al., 1998;Morse et al., 2002;Weber et al., 2004). The disadvantage of all these methods is that the mitigation processing builds upon the Doppler spectra (either before or after spectral integration). Given the highly non-stationary characteristics of the intermittent clutter signal, it is necessary to deal with the problem before the Doppler spectrum is estimated, because Fourier methods are generally inadequate for nonstationary signals. In other words, the necessary nonlinear filtering has to be performed in the time domain. This approach was first suggested by Jordan et al. (1997) and further by Lehmann and Teschke (2001), who suggested wavelet decomposition and wavelet coefficient thresholding, to remove the clutter part of the signal. However, the a-priori unclear choice of the mother wavelet and -at least for the dyadic wavelet trans-form -a suboptimal signal separation in the wavelet domain, especially near zero Doppler shift, makes an efficient separation of clutter and signal difficult.
Ideally one would like to have an intermittent clutter suppression algorithm that reduces the clutter part of the signal as best as possible, given the sampled data and that further quantifies its degree of contamination by providing some measure of clutter energy for quality control purposes. Furthermore, the algorithm must not degrade both data quality and availability in the no-clutter case, but it should perform as well as the proven standard processing methods. This requirement is more stringent than it may appear at first glance. In this paper, we propose a new signal-clutter separation method that attempts to meet these objectives. It is based on a redundant frame decomposition of the time series followed by the statistical filtering approach suggested by Merritt (1995).
The paper is organized as follows: Sect. 2 gives an overview of RWP signal characteristics and signal processing and identifies shortcomings of the currently used methods when intermittent clutter signals are present. Section 3 reviews basic results of the mathematical theory of frames, which deals with linear discrete signal representations. The goal is here to find a signal representation, that achieves optimal separation between the atmospheric and the clutter part of the signal. This is achieved by the discrete Gabor representation, which is discussed next. Section 4 focuses on a statistical approach to objectively identify the atmospheric signal component, based on well-justified statistical assumptions. A comparison of the new algorithm with the previously used signal processing techniques is shown in Sect. 5. The data used were obtained during routine operation of a 482 MHz wind profiler radar of the Deutscher Wetterdienst at Bayreuth, Germany in the fall of 2005. Finally, a summary and conclusions are given in Sect. 6.

General properties of the received signal
The relationship between the signal received by the radar and the scattering medium is the topic of radar instrument theory, which basically describes how atmospheric properties are mapped to the measurable function at the radar receiver output (Woodman, 1991;Muschinski, 2004). It is known that models for the scattering processes and the technical properties of the radar system must be considered here, which makes the task quite formidable. However, for the problem at hand it is not required to consider such theories in detail, because we are only interested in some rather general properties of the received signal, like statistical stationarity. For a pulsed RWP, the received signal at the antenna output has the following properties: Ann. Geophys., 26, 759-783, 2008 www.ann-geophys.net/26/759/2008/ 1. Continuous real-valued random voltage signal: Every measurable physical quantity is real. The randomness is the result of the random nature of the scattering process.
2. Intrinsically nonstationary: This is due to the impulsive character of the transmitted signal and the inhomogeneous vertical structure of the atmosphere.
3. Multi-component: Beside the ubiquitous noise, there may be signal contributions from several independent scattering processes, like Bragg scattering at fluctuations of the refractive index, Rayleigh scattering at precipitation and scattering at various clutter targets.
4. Narrowband: The signal is band-limited, with a maximum width that is largely determined by the bandwidth of the transmitted pulse.
5. Large dynamic range: The signal varies easily over many orders of magnitude, which is typical for all radar systems.
After a linear low-noise amplification, the first processing step is a (digital) quadrature demodulation of the analog band-limited signal. This leads to a complex baseband representation, where the signal is described through the time series of its in-phase (I) and quadrature-phase (Q) components. Property 1 is thus modified, because the signal has now become complex. Furthermore, uniform sampling for N fixed delay times (after pulse transmission, corresponding to N fixed ranges) at multiples of the radar inter-pulse period is then applied to generate N quasi-stationary sequences from the nonstationary signal. This stationarity assumption is usually valid for atmospheric scattering, ground clutter and noise, provided the scattering medium at a fixed height does not change its properties significantly over the length of the time series (Woodman, 1991). It is one of the basic assumptions of signal processing for atmospheric radars (Keeler and Passarelli, 1990). The process of generating the N sequences is called range-gate sampling and thereby, property 2 is modified. The remaining signal properties 3-5 are preserved for the N discrete data sequences, provided processing is linear. Finally, matched filtering of the band-limited signal is performed to achieve an optimal signal-to-noise ratio.

Classical signal model and its limitations
The classical RWP signal model assumption is that the demodulated discrete voltage sequence at the receiver output can be written as where I[k]∼N (0, σ 2 I ) and N[k]∼N(0, σ 2 N ) are independent complex zero-mean Gaussian random vectors describing the atmospheric signal and the receiver noise, respectively (Zrnić, 1979), t is the sampling interval of the sequence and ω the mean Doppler frequency. Furthermore I[k] is narrowband compared to the receiver bandwidth and |ω|≤π/ t (Nyquist criterion). Because S[k] is the result of the demodulation of a real valued zero-mean and stationary Gaussian random process, the resulting Gaussian complex random process is also wide-sense stationary and zero-mean. Furthermore, the sequence has a vanishing pseudo-covariance, that is we have E(S[k]S[l])=0. Such a process is usually called proper, circular or phase-invariant (Neeser and Massey, 1993). We will use this property later in connection with a moments theorem for these processes (Reed, 1962).
Because S[k] is Gaussian, it is completely characterized through its covariance matrix R with entries where is specified below. Furthermore, stationarity is assumed over typical dwell-times of O(1 min). While this is a classical assumption in radar signal processing (Zrnić, 1975(Zrnić, , 1979Woodman, 1985;Frehlich and Yadlowsky, 1994;Lottman and Frehlich, 1997), it is unknown for which maximal time series length this assumption can be made safely. We found that bird clutter signals are significantly nonstationary over typically used dwell times of about 30 s to 60 s. This is in sharp contrast to observed atmospheric signals, which exhibit a high degree of stationarity on that time scale, well in line with the classical assumptions. Therefore we get the following expression for the autocovariance function where we set (The sequence ρ will be of importance when constructing adequate mean and variance estimators.) Finally, the autocorrelation function [k] is often assumed to follow a Gaussian correlation model, which corresponds to a Gaussian signal peak in the power spectrum. If the spectral width of the signal is w, then we have (Zrnić, 1979;Frehlich and Yadlowsky, 1994) [k] = e −2π 2 w 2 k 2 t 2 . (3) Note that this Gaussian correlation model must not be confused with the characterization of the random process as Gaussian, which covers a much wider class of signals. The assertions are normally very well justified and therefore often used in simulations of the radar signal (Zrnić, 1975;Frehlich and Yadlowsky, 1994;Muschinski et al., 1999).
www.ann-geophys.net/26/759/2008/ Ann. Geophys., 26, 759-783, 2008 In reality, however, there is sometimes a third component contributing to the signal, namely clutter (Muschinski et al., 2005), so that the signal model must be written as: Clutter is the totality of undesired echoes and interfering signals, therefore it is impossible to generalize the properties of C[k]. In the case of RWP, clutter includes in particular echoes from airborne objects such as aircraft and birds as well as returns from the ground. Interfering signals may be caused by other radio transmitters operating in the RWP receiver band. In the remainder of the paper, we restrict ourselves to intermittent clutter signals.
While the properties of the intermittent clutter component have not been systematically investigated, it is instructive to take a look at a few examples. Such have been presented by various authors: Wilczak et al. (1995) described the distinct characteristic of bird contaminated I and Q data when seen in an A-scope display, but the shown time series taken with a 924 MHz RWP is only 0.5 s long, which is too short to see its essential characteristics. Jordan et al. (1997) show an example of a 30 s long time series taken with a 915 MHz RWP during bird migration, which exhibits a variation in the envelope of the signal due to modulation of signal amplitude by the antenna beam pattern. Another example of intermittent clutter caused by airplanes and a simple theoretical model is given by Boisse et al. (1999). The most distinct feature here is also the time-dependent amplitude of the signal. A 19 s time series of a 482 MHz RWP containing an airplane echo is discussed in Muschinski et al. (2005).
In the fall of 2005, time series data of the coherently integrated I/Q signal of the RWP at Bayreuth, Germany were saved in the wind low mode to get a unique dataset for the investigation of bird migration. For 13 October, it was subjectively judged that the data showed a maximum of bird echoes. We have therefore selected this day for demonstration of the proposed algorithm. One particular dwell is shown in Fig. 3. The time series has a length of about 35 s and its nonstationarity is striking.
When data containing intermittent clutter components are compared to both clear air and ground clutter signals (see Muschinski et al., 2005, for an example), it is very obvious, that the main difference is the transient character of the intermittent clutter signal component. Following Friedlander and Porat (1989), we define a transient signal as a signal whose duration is short to the observation interval, in our case the dwell time. Such a behavior clearly reflects a nonstationarity of the underlying scattering process. It is not the sinusoidal signature that makes the difference, as a sufficiently strong clear air signal also exhibits a sinusoidal nature (see Figs. 1 and 2 in Muschinski et al., 2005) -the most distinct property of intermittent clutter is its nonstationarity.

Consequences for signal processing
Signal processing can be regarded as the art of extracting the maximum amount of information from a given measurement. This obviously means that the general properties of the signal determine the optimal mathematical processing methods. A stationary Gaussian stochastic process is without loss of information described by its time-independent second-order properties, that is the autocovariance function or, equivalently, the power spectrum. This assumption holds when Eq. (1) is valid, and the classical way to process RWP data is then based on a non-parametric estimation of the power spectrum using a discrete Fourier transform of the (usually coherently integrated) raw signal over the dwell-time. The power spectrum is commonly called the Doppler spectrum. Its first three moments are estimated after the noise contribution to the spectrum has been subtracted, to describe the basic properties of the atmospheric signal (Woodman, 1985). However, we have seen that the clutter contribution can be highly nonstationary. If the signal S[k] contains nonstationary components, then the Doppler spectrum is no longer an adequate representation of the stochastic process because information regarding time dependency is already lost. So it cannot be expected that a successful intermittent clutter filtering strategy can be developed based on the Doppler spectrum. Therefore it is tempting to try methods that were developed in the framework of nonstationary signal processing. A necessary condition is obviously a separation of C[k] from the stationary components I[k]e iωk t +N[k]. To achieve this, we look for a representation of the signal in which we are able to discriminate between stationary and nonstationary signal components. This is the goal put forward in Wilczak et al. (1995): Clearly, a superior technique would be one in which the bird signal and atmospheric signal could be differentiated from each other and processed independently.
So far we have considered either a pure time representation of the signal, namely its discrete time series, or its complex Fourier transform as a pure frequency representation. Both are not optimal for transient phenomena, although they are complete representations of the same information. Therefore we look for an intermediate representation that aims at the joint time-frequency structure of the signal, so it needs to depend on both time and frequency. This is the topic of the next section. If we are able to separate stationary and nonstationary signal components in such a representation, then we might be able to suppress the nonstationary clutter part while leaving the stationary signal component essentially intact.

The windowed Fourier transform and the timefrequency plane
Let us consider continuous signals first, although in practice we are always given a discretized signal. A quite natural way to analyze a continuous signal simultaneously in time and frequency is provided by the windowed Fourier transform (WFT), see Gabor (1946); Daubechies (1992); Kaiser (1994);Mallat (1999). It is essentially an extension of the well-known Fourier transform, where time localization is achieved by a pre-windowing of the signal with a normalized window function h∈L 2 (R). For any given function S∈L 2 (R), the WFT is defined as The operator V h maps isometrically between L 2 (R) and L 2 (R 2 ), that is a one-dimensional function/signal is with no loss of energy transformed via the WFT into a twodimensional function depending on both time τ and frequency ω. The (τ, ω)-plane is called the time-frequency (TF) plane or briefly the phase space. This representation was suggested by Gabor (1946) to illustrate that both time and frequency are legitimate references for describing a signal. The squared modulus of V h S is called the spectrogram, denoted by and provides a measure for the energy of the signal in the time-frequency neighborhood of the point (τ, ω) and thus insight about the time-frequency structure of S. However, due to Heisenberg's uncertainty relation, there is no arbitrary resolution in time and frequency simultaneously, i.e. a pointwise frequency description in time domain and a point-wise time description in frequency domain is impossible. Formally, one considers in the uncertainty context for some centralized signal h with h =1, time and frequency variances (7) for which the Heisenberg uncertainty relation yields It can be shown, that equality in Eq. (8) is achieved when h is a translated, modulated or scaled version of the Gaussian function (equality means achieving optimal resolution in the time-frequency plane). Their time-frequency spread is visualized through a rectangle with widths σ t and σ ω in the TF-plane, this is called a Heisenberg box -see Fig. 1. This optimality result shall be used later on, when elaborating a 18 V. Lehmann and G. Teschke: Advanced   discrete version of Eq. (5). Since the WFT is an isometry, the inversion of V h can be performed by its adjoint, Hence, in the continuous setting we still have signal analysis, transform Eq. (5), and signal synthesis, transform Eq. (9), in some straightforward way available and therefore timefrequency signal filtering can be performed in three simple steps (see e.g. Hlawatsch and Boudreaux-Bartels, 1992): 1. Analysis: Computation of the WFT using Eq. (5).

From windowed Fourier transform to Gabor frame expansions
For discrete signals, continuous transforms (5) and (9) are not suitable and would create very redundant representations of the signal. A first adjustment can be achieved when Eqs. (5) and (9) are approximated by discrete sums. Discretizing Eq. (9) means taking only values of the WFT at some discrete lattice in phase space. As it was pointed out, e.g. in Daubechies (1992), the sampling density in phase space www.ann-geophys.net/26/759/2008/ Ann. Geophys., 26, 759-783, 2008 plays a significant role for the existence and stability of a reconstruction formula, i.e. of a discrete version of Eq. (9). Assume we are given some discrete subset (to be specified below) of the TF-plane, then a naive discrete version of the inversion formula (9) would be where the parameter T controls the discrete linear shift mT along the time axis and the sampling shift k in the frequency domain. In order to verify whether Eq. (10) indeed exhibits a reconstruction formula, we first observe that for a family of elementary signals or so-called atoms {h m,k } (m,k)∈ that is complete in L 2 (R) any S∈L 2 (R) can be represented by a linear expansion of the form But only in very specific cases, e.g. when {h m,k } (m,k)∈ forms an orthonormal basis, and then Eq. (10) would indeed be an equality, In general, this is not the case, i.e. we only have where the operator F * F and its properties are briefly discussed in Appendix A. For a detailed analysis and discussion on this subject we refer the interested reader to, e.g., Daubechies (1992). To reconstruct S (i.e. to invert F * F ), special properties on and on the analyzing atoms (the dual functions to h) are required. In what follows, we shall focus on the practically relevant biorthogonal case, in which the construction of the analyzing atoms becomes simple and, moreover, numerically stable. To this end, suppose there is some auxiliary family g m,k (t)=g(t−mT )e ik t (yet unknown) available that serves as a reservoir of analyzing atoms used to compute the Gabor coefficients a m,k via Eq. (5), This approach was originally proposed by Bastiaans (1980). Inserting now Eq. (12) into Eq. (11) yields Equality in the latter equation is assured as long as Condition (13) is called the biorthogonality relation and restricts the choice of g in dependence on the preassigned function h. The particular choice of the window function h (e.g. its variance σ h ), the time shift T and the frequency shift directly controls the existence, uniqueness, convergence properties and the numerical stability of the Gabor expansion (11), which exists for arbitrary signals S(t) only if T ≤2π; this is a frame theoretical result, see Daubechies (1990); Mallat (1999). The physical meaning of this inequality is nothing but the Nyquist sampling criterion and represents the sampling density. T =2π is called critical sampling. This was Gabor's original suggestion, as he was aiming at elementary signals conveying exactly one datum or one "quantum of information". In other words, there was no interest in any redundancy. Gabor (1946) called the sampling density an information diagram. In his attempt to derive a theory of communication, each area represents one elementary quantum of information which Gabor proposed to call a logon. Although conceptually simple and appealing, the Gabor expansion at minimal sampling density in the TF-plane ( T =2π) has no nice mathematical structure. In particular, it does not form a basis with the basis functions localized in time and frequency. A relaxation of the equality T =2π is therefore required and generates a crucial degree of freedom in the Gabor expansion, this at the expense of oversampling and a possible non-uniqueness. For T >2π the stability of the expansion is lost.

Gabor frame expansions for discretely sampled signals
So far we have discretized Eq. (9) resulting in the Gabor frame expansion (11) for S∈L 2 (R). But when it comes to real applications, only finitely many discretely sampled values of S are available; namely S[n]=S(n t), n=0, . . . , N −1. Therefore it becomes necessary to develop a fully discrete concept for evaluating the Gabor coefficients (12). Moreover, the discrete subset in Eq. (11) is in general infinite and hence also not suitable for a numerical implementation. The sum needs to be appropriately truncated and, in addition, a discrete version of the dual function g needs to be derived.
We now illustrate how to proceed for discrete data S. More details can be found in the original paper by Wexler and Raz (1990) and Appendix B. Assume we are given some discrete and finite time (periodic) signalS with sampling points n=0, . . . , N −1, that isS[n]=S [n+N]. We therefore have to periodize the analysis and synthesis windows as well, Slightly abusing the notation, we omit the tilde denoting periodic (finite) functions in the following. The signal S can be discretely represented by whereas the Gabor coefficients can be derived from Introducing integers M and K and the toral component , the discrete analysis and synthesis windows can be rewritten as As can be seen, M denotes the time and K the frequency step size. They correspond to T and . In our setting they are constrained by M·M= K·K=N. The reconstruction formula becomes where we have assumed that the following discrete version of biorthogonality relation (13) for the sequences h and g is fulfilled, It can be shown (for a proof see Appendix B) that the biorthogonality relation is satisfied if for 0≤p≤ M−1 and 0≤q≤ K−1. System (16)  the vector representing the discretely sampled dual frame, and let A be the matrix of size M K×N with entries A (p,q),j =h(j +qK)W jpM N , then the dual frame atom g is the solution of the linear system For oversampling M K<N, system (17) is underdetermined, and the solution is no longer unique and therefore there is a variety of possible dual frame atoms g.
3.4 On the choice of the analysis and synthesis atom and the TF-plane lattice As we have seen, there is a high degree of freedom when constructing a frame representation of some signal S. In particular, i) the choice of the synthesis window h ii) the choice of the time-frequency sampling grid , i.e. the choice of M and K, which specifies the redundancy/non-redundancy and therewith the nonuniqueness/uniqueness of the Gabor frame expansion (14) iii) the choice of g in case of M K<N, i.e. in the oversampling situation, one may add further desirable constraints on the solution g of system (17), e.g. minimum energy-norm.
These three aspects shall now be discussed: At i): Any absolute and square integrable function h is appropriate. However, as mentioned above, Heisenberg's uncertainty relation (8) requires for optimal time-frequency resolution a Gaussian function. Therefore, we choose where the scaling parameter σ h (determined below) shall allow either a better resolution in time or in frequency. As we shall see in iii), the time-frequency localization properties of synthesis function h carry over to analysis function g.
At ii): The most important parameters that control the sampling density in the TF-plane are K and M. Together with the specification σ h they fully determine (up to noncanonical choices of g) the discrete Gabor representation of some given function. In principle, the only requirement is K M≤N.   are due to Nyquist's law automatically spaced with resolution 1/T d within [−1/2 t, 1/2 t]. Through the flexibility of the Gabor representation we may individually set up the time and frequency spacing. Let us consider to this end the Heisenberg box size, i.e. the time and frequency variances (Eq. 7), which take for our particular h the form σ 2 t =σ 2 h /2 and σ 2 ω =(2σ 2 h ) −1 . If we restrict the spacing of the TF-plane to this box size (essentially smaller would produce an overlapping of the boxes), i.e. setting τ = M t=σ 2 t and ω= K/T d =σ 2 ω , Heisenberg's uncertainty principle (Eq. 8) and the solvability of Eq. (17) yields The right inequality in Eq. (19) represents an upper sampling bound, which prevents an unnecessary Heisenberg box overlapping. If now an application requires a time resolution τ in the Gabor representation, we immediately obtain in the context of Heisenberg's uncertainty principle the optimal scaling factor for the synthesis (and therewith for the analysis) atom, and a suggestion for the sampling density in time and frequency, At iii): In the oversampling situation ( M K<N), the nonuniqueness can be used to add desirable constraints to the solution, for example minimum energy. This was discussed in greater detail in Qian and Chen (1993) and Qian et al. (1992): Since A is underdetermined, we may rewrite Eq. (17) by applying the QR decomposition to its transposed form as Since h is in the range (Q x ) and because range (Q x )⊥ range(Q y ), one has Q T y h=0 (which is of interest below). Moreover, we observe that the analysis window g is the sum of two orthogonal vectors with g 2 = x 2 + y 2 . Due to Eq. (17), Q x x=Q x (R T ) −1 v, but Q y y may depend on other constraints. When searching for the minimum norm solution, we simply set Q y y 2 = y 2 =0 and obtain which is nothing than g min = A T (AA T ) −1 v. However, for a meaningful interpretation of the Gabor expansion, we would prefer an analysis window g that is locally concentrated in the TF-plane. The design of such a function g when the synthesis function h as well as K and M are given is a nontrivial problem, which was addressed in Qian and Chen (1993) and Qian et al. (1992). The problem can be formulated as follows: Given an optimally concentrated function h (e.g. the preassigned synthesis function), find its biorthogonal function g whose shape best approximates time and frequency shifted versions of h, i.e. minimize while Ag=v. For fixed a and b, the optimal vector y in the representation for g (x is still fixed through the biorthogonality relation) is given by Choosing h a,b =h yields Q T y h a,b =0 (see above) and thus y=0 and consequently, g=g min , i.e. the shape of g min best approximates the shape of h. Therefore, the TF-plane localization properties of h carry over to g in this case. But note that in principle any function h a,b is allowed and thus there is a large variety of possible analysis atoms g.

Gabor representation of two examples
To illustrate the signal separation property of the discrete Gabor expansion for a single dwell, we consider two examples.   The method of Zrnić (1975) was first used to simulate a signal in line with the classical signal model, which contains only noise and a stationary atmospheric component. In the frequency domain, the atmospheric signal peak is assumed to be a Gaussian centered at f d =ω/2π=−10.9 s −1 and with a spectral width of w=0.9 s −1 . The discrete spectrogram of this signal is shown in Fig. 2. The atmospheric signal component is represented as a horizontal line (stationarity) centered at the prescribed Doppler frequency. Noise is spread over the complete TF plane. Now lets take a look at measured time series data containing an additional intermittent clutter component. This dataset is further discussed in Sect. 5. The original I/Q data is shown in Fig. 3. Clearly, this time series is not stationary but contains transient components due to migrating birds. Assuming that a time resolution of O(1s) is sufficient to resolve these transients, we select a time resolution of about 0.5 s for the Gabor expansion. This corresponds to a frequency resolution of about 2 Hz. An appropriate sampling density in the TF-plane is achieved with M=64 and K=64. Setting M=128 and K=128, we get an oversampling factor of 3.5; the optimal scaling is given by σ 2 h ≈1. In contrast to the simulated case, the spectrogram of the measured signal shown in Fig. 4 shows additional nonstationary signal components, which are a typical signature of contamination by intermittent clutter. Taking    of the signal it is difficult to identify the separate transients, which show up as maxima of the envelope of the I/Q signal. However, Fig. 4 shows the same signal, but now its Gabor phase-space representation. This clearly provides a far better picture of the signal transients, even if the spectrogram shows only the modulus of the Gabor coefficients (the Gabor coefficients itself are complex). Visible are three distinct transitory bird-events. Two of them overlap in time and can therefore not easily be distinguished in the time representation. All bird signals are much stronger in amplitude than the atmospheric signal of interest. The latter can be seen as a line at quasi-constant frequency, centered at about 3 Hz. By comparing Fig. 2 with the real data shown in Fig. 4, the goal of the filtering process becomes evident.

Motivation for the statistical approach
With the tool of the Gabor representation at hand, the next step is to derive an appropriate filtering strategy for removal of the transient clutter signals. Our intention is to use the available a-priori knowledge about the signal components (atmosphere, noise, clutter) to construct an objective decision process aiming at a proper signal component separation. It is well-justified that both the atmospheric and the noise signal component are stationary Gaussian random processes. The atmospheric signal has a bounded spectral width much smaller than Nyquist interval, whereas noise is white and spread over the full TF plane. Not much is known in contrast about intermittent clutter, only the non-property that this signal component is nonstationary over typical dwell-times. We make use of this a-priori information to derive a filter that has a pass-characteristics for realizations of wide-sense stationary random processes and a stop-characteristics for all nonstationary processes. That is, signals looking like the simulated example shown in Fig. 2 should not be affected by the filtering process. The goal is thus to derive an objective procedure, which modifies the Gabor phase space representation of signals in such a way, that stationary Gaussian signal components are preserved.
One can imagine several strategies for implementing such a filter. For instance, this could be based on image processing techniques or a fuzzy-logic approach similar to the one used by Cornman et al. (1998). We follow a simpler statistical approach, which has first been used by Merritt (1995) for the same problem, where it is applied to the temporal sequence of Doppler spectral coefficients at fixed frequency bins. The goal is to construct a similar test, but this time in Gabor phase space. We therefore need to analyze the statistical properties of the Gabor coefficients with respect to the different signal components, in order to distinguish between clear air and clutter return. This immediately leads to the question of how the properties of Gaussian stationary processes are mapped to the Gabor coefficients a m,k or |a m,k | 2 . This problem is discussed in the next paragraph.

Mean and variance estimator for Gabor spectrogram coefficients
Since we aim at constructing a statistical test (see the next section below) which is based on the expectation and the variance of the individual Gabor spectrogram coefficients |a m,k | 2 , we need to define adequate estimators for the expectation and the variance based on our observations (given through S). First, to simplify the notation, we introduce a λ as a shorthand notation of a m,k , i.e. in what follows we set λ=(m, k). Then, the Gabor spectrogram coefficients take the form As mentioned in the previous section, we assume that the data sequence S satisfies for all n=0, . . . , N −1,

ES[n] = 0 and ES[n]S[n
With these two assumptions, the expectation and the covariance of the Gabor spectrogram coefficients are given by , which is shown in Appendix C (Lemma 3 and Lemma 4). The " * "-symbol stands here for the discrete convolution. The latter two formulas show the influence of the dependency of S and the redundancy of the Gabor frame expansion. In case S would be i.i.d. (i.e. ρ[l] = δ l,0 ), it follows E|a λ | 2 = σ 2 and Cov(|a λ | 2 , |a η | 2 ) = σ 4 | g λ , g η | 2 .
If, moreover, {g λ } λ∈ forms an orthonormal system, the covariance matrix becomes diagonal; i.e. as long as we deal with a redundant frame, the Gabor spectrum is always correlated with a range of dependency described by the decay of the Gramian matrix of {g λ } λ∈ (up to the convolution with ρ). The essential observation for our purpose is Consequently, which holds true for independent as well as dependent samples S[n] that follow a distribution which is determined by its moments. As property (20) constrains only the first two moments, it may hold true for a much richer class of distributions (in particular, it holds true for normally distributed random variables). In order to construct a statistical test that verifies property (20), we have to find optimal estimators for E|a λ | 2 and Ann. Geophys., 26, 759-783, 2008 www.ann-geophys.net/26/759/2008/ Var|a λ | 2 that are based on a finite number of observations. To this end, we introduce an index subset λ ⊂ containing λ and L−1 further different indices η, i.e. | λ |=L. As an estimator for E|a λ | 2 =σ 2 ρ * g λ , g λ , which is based on L neighboring observation variables, we definê where the constant is given by For i.i.d. samples S[n], the correcting multiplier in estimator (21) reduces to C λ =| λ |=L, and therefore Eq. (21) is then nothing but the well-known mean estimator, Assuming there exists some small ε>0 with Lemmas 5 and 6 (see Appendix C) verify that Eq. (21) is a consistent estimator for E|a λ | 2 , i.e.
By the same reasoning, we define an estimator for variance, where the constant is defined by Similar as before, it is shown (see Lemma 7 in Appendix C) that estimator (22) is unbiased (and certainly consistent, but the proof is omitted). Switching to the i.i.d. case yields and therefore Eq. (22) simplifies tô which can be easily seen with the help of formula (C1). If, moreover, {g λ } forms an orthonormal basis, we end up with the classical variance estimator

A statistical test performing signal identification
After having established estimatorsÊ( λ ) andV ( λ ), we aim now for the construction of a test that identifies Gabor coefficients associated with intermittent clutter returns. Typically, an atmospheric return is stationary and assumed to follow a Gaussian distribution, i.e. a test on the first two moments of the signal will give us some indication if this is true.
The basic idea goes back to Merritt (1995), who statistically tested a sequence of single (non-averaged) Doppler spectra in to order decide whether a particular Fourier power spectrum coefficient was due to a Gaussian or non-Gaussian signal. For this, he used the classical test of Hildebrand and Sekhon (1974) in a modified way. Following this approach, we consider the squared modulus of the Gabor phase space coefficients, |a m,k | 2 . Because we are interested in stationary signal components, we consider the sequence |a m,k | 2 for fixed frequency bins, i.e. we just pick individual rows and let only the time index change. For a fixed frequency index k, we have to test the elements |a m,k | 2 to be of stationary Gaussian type. Typically, we assume the observed sequence |a m,k | 2 to be possibly affected by non-stationary intermittent clutter. Then, we will getÊ( λ ) 2 /V ( λ )<1. We also make use of the fact that intermittent clutter signals are almost always stronger than the (clear air) atmospheric return.
To identify intermittent clutter practically, we proceed as follows: in a first step, define the index set representing the k-th row, which we denote by k ={(m, k): m=0, . . . , M−1}, and sort for each k the sequence {|a m,k | 2 } (m,k)∈ k in decreasing order. That is, we derive the order statistic of {|a m,k | 2 } (m,k)∈ k which we denote by {|[a] m,k | 2 } (m,k)∈ k ([·] stands for the order statistic map). Therefore, we have |[a] m,k | 2 ≥|[a] m+1,k | 2 for all (m, k)∈ k . For l=0, . . . , M−1, we define subsets k (l)={(m, k) : m=l, . . . , M−1}. The largest coefficients are stepwise discarded, which has the goal of eliminating the clutter signal component. Using the quantitiesÊ( k (l)) and V ( k (l)) of the subset, the test statistics ϑ is computed for l=0, . . . , M−1 as long as holds. The largest coefficient of the first subset for which the test (positive for clutter) is not satisfied (a clutter-free subset) is then taken as a threshold for a frequency-dependent identification of the clutter component.

Signal separation through Gabor coefficient thresholding
All coefficients |a m,k | 2 greater than the threshold determined in the last section are regarded as clutter. One problem exists, if the subset k (l) becomes too small in this iterative process. Then the statistical estimate will become unstable and the estimation of a local threshold is no longer meaningful. This should not happen if the dwell time is sufficiently long, but it is not always known how long the dwell time must be for various types of intermittent clutter. Further investigations are needed to clarify this question. However, it might nevertheless be attempted to clean data sets regardless of the dwell time used. In such cases it can happen that some nonstationary components have a duration on the order of the dwell time. Then it can be useful to replace the local threshold with a non frequency-dependent global threshold, which could be derived from stable estimates of local thresholds at other frequencies k. Such a global threshold should be constructed in such a way, that it reflects the noise level in the Gabor representation. For instance, it could be estimated by averaging over a certain number of the smallest local thresholds. This method, however, has a risk of clipping also the atmospheric (clear-air) signal component. Leaving this problem aside, we can formulate the filtering procedure as follows: A coefficient |[a] l,k | 2 for which ϑ(|[a] l,k | 2 )≥1 holds is associated with clear air return. Based on the test, we introduce a clutter index set as The coefficients a m,k ∈ c k are finally set to t k e i arg a m,k , where t k is the average value of the remaining coefficients, The main result of this paper -the nonlinear filtering -is now formulated in the following: Let S be the given RWP signal. Based on our model assumptions, the filtered component is given by Finally, we discuss a practical aspect of the filtering method: The evaluation of the clutter index set c k requires the computation of the modified variance estimator. However, the computational overhead involved in calculating the modified variance estimator is obviously greater than in the case of the classical variance estimator. Our experience has shown that the variance estimates obtained with the two methods usually do not differ much. It may therefore be appropriate to use the classical variance estimator, if a saving of processing power is necessary for a real time implementation of the algorithm. This is left for a future study.

Data set
Now let us illustrate the performance of the proposed filtering algorithm by applying it to RWP data obtained with the 482 MHz wind profiler at Bayreuth, Germany on 13 October 2005. This radar is one of three operational systems that the Deutscher Wetterdienst currently uses in its aerological network. The technical characteristics are summarized in Table 1. More details and an overview of the standard signal processing steps are given in Lehmann et al. (2003). For wind measurements, the system is running in a fourbeam Doppler beam swinging configuration using two different pulse widths of 1667 ns (low mode) and 3333 ns (high mode). The averaging time for wind measurement is 26 min, another 4 min are used for RASS measurements of the virtual temperature. For the investigation of bird migration we consider only low mode data. The relevant sampling parameters are given in Table 2. Of interest are further the resolution of the time series t=0.007708 s, the number of data samples N =4608 and the total length or dwell time T d =N t=35.518464 s.
During the bird migration season in October of 2005, full time series data of the coherently integrated demodulated receiver voltage signal were saved in the wind low mode. Both wind and spectral data were manually reviewed to identify days with significant bird migration. It is well known, that a human expert can easily detect bird migration events by searching for typical patterns in the wind measurements (northeasterly directions in fall, discontinuities at sunrise and sunset), which are additionally accompanied by irregular and wide, sometimes multiple peaks in the Doppler spectra. In contrast to most clutter-free situations, those peaks often exhibit a poor time and range gate continuity. Time-height plots of the estimated moments (power, radial velocity and spectral width) are helpful to get a quick overview of potentially interesting cases, and a closer look into the time series data then typically confirms the conjecture of bird migration. Particulary significant bird migration was noted on 13 October and we therefore selected this day as a test case for the new bird mitigation algorithm. A significant fraction of this data was contaminated with bird returns; the effect is best seen in Fig. 10. Here, the winds have been computed without any intermittent clutter removal algorithm.
The consensus method is normally not able to suppress the effect of the bird echoes because of their frequent occurrence. The operationally used intermittent clutter removal algorithm (ICRA), a particular implementation of the statistical averaging method proposed by Merritt (1995), could only alleviate the problem, see Fig. 11. Also, the operational quality control (Weber-Wuertz continuity check, not shown) was only able to flag a small percentage of the contaminated data, because the erroneous wind data exhibited the typical intrinsic consistency.

Processing details and results
A software was developed for reading and writing of the profiler time series data using the proprietary binary data format. This made it easy to process the data using the Gabor filter and to save them again in the original file format. The reprocessed data could therefore be seamlessly integrated in the off-line version of the operational wind profiler software, to compare the performance of the different algorithms.  Fig. 36. Same as in Figure 33, but for the cleaned signal obtained from the filtered Gabor representation shown in Figure 35.   As an example, we consider again the measurement taken in the south beam of the profiler at range gate 9 (1.6 km height a.g.l.), with a start time of the dwell at 00:09:45 UTC. This time series was already discussed in Sect. 3. As described in Sect. 4, local (constant frequency) thresholds were estimated to separate the clutter part of the signal from the stationary components atmosphere and noise. During processing of the complete dataset it was revealed that the dwell time of about 35 s was apparently rather short to guarantee that every observed intermittent clutter signal exhibits a clear transient behavior. Sometimes the duration of the clutter signal component was on the order of the dwell time instead. If this is the case, then the estimation of the local threshold may become unstable and signal separation can partly fail with the result that clutter energy leaks through the filter. One way to remedy this problem is to replace local thresholds with a global threshold as described above. For the example data, this was done if more than 30 percent of the Gabor coefficients at a particular frequency were classified as clutter. The global threshold was then computed as the median over 15 percent of the smallest local thresholds, to get an estimate for the noise level. Another way to handle this situation would be to either flag this range gate as suspect or to replace the data with random white noise. Further research is needed to learn more about typical intermittent clutter characteristics and to optimize both the data sampling and the performance of the filter. The method described in this paper should be a useful tool for such investigations.
Application of the filtering strategy yields a filtered Gabor phase-space representation, which is shown in Fig. 5. Here, the moduli of the coefficients a m,k representing the transient (bird) contributions have been replaced by an estimation of the stationary signal component at that frequency (either noise or atmospheric signal). The reconstructed I/Q time series after back-transformation from the Gabor phase space domain is presented in Fig. 6. The nonstationary signal components have been suppressed and also the amplitude has been significantly reduced. It is easy to measure the reduction of total power by computing the difference in variance between the unfiltered and the filtered data, to get an information about how much clutter energy was removed by the filter.   Gabor filtering was performed for the complete dataset and the resulting bird-cleaned time series data were used for reprocessing of the whole day. This was compared with two additional processing methods: Method 1 used no intermittent clutter filtering algorithm, whereas method 2 used the routinely employed Intermittent Clutter Algorithm (ICRA), an implementation of the Statistical Averaging Method (SAM) originally proposed by Merritt (1995). The results for all range gates for the dwell taken at 00:09:45 UTC (stacked Doppler spectra) are shown in Figs. 7 (no filtering), 8 (ICRA filtering) and 9 (Gabor filtering). Without filtering, the lowest 17 range gates show spurious peaks and also large spectral widths due to the transient bird echoes. Note especially the discontinuity in height of the location of the estimated signal peak (derived Doppler velocity). With ICRA processing, the effect of the birds has been drastically re-duced, but there are still range gates which show spurious peaks. This indicates that ICRA was unable to reduce the clutter energy completely. Figure 9 shows the processing results of the newly suggested filtering algorithm. The spurious remnants of the bird clutter are almost completely gone, although range gates 15 and 16 (2.49 and 2.64 km height agl) show apparently some bird clutter energy leaking through. This is also reflected in the somewhat larger spectral width at these heights. However, the spectral peak is now continuous across all heights and the spectral width estimates are mostly unaffected by the clutter.
Finally, the horizontal wind vector data derived through the three different processing methods are shown in Figs. 10 (no clutter filtering), 11 (ICRA processing) and 12 (Gabor filtering), respectively. The color coding is due to the wind speed (magnitude of the horizontal wind vector). Obviously,  Fig. 38. Same as in Figure 37, but Doppler spectra were estimated using the operational bird-mitigation algorithm ICRA. Bird contamination below 3.0 km height is reduced compared to Figure 37, but still significant. Fig. 7, but Doppler spectra were estimated using the operational bird-mitigation algorithm ICRA. Bird contamination below 3.0 km height is reduced compared to Fig. 7, but still significant.

Fig. 8. Same as in
the clutter contamination has been drastically reduced by the new algorithm.

Conclusions
We have dealt with wind profiler signals obtained during bird migration and shown, how the signals can be decomposed into a time-frequency representation. A Gabor frame representation was used for a time-frequency analysis of the data and turned out to be a good method for signal-clutter separation. Previous attempts for intermittent clutter filtering have made use of the wavelet transform (WT) and its discrete versions (Jordan et al., 1997;Boisse et al., 1999;Lehmann and Teschke, 2001), so it is interesting to briefly discuss the difference between the dyadic wavelet and the Gabor approach, and to point out why we favor the Gabor method in comparison. The dyadic WT is another way of analyzing nonstationary signals. The difference lies in the tiling of the TF-plane by the elementary signals (or time-frequency atoms). In the Gabor (WFT) approach, the tiling is uniform with fixed resolution. This is in contrast to the wavelet approach, where the tiling is generally variable. For example, an orthonormal wavelet basis decomposes the frequency axis in dyadic intervals whose sizes grow exponentially. In other words, the frequency resolution gets worse the more the time resolution is improved. This is wanted if the signals under investigation have high-frequency components of short duration embedded within low-frequency components of slow temporal variation. For the RWP signals however, we found no evidence for such a behavior. The intermittent clutter components occur at nearly all frequencies within the typical Nyquist range, with no obvious difference in temporal characteristics. In particular, they can occur close to zero frequency where the temporal resolution of the WT is the worst. Especially in this case, the WT seems not to be ideal for resolving the transient nature of the intermittent clutter.  Fig. 39. Same as in Figure 37, but Doppler spectra were estimated after statistical Gabor filtering of the original time series. Only minor remnants of bird contamination can be seen in range gates 15 and 16 (at 2.5 and 2.6 km height). Fig. 9. Same as in Fig. 7, but Doppler spectra were estimated after statistical Gabor filtering of the original time series. Only minor remnants of bird contamination can be seen in range gates 15 and 16 (at 2.5 and 2.6 km height).
Examples of clutter signals in both representations shown by Justen and Lehmann (2003) illustrate this quite clearly. So the argument which is often used against the WFT, namely its constant time-frequency resolution, turns out to be advantageous. Additionally, the Gabor expansion using a Gaussian window achieves the best possible time-frequency resolution by reaching the lower limit of the Heisenberg uncertainty constraint.
To identify the clutter contribution to the signal, we make use of the a-priori information that the atmospheric signal component of interest can be adequately modelled as a stationary, proper complex Gaussian random process. Using this assumption, a test statistic is constructed to serve as a criterion for the discrimination between stationary and nonstationary signal components. This follows the approach first suggested by Merritt (1995). However, in case of the redundant Gabor transform it turns out, that the variance estimator has to be modified to guarantee its unbiasedness and consistency. Proofs for the necessary modifications are given in detail.
Finally, the algorithm has been applied to a dataset obtained with a 482 MHz wind profiler during bird migration. It could be demonstrated that the performance of the new algorithm was superior to the performance of the operationally used intermittent clutter reduction algorithm, without obvious negative side effects. Application of the algorithm has shown, that sampling settings of the wind profiler play an important role in the clutter mitigation capabilities of the algorithm. This is not unexpected, since both the sampling period and the dwell time determine the resolution of the Doppler spectrum and obviously also the resolution of any time-frequency representation. Furthermore, longer dwell times may ease the identification of transient clutter signals and the stable estimation of the thresholds for noise and the  . The x-axis shows time and the y-axis denotes height. Data have been color coded by wind speed. The signal processing was using no bird mitigation algorithm. Relatively strong northeasterly winds below about 3.5 km indicate strong bird migration, this can be seen between 00 and 05 UTC at heights around 1000 m and above 1600 m and especially after 18 UTC from the lowest gate to about 3500 m. . The x-axis shows time and the y-axis denotes height. Data have been color coded by wind speed. The signal processing was using no bird mitigation algorithm. Relatively strong northeasterly winds below about 3.5 km indicate strong bird migration, this can be seen between 00:00 and 05:00 UTC at heights around 1000 m and above 1600 m and especially after 18 UTC from the lowest gate to about 3500 m. stationary atmospheric component. This is especially important for cases where atmospheric and clutter signal overlap in frequency.
Future work is suggested for a better quantitative characterization of intermittent clutter signals during dense bird migration. This should allow to optimize both sampling and processing settings for operational wind profiler systems. A long-term evaluation of the new algorithm would be useful to determine its limits and to estimate the performance improvements of the new methods, in comparison with previously used algorithms. This would be facilitated by an online-implementation of the method and a means to compare the profiler wind measurement with independent data, e.g. radiosonde measurements.  Figure 310. The signal processing was using the standard ICRA algorithm. Bird contamination has been reduced compared to Figure 310, but is still significant after 19 UTC. A few other northeasterly wind barbs around 02 UTC are affected by intermittent clutter echoes. Fig. 11. Same as in Fig. 10. The signal processing was using the standard ICRA algorithm. Bird contamination has been reduced compared to Fig. 10, but is still significant after 19:00 UTC. A few other northeasterly wind barbs around 02:00 UTC are affected by intermittent clutter echoes.
The map, F : H → 2 , defined via F : f → { f, h λ } is usually referred to as the frame operator (analysis operator). So the signal is characterized by inner products with the frame. To answer the question of how f can again be synthesized from the inner products { f, h λ }, we consider the adjoint frame operator given by F * c= λ∈ c λ h λ . This allows us to write If F * F equals the identity Id, F * performs a perfect reconstruction. This is the case when {h λ } forms an orthonormal basis. However, in general one has to apply (F * F ) −1 to Eq. (A2). This is possible since the inverse exists and is bounded because of Eq. (A1), A · Id ≤ F * F ≤ B · Id and thus Since (F * F ) −1 is self-adjoint and denoting (F * F ) −1 h λ =g λ , one consequently has In frame lore, g λ is referred to as the canonical dual frame with respect to h λ . In general, (F * F ) −1 cannot be explicitly computed but must be approximated by an iterative approach. However, the situation can be essentially relaxed when assuming that the frames {h λ } and {g λ } form not a primal-dual, but a biorthogonal frame pair, i.e. h λ , g η =δ λ,η . IfF denotes the frame operator with respect to g λ , thenF =F (F * F ) −1 and one may write ( h λ , g η ) λ,η∈ =F F * , which is diagonal. Therefore,F is an analysis and F * a synthesis operator yielding perfect reconstruction (and vice versa, i.e. exchanging the roles ofF and F * ). If now the bi-orthogonality relation   Figure 310. The signal processing was using the new Gabor filter algorithm. Bird contamination has again been reduced compared to Figure 311. There are no indications of bird migration between 00 and 05 UTC, and only a few obvious outliers and missing data after 19 UTC.  Fig. 10. The signal processing was using the new Gabor filter algorithm. Bird contamination has again been reduced compared to Fig. 11. There are no indications of bird migration between 00:00 and 05:00 UTC, and only a few obvious outliers and missing data after 19:00 UTC. yields a way to derive g λ , the inverse of F * F needs not to be computed.

Biorthogonal discrete Gabor frame expansion
The following lemma can be retraced to its original form in Wexler and Raz (1990), it gives an explicit proof of the biorthogonality relation. Cov(a λ , a η )=σ 2 ρ * g λ , g η , where " * " denotes the discrete convolution.
The latter lemma states that the Gabor coefficients a λ turn into dependent random variables (even when ρ is a delta sequence, i.e. for independent samples of S). The range of dependency is determined by the sampling density in the time-frequency space and the range of dependency of S. In case S is a sequence of i.i.d. random variables, the dependency of a λ is fully characterized by the reproducing kernel g λ , g η .
After having verified the basic properties of the Gabor power coefficients, we prove that estimator (21) is consistent and that estimator (22) is unbiased (The proof of consistency is omitted, because this requires the computation of the 8thmixed moment).