## Abstract

Intermittent clutter signals are frequently observed by radar wind profilers during the seasonal bird migration. A novel statistical filtering algorithm based on a simultaneous time–frequency analysis of the profiler’s raw data was recently proposed to address shortcomings of existing methods. The foundation of this method is a Gabor frame expansion of the complex time series of the demodulated receiver voltage. In this paper, two objective criteria are suggested to obtain an optimal setup for the discrete Gabor frame expansion from the multitude of possibilities: first, the choice of almost-tight frames for a predefined maximum redundancy and second, the requirement that the analyzing bandwidth of the used Gaussian window function should provide a simultaneously sparse representation of both atmospheric signal components and intermittent clutter. The question of optimal sampling settings, especially dwell time, for a maximum reduction of bird interference is also discussed. Using data obtained during intense bird migration events it is shown that a combination of filtering and quality control of the result is required to prevent the occurrence of significant systematic and correlated errors in the final wind measurement.

## 1. Introduction

Radar wind profilers (RWPs) have become standard instruments for measuring wind in the atmosphere. Their data are used in numerical weather forecasting (Monna and Chadwick 1998; Bouttier 2001; Benjamin et al. 2004a; St-James and Laroche 2005; Ishihara et al. 2006), air quality monitoring (Dabberdt et al. 2004; White et al. 2006), and scientific studies and field campaigns (Fukao 2007). Reviews of the theoretical foundations of RWP have been given by Röttger and Larsen (1990), Woodman (1991), Gage and Gossard (2003), and Muschinski (2004). While further work is still necessary for a full theoretical understanding of the measurements, there are also difficulties outside of the scope of existing instrument theory, which need to be dealt with in practice. The most prominent issue here is certainly the clutter problem (Muschinski et al. 2005).

For UHF and L-band profiler systems, intermittent clutter echoes from birds during their seasonal migration in spring and fall can lead to significant correlated errors in the data. If such erroneous measurements are assimilated, then the quality of the forecasts is negatively impacted (Semple 2005; Cardinali 2009). This problem has been known for many years, see, for example, Ecklund et al. (1990), Barth et al. (1994), Wilczak et al. (1995, 1996), Engelbart et al. (1998), Richner and Kretzschmar (2001), Benjamin et al. (2004b), and Lehmann and Teschke (2008).

Meanwhile, the increased performance of radar processors allows for the application of rather sophisticated clutter detection and filtering algorithms. Recently, Lehmann and Teschke (2008) have suggested a Gabor frame decomposition of the raw data of the measurement, which is the time series of the demodulated (and, to some extent, coherently averaged) receiver voltage, followed by a statistical classification and filtering based on standard RWP signal model assumptions (Zrnić 1975, 1979; Woodman 1985; Lottman and Frehlich 1997). However, the Gabor frame expansion for a given dataset is not unique, with the main parameters being selectable subject to some constraints. The question of how to set up the algorithm for optimal filtering results was not fully addressed in Lehmann and Teschke (2008), and it is therefore the main subject of this article.

The paper is organized as follows: section 2 briefly reviews the Gabor filtering algorithm. In section 3, strategies for finding optimal discrete Gabor frame expansion parameters are discussed. This includes the selection of the discretization lattice, the bandwidth of the analyzing Gaussian window function, and the optimal dwell time. In the fourth section, a proposal is made for an additional quality control step after filtering, which is based on available a priori information of the desired atmospheric signal component. The goal is here to identify cases where filtering fails because of a breakdown of assumptions. Case studies illustrate the approaches discussed before. Finally, conclusions are given in section 5.

## 2. The Gabor filtering algorithm

### a. Gabor frame expansion

Frames are a generalization of the concept of bases in separable Hilbert spaces , for example, the space of functions of finite energy denoted with (Daubechies 1992; Allen and Mills 2004; Christensen 2008); they are useful for flexible linear decompositions of signals of the form

The set of expansion coefficients {*a _{i}*} for a particular basis or frame {

**b**

*} often provides a simpler and more concise description of a signal*

_{i}**f**, especially when only a few of the coefficients {

*a*} have significant absolute values (e.g., greater than zero). This is called a sparse representation (Mallat 2009). The discrete Gabor frame expansion has turned out to be useful for expanding radar wind profiler signals (Lehmann and Teschke 2008). Further details can be found in Wexler and Raz (1990) and Qian and Chen (1993).

_{i}Let **S** be a discrete and finite time signal with sampling points separated in time by Δ*T* at *n* = 0, … , *N* − 1, that is, *S*[*n*] = *S*(*n*Δ*T*). The signal is assumed to be *N* periodic, that is, . In finite dimensional signal spaces, the Gabor frame expansion leads to the discrete Gabor transform (DGT). It is essentially a discretized version of the continuous windowed Fourier transform, where the time–frequency plane (Gabor phase space) is sampled finitely at many points of a discrete lattice Λ (Daubechies 1990; Gröchenig 2001). Only regular rectangular grids described by the integer lattice parameters Δ*M*, Δ*K* (must be divisors of *N*) are considered, where Δ*M* denotes the time and Δ*K* is the frequency step size. If the set of vectors of translated and modulated windows **h**, with

for *m* = 0, …, *M* − 1, *M* = *N*/Δ*M*, and *k* = 0, …, *K* − 1, *K* = *N*/Δ*K* fully spans the vector space , then the discrete Gabor family is a frame for . The number of elements of the discrete Gabor family is

where *r* is called the redundancy factor of the frame. Here cannot be a frame for if *r* < 1, so *r* ≥ 1 is required. Only the oversampling case *r* > 1 is used in practice because of the Balian–Low obstruction theorem (Mallat 1999).

The discrete Gabor expansion of the signal **S** (synthesis equation) can be written as

where the Gabor coefficients *a*_{m,k} are obtained through inner products of the signal vector with the *MK* dual-frame vectors **g**_{m,k} (analysis equation),

The overbar denotes complex conjugation. Both relations are mappings between the transformed signal space, that is, the vector space of Gabor coefficients and the signal space . Inserting the analysis Eq. (5) in the synthesis Eq. (4) leads to the condition that the signal **S** can be recovered if

This is a discrete biorthogonality relation for the sequences **h** and **g**. This condition is satisfied if the dual-frame **g** fulfills the Wexler–Raz relation (Wexler and Raz 1990; Lehmann and Teschke 2008)

for 0 ≤ *p* ≤ Δ*M* − 1 and 0 ≤ *q* ≤ Δ*K* − 1. System (7) can be rewritten in matrix form: let **v** = [*N*/(*MK*), 0, … , 0]^{T} be a vector of length Δ*M*Δ*K* and **g** = [*g*[0], … , *g*[*N* − 1]] be the vector representing the discretely sampled dual frame, and let be a matrix with entries

then, the dual frame atom **g** is the solution of the linear system

Because of the restriction to the oversampling case *r* > 1, Δ*M*Δ*K* < *N.* As a consequence, the linear system (9) is underdetermined. A particular solution is obtained through the pseudoinverse of (Qian and Chen 1993); it is of minimum norm and is as close as possible to the primal window **h** (Qian et al. 1992; Lehmann and Teschke 2008),

The superscript *H* denotes Hermitian transposition. The window **h** is selected as a discrete Gaussian, because it provides the best possible localization in time and frequency (see, e.g., Mallat 1999). In the discrete Gabor frame expansion it is necessary to periodize this window (Wexler and Raz 1990, their appendix D) so

The width of this window controls the actual time–frequency resolution of the DGT through the parameter *s.* Furthermore, it is normalized to ensure .

### b. Statistical signal classification

Under the classical RWP signal model assumption, the raw signal *S* at the receiver output is the result of the demodulation of two independent real valued zero-mean and stationary Gaussian random processes describing the atmospheric signal and the receiver noise,

where and are independent complex (circularly) stationary Gaussian random vectors with zero mean and covariances **R*** _{I}* and

**R**

*(Zrnić 1979; Neeser and Massey 1993). The sequence has a vanishing pseudocovariance and is therefore called proper, that is, the expectation value*

_{N}*E*(

*S*[

*k*]

*S*[

*l*]) = 0. If

*S*[

*k*] is stationary, then the sampling interval is Δ

*T*and

*ω*denotes the mean Doppler frequency, and the covariance matrix has entries

Typically, a Gaussian correlation model is presumed for ϱ[*k*] (Zrnić 1979; Frehlich and Yadlowsky 1994), where *w* is the spectral width of the signal,

For identification of nonstationary signal components, which are assumed to represent intermittent clutter echoes, a statistical test based on the expectation and the variance of the individual Gabor spectrogram coefficients |*a _{m}*

_{,k}|

^{2}is constructed. With the shorthand notation

*λ*= (

*m*,

*k*), the Gabor spectrogram coefficients

*a*:=

_{λ}*a*

_{m}_{,k}take the form

Assuming that the data sequence **S** has zero mean *E*(*S*[*n*]) = 0 and autocovariance , the expectation and the covariance of the Gabor spectrogram coefficients are given by

Here, the asterisk symbol stands for the discrete convolution and the angled brackets denote the inner product. The essential observation is , so that

For a statistical test that verifies property (18), unbiased and consistent estimators for *E*(|*a _{λ}*|

^{2}) and Var(|

*a*|

_{λ}^{2}) based on a finite number of observations need to be constructed (Lehmann and Teschke 2008).

The test goes back to an idea of Merritt (1995), and it indicates whether the signal is stationary and Gaussian, and thus is likely of atmospheric origin. First, index sets denoted by Ω* _{k}* = {(

*m*,

*k*):

*m*= 0, … ,

*M*− 1} are introduced, and the sequence is sorted for each Ω

*in decreasing order. Because an intermittent clutter is almost always stronger than the (clear air) atmospheric return, the clutter part of the signal is removed by discarding the sorted Gabor coefficients (starting with the largest) in a stepwise fashion until*

_{k}holds for the test statistics. This procedure identifies the clutter subset . Then, |*a _{m}*

_{,k}| in the clutter subset is replaced with the average value in the clutter-free subset ,

The complete filtering operation can be written as

As a consequence of the nature of the statistical clutter identification it cannot be avoided that a few randomly large Gabor coefficients are erroneously identified as clutter even in a clutter-free dataset. The filtering step, which is always applied in the practical application, then marginally reduces the signal energy by about 0.3 dB (depending on the dataset).

The computational complexity of both the DGT and the Gabor filtering approach makes an online implementation feasible. For *N =* 8192, with a frame redundancy of 4 and 50 range gates, the full calculation takes only about 25% of the dwell time of 41.6 s on a standard PC with a 2.4-GHz clock frequency dual-core CPU. The implementation could be optimized further for speed through parallelization, because every range gate can be processed independently.

## 3. Optimality considerations

In the setup of the Gabor filtering a selection has to be made for the lattice parameters Δ*M*, Δ*K*, the width parameter *s* of the primal window, and the length *N* of the signal (or, equivalently, the dwell time *T _{d}*). It is not a priori clear how these parameters should be chosen to obtain the best possible results with the Gabor filter. This will be discussed in the following.

### a. Discretization lattice parameters

For a given discrete signal **S** of length *N*, a discrete Gaussian primal window **h** with fixed parameter *s* (the choice of an optimal *s* is discussed later), and the selection of the canonical dual window **g**, the remaining choice is the selection of the discretization lattice Λ. The principal requirements are that *r* > 1 and Δ*M*, Δ*K* must be divisors of *N.* Furthermore, an upper bound for the redundancy might be useful for the sake of saving computational costs (Lehmann and Teschke 2008).

The number of possible Gabor systems is finite for fixed *N* and *s*, and it is possible to use the additional constraint

for selecting the Gabor system that provides this minimum globally, that is over all possible lattices. Then, the analysis (or dual) window **g** is as close to **h** as possible, with the good localization properties of the primal (Gaussian) window **h** carrying over to the dual window. This can be achieved with frames that are almost tight (Daubechies 1992). If *A* and *B* are the frame bounds, then the following relation between primal and dual window is observed (Daubechies 1991):

By definition, *A* ≈ *B* for an almost-tight frame, and the redundancy of the frame is (*A* + *B*)/2 ≈ *A* for normalized frame vectors (Daubechies 1992; Mallat 1999). Then, **g** ≈ **h**/A, with and the analysis Eq. (5) can be written as

Inserting Eq. (24) into (4) leads to the orthogonal-like DGT (Qian and Chen 1993), where *a*_{m,k} would indeed be a measure for the similarity between the signal **S** and “the basis” **h**_{m}_{,k},

Although it is practical to limit the computational load through the requirement *r* ≤ *r*_{max}, a higher frame redundancy is useful, of course, in reducing effects of numerical noise in the Gabor coefficients (Mallat 1999). Selecting the highest affordable redundancy is therefore advantageous. Setting a maximal redundancy limits the number of possible canonical Gabor systems to be tested for finding the best lattice parameters Δ*M*, Δ*K*. The optimal pair can be easily found by computation of *E*(**g**, **h**) for all possible lattices. Note that *E*(**g**, **h**) does not depend on the data.

As an illustrative example, again take *N =* 8192, *s =* 0.5 (the selection of this parameter will be discussed in the next section), and *r*_{max} = 4. With the useful restrictions 2 ≤ *M*, *K* ≤ *N*/2, and *r* > 1, there are 21 different pairs (Δ*M*, Δ*K*) that generate admissible oversampled discrete Gabor systems . The minimum of *E*(**g**, **h**) is achieved for Δ*M* = 64 and Δ*K* = 32. In general, the optimal pair tends to have maximum redundancy, that is, *r* = *r*_{max}, although not always. For the RWP time series shown in Fig. 1, the corresponding Gabor spectrogram for the optimal lattice parameters (*M* = 128, *K =* 256) is shown in Fig. 2. Primal and dual-window functions are plotted in Fig. 3. The similarity in the shape of both windows is obvious, while the difference in absolute values reflects the redundancy of the frame [see Eq. (23)]. For a nonoptimal choice of *M* = 512, *K =* 64 (with all other parameters being equal), the dual window is no longer well localized (see Fig. 4). As a consequence, the spectrogram (Fig. 5) becomes blurred and the strong transient is smeared in time. While the essential features of the signal are still visible, it is evident that the selection of the optimal lattice parameters is much better adapted to the problem of separating signal and clutter.

### b. Analyzing window bandwidth–sparse representations

From a physical point of view, the most important parameter of the DGT is the joint time–frequency resolution controlled by the parameterized width *s* of the primal window. Wexler and Raz (1990) have suggested characterization of the time resolution through an effective width parameter defined as

where *h*_{max} is the maximum value of **h** with . The essential support of the primal window could also be used as a measure of time resolution (Gröchenig 2001), but any exact number is solely a matter of definition.

Intuitively it is clear that the time resolution must be “high enough” (the window “short enough”) to resolve the transient clutter echoes. However, a high time resolution comes at the expense of frequency resolution. Typically, the atmospheric signal has a finite spectral width *w* according to (14). Frequency resolution must therefore remain sufficient to resolve the atmospheric signal adequately.

As an example, consider the data shown in Fig. 1 again. The effects of a short primal window (high time resolution) with *s* = 10.0 corresponding to *T*_{1} = 0.15*s* is clearly visible in Fig. 6. What is remarkable is the appearance of a modulation signature in the clutter transient, which is likely due to the wing beat of the bird (Bruderer 1997; Zaugg et al. 2008). Another spectrogram was computed using a longer primal window (with low time resolution) with the parameters *s* = 0.01 corresponding to *T*_{1 }*=* 4.6*s* (see Fig. 7). While the wing-beat signature of the bird echo can no longer be seen, the better frequency resolution makes it easier to detect the clear-air atmospheric signal near a frequency of *f _{d} = −*2

*s*

^{−1}.

Obviously, the resolution in time and frequency has to be balanced for an optimal simultaneous representation of both atmospheric signal and clutter. If Ω ⊆ Λ denotes the subset of the lattice points that is influenced by the signal, then the ideal case would be a clear partition of Ω into the subsets Ω* ^{a}*, Ω

*, describing atmospheric signal and clutter components in a way that*

^{c}A nonoverlapping of both index sets describes the ideal situation, in which a perfect filtering is achievable. In such situations, the atmospheric signal would not be affected at all. Of course, this is merely a thought experiment, because it is impossible to investigate this question for real data: there is no objective method to classify whether a Gabor coefficient *a _{m}*

_{,k}indeed belongs to the atmospheric subset or the clutter subset or to both. Although this ideal situation can hardly be obtained in practice, it leads to the idea of selecting the Gabor representation where both the atmospheric signal and clutter are described with only a minimum of coefficients

*a*

_{m}_{,k}, or, in other words, a

*sparse representation*(Allen and Mills 2004; Mallat 2009).

Typically, sparse representations are constructed with a dictionary of elementary signals or atoms (see Chen et al. 2001). In this paper, only the Gabor frame representation is considered because it appears to be quite well adapted to the observed signals. It remains an open question whether it is possible to find signal representations that could achieve a better sparsity for the intermittent clutter problem.

The goal is then to find a sparse Gabor frame representation for both atmospheric and clutter components with a minimum of nonnoise Gabor coefficients *a*_{m,k} to reduce the chance of overlapping Ω* ^{a}* and Ω

*Note that no sparse representation can be found for noise.*

^{c}.The sparsity of a signal representation can be measured in a variety of ways (Hurley and Rickard 2008). A simple approach is to count the number of nonzero coefficients in the series expansion. However, in the presence of noise a slightly modified measure must be used,

where the number of all Gabor coefficients greater than a threshold *ε* is counted. Furthermore, is normalized by the number of frame vectors. The threshold is selected sufficiently large to avoid the counting of noise coefficients. The precise selection of this threshold is not critical as long as there is a significant difference between the noise and signal. For this reason, the following examples consider only wind profiler signals with a sufficient signal-to-noise ratio. This is typically the case for the 482-MHz RWP system at Bayreuth, Germany (Lehmann et al. 2003), at heights below 4 km. The estimate for *ε* was chosen to be one order of magnitude larger than the median value of the smallest 25% of the *MK* Gabor coefficients *a*_{m,k}.

Examples of values are given (see Figs. 2, 6, 7, and 11). Smaller values indicate a better adapted signal representation, represented by a better sparsity. To get an estimate for the average sparsity over a longer measurement period, the average value of was calculated for the following two datasets:

0100–0117 UTC 31 December 2008: clear-air signal and noise, but no clutter, and

2346–2354 UTC 29 March 2008: bird echoes during nocturnal migration in spring, with clear air and noise.

The sampling parameters of the radar are given in Table 1. Both datasets were manually reviewed: the first set contains no intermittent clutter returns (e.g., airplanes) in a cloud-free winter night, while the second set exhibits plenty of bird echoes confined to the lowest 25 range gates (with heights below 4 km AGL) during a nocturnal spring migration peak. The average value of the sparsity has been calculated for a primal window parameter range from *s =* 2^{−8} to *s =* 2^{5} in dyadic steps. The result is shown in Fig. 8 for the clear-air dataset and Fig. 9 for the clutter dataset. For the clear-air data, is constant for small width parameters *s* (for a long window with high-frequency resolution) until about *s =* 2^{−5}, and then increases monotonically toward larger values. This is an indication that the Doppler spectrum provides a sparse representation for stationary signals as long as the spectral peak is sufficiently resolved. Because of the finite spectral width of the clear-air signal, sparsity is not improved by a higher spectral resolution. A reduction of frequency resolution reduces sparsity, because the clear-air signal is stretched over more Gabor coefficients along the frequency axis. In contrast to this simple relation, the sparsity for the clutter case has a pronounced minimum for *s =* 2^{−1}. This indicates the time–frequency resolution, which achieves the optimally sparse representation in an average sense. Further investigations are necessary to see if these results are universal. Although it would be possible, in principle, to calculate the optimal *s* adaptively for a given dataset, the computational requirements are currently prohibitive.

### c. Sampling settings: Dwell time

An important RWP sampling parameter is the length of the time series, or dwell time, *T _{d}* =

*N*Δ

*T*, which essentially balances the time and frequency resolution of the Doppler spectrum. It is therefore appropriate to ask whether the Gabor filtering results can be optimized by a judicious setting of

*T*.

_{d}There are two conflicting requirements for dwell time. On one hand, it should be as short as possible to get a high time resolution of the measurement. Indeed, meaningful and valid measurement of the Doppler shift can be obtained with dwell times as short as 1 s, as discussed in Muschinski (2004). On the other hand, a better spectral resolution than the corresponding 1 Hz requires longer dwell times. Also, the statistical uncertainty of the estimation in the presence of noise (Woodman and Guillen 1974; Zrnić 1979; Woodman 1985) is reduced and the detectability of weak signals is improved by longer dwells. (Gage and Balsley 1978). In most wind profiling applications, dwell times of around 30 s are typical.

Theoretically, the dwell time is bounded from above by the requirement that the contribution of the dwell time to the estimated spectral width (resulting from intradwell frequency changes) should be negligible (Gage and Balsley 1978; Strauch et al. 1984; Gossard et al. 1998; White et al. 1999). Furthermore, it is bounded from below to allow for a high enough spectral resolution to adequately sample the spectral shape of the atmospheric signal peak (Woodman 1985; Wilfong et al. 1999) and other stationary signal peaks in the Doppler spectrum.

From the perspective of intermittent clutter detection, *T _{d}* must obviously be long enough so that the individual clutter transients can be identified as such. Merritt (1995, pp. 992, 994) has suggested that

The radar must dwell on each antenna beam long enough to allow moving objects sufficient time for their signals to change Doppler bins, angular position, and/or range gates. Therefore, a conservative approach would be to use the largest possible radar dwell time. (…) Observations to date … demonstrate that radar contamination from migrating birds is effectively reduced using somewhat longer dwells (1–2 min) than have been traditionally used (20–30 s).

During the spring migration season 2009, the Lindenberg 482-MHz wind profiler (Steinhagen et al. 1998) was configured to use a rather long dwell time of 147 s for the purpose of testing this assertion. The complete radar settings are given in Table 1. This dwell time should be sufficient to resolve the transient clutter echoes of single fliers in time–frequency space. Assuming that the horizontal speed of a bird is *O*(10 m s^{−1}) (Bruderer and Boldt 2001), then a typical (one way, half power) beamwidth of 3° corresponds to a lateral width of the beam from about 50 m at a height of 1 km to about 250 m at a height of 5 km. The crossing time of a single bird through the main scattering volume is therefore only *O*(10 s). Scattering in the sidelobes of the radar will only extend this time by an order of magnitude at most.

An instructive example from this dataset is shown in Fig. 10. The corresponding Gabor spectrogram (Fig. 11) shows quite a few energy maxima corresponding to several distinct clutter transients. This illustrates that a longer dwell merely increases the probability of having more than one flier in the radar scattering volume. While the detection of clutter transients is only a necessary condition for a successful filtering, the dwell must also contain a sufficiently long clutter-free time interval to allow for an unbiased estimate of *t _{k}* for all

*k*[see Eq. (20)]. Following an argument of Hildebrand and Sekhon (1974), about 50% of the Gabor coefficients

*a*

_{m}_{,k}for every fixed

*k*need to be uncontaminated to avoid erroneous estimates. Contrary to that, there is no contiguous time interval of more than 1 s without a bird signal component in the example data. The frequency of occurrence of clutter transients obviously depends on the density and distribution of birds in the radar scattering volume. These parameters can only be estimated in an elaborate way by using special radars and/or thermal infrared cameras (Holleman et al. 2008; Dokter et al. 2010; O’Neal et al. 2010). Such measurements were not available for the data used in this study, and an investigation of the quantitative relation between bird density and Gabor filtering efficiency between must therefore be left to future work.

For this particular example it could be surmised from the range gates above and below (not shown) that the true atmospheric signal was located at a frequency of about 3 Hz, for reasons of vertical wind profile continuity. However, such a signal is visible neither before nor after filtering (see Figs. 11 and 12). Although the energy of the strong transients has been significantly reduced, the signal energy in the frequency range from about −30 to +45 Hz is still larger than the energy of the atmospheric signal. This is due to a leakage of clutter signal energy through the filter. The example thus illustrates the limits of the proposed filtering method. Note, however, that this does not rule out a successful clutter detection. This will be discussed in the next section.

Longer dwell times are therefore no panacea to improve intermittent clutter filtering, at least not during dense bird migration. Consequently, the dwell time only needs to be long enough to resolve the nonstationary character of the clutter transients.

## 4. Quality control of the filtering

### a. An indicator of nonstationarity

As shown in the previous section, Gabor filtering is not always able to completely remove the bird contamination during extreme bird migration episodes. Filtering fails in cases where either the number of intermittent clutter transients is very high or where the dwell time is short compared to the duration of the transients. This makes it necessary to think about an additional quality control step. Fortunately, the algorithm can provide an indication about the nonstationarity of the signal, which would otherwise not easily be available. For each frequency index *k* in Gabor phase space, the test identifies an index subset that is assumed to be affected by intermittent clutter,

For fixed *k*, this condition will be fulfilled for *m* = 0, … , *m _{c}* − 1. Here,

*m*is the index of the largest clutter-free Gabor spectrogram coefficient in . For each range gate, a vector

_{c}**m**

*= [*

_{c}*m*(0),

_{c}*m*(1), … ,

_{c}*m*(

_{c}*K*− 1)]

^{T}of such indices can be derived. Information about the degree of nonstationary contamination can then be obtained through the normalized quantity

The parameter *β* tends to be very small for stationary signals. On the other hand, *β* → 1 is an indication that the statistical test (19) was not fulfilled because of nonstationarity. As a consequence, the filtered time series may still be affected by intermittent clutter, and additional a priori information must be used to quality control the data.

### b. A priori information about atmospheric signals

Because of the assumptions inherent in the signal model of Eqs. (12)–(14), moments of higher orders than two are usually not considered in profiler signal processing. In contrast, bird contamination leads to signal peaks that deviate significantly from the Gaussianity assumption (Wilczak et al. 1995; Kretzschmar et al. 2003). A test of the deviation from Gaussianity should therefore give additional information about the origin of the particular signal peak. Ideally, it should include higher-order moments like skewness and kurtosis because they should be good indicators of non-Gaussian signal peaks. For a first test, however, only the first three moments are considered.

Out of these, the spectral width *σ _{υ}* (Woodman 1985; Carter et al. 1995; Griesser and Richner 1998) contains the most promising information for the identification of intermittent clutter, because it is significantly enlarged by the bird transients. Otherwise,

*σ*is dependent on the spatiotemporal structure of wind and refractive index field for clear-air (Bragg) scattering, the size and velocity distribution of scattering particles for Rayleigh scattering, and the geometry (size) of the radar sampling volume. Unfortunately, statistical information about characteristic values of

_{υ}*σ*seems to be restricted to either smaller datasets or case studies (Williams et al. 1995; Ralph et al. 1995; White et al. 1996).

_{υ}### c. Combining spectral width with the nonstationarity indicator

The bird detection algorithm that is operationally used in the National Oceanic and Atmospheric Administration (NOAA) wind profiler network employs a spectral width–based thresholding (1 m s^{−1}) with considerable success (van de Kamp 1996). However, a large value of *σ _{υ}* is not, per se, an indication for intermittent clutter. Rayleigh scattering resulting from precipitation is known to generate wide signal peaks in Doppler spectra, depending on the actual drop size distribution (Wakasugi et al. 1986; Ralph et al. 1995; Williams et al. 1995; McDonald et al. 2004). This suggests a combination of

*σ*thresholding with the nonstationarity indicator

_{υ}*β*. A large value of the spectral width is considered to be indicative of bird clutter if and only if

*β*is above a certain threshold, that is, when the Gabor filtering step has indicated a high degree of nonstationarity of the time series. Through this combination, signal peaks with large values of

*σ*as observed in precipitation are not generally discarded.

_{υ}Results from a single case are presented in the following. To get an idea of clutter contamination, the winds are first processed without any attempt of clutter filtering. In the following it is shown that bird clutter has affected the height ranges below 5 km, with an artificial local wind maximum in a height band extending from about 1000 to 3000 m.

It would be useful to compare the profiler wind measurements with independent data with nearly the same vertical and temporal resolution, but unfortunately such data were not available. In such a situation, numerical weather prediction model data are often used in lieu of independent upper-air wind measurements to estimate the quality of a wind profiler (Steinhagen et al. 1994; Panagi et al. 2001; Hooper et al. 2008). A very short range numerical forecast of Deutscher Wetterdienst’s (DWD’s) Consortium for Small Scale Modeling-Europe (COSMO-EU) model (Steppeler et al. 2003; Schättler et al. 2011) is therefore taken to get some estimate of the “true” wind field and its extremes. Figure 13 shows the COSMO-EU wind forecast for Bayreuth plotted on exactly the same height and time scale as the profiler measurements. Model wind profiles were only available at hourly intervals and the vertical resolution is coarser than the corresponding resolution of the wind profiler. Although there are differences between model and measurement, as expected, it is nevertheless obvious that the model shows no local wind maximum below 5-km height in the time range of interest.

A comparison of the winds obtained without intermittent clutter filtering (Fig. 14) with the results from the operational intermittent clutter reduction algorithm (ICRA; Fig. 15) shows that the effects of the birds have been greatly reduced. However, the wind speed maxima below about 3000 m are still due to residual clutter. Gabor filtering (Fig. 16) has further reduced the clutter effects. However, artificial wind speed maxima also are visible in the last three profiles. The issue becomes obvious when the Doppler spectra of the filtered data (not shown) are inspected.

Finally, the results of Gabor filtering with the additional quality control step are shown in Fig. 17. The threshold for *σ _{υ}* was set to 2 m s

^{−1}and

*β*was selected as 0.5. Comparing this plot with Fig. 14 indicates that only data from the height band with the strongest clutter echoes were excised. Artificial wind speed values resulting from bird echoes are no longer obvious, but of course this comes at the expense of a decrease in data availability.

## 5. Conclusions and further work

The discrete Gabor frame expansion is a method for decomposing wind profiler raw data simultaneously in time and frequency, which allows a separation of stationary from nonstationary signal components. A statistical filtering method can then be constructed to identify and remove intermittent signal components from the data. This removal is generally independent of the acting scattering mechanism. While most of these nonstationary signal components are caused by intermittent clutter from birds, it has to be noted that either time-varying clear-air signals (e.g., in a strongly convective boundary layer) or intermittent precipitation will also be (partly) removed by the algorithm. This is not regarded as a disadvantage of the method, because nonstationary processes are not adequately described by a Doppler spectrum, so that relevant information is also lost in the case of no filtering.

The filtering results can be optimized by a judicious selection of the discrete Gabor frame expansion parameters in order to achieve the best possible separation in Gabor phase space between the atmospheric signal component and the intermittent clutter component. This can be realized by first selecting a discretization lattice for a prescribed maximum oversampling ratio in such a way that the resulting discrete frame is almost tight. Then, the dual-window function has a nearly identical shape compared to the primal (Gaussian) window, and both the analysis and the synthesis window have superior localization properties in the time–frequency plane. Second, the analysis bandwidth of the Gaussian window must be adapted to the data for obtaining a simultaneous sparse representation of both the atmospheric signal and clutter. For the test dataset it turned out that a maximally sparse representation could be achieved for a time resolution *T*_{1} of the primal window of about 0.5 s. Both optimizations improve the separability of the signal components, and thus make the filtering step more efficient.

There also appears to be a critical bird density beyond which clutter filtering fails and a retrieval of the clear-air atmospheric echo signals seems to be impossible. Such cases were observed during extreme migration events. An additional simple quality control step was therefore added to identify a breakdown of the assumptions inherent in the Gabor filter. This test can likely be improved through the consideration of moments of a higher order than two. Furthermore, a reduction of the size of the radar scattering volume should help to reduce the number of bird transients per time interval in the time series data.

A comprehensive objective evaluation of the algorithm performance is desired before the Gabor filtering method can be introduced into operational RWP systems. This would benefit from high-resolution reference data, ideally from other measurement systems like Doppler lidar (Grund et al. 2001; Shen et al. 2008; Pearson et al. 2009). A long-term statistical validation of the final wind estimates, after Gabor filtering, with independent data appears to be the best approach to assess the performance of the proposed method in the real world.

## Acknowledgments

The author wishes to express his appreciation to Gerd Teschke for help with the mathematical details of this method. Further thanks go to Raisa Lehtinen. Part of this work was supported by the European COST Action ES-0702, “European Ground-Based Observations of Essential Variables for Climate and Operational Meteorology.” The author benefitted from numerous discussions with many colleagues working on wind profiler radars over the last decade, in particular, Andreas Muschinski, Lutz Justen, Earl E. Gossard, Richard G. Strauch, Bob L. Weber, David A. Merritt, James R. Jordan, and John W. Neuschaefer. Special thanks go to Thomas Foken for his continuous support. The comments of the anonymous reviewers led to several improvements of the manuscript and are also gratefully acknowledged.