Abstract

Analyzed data for numerical prediction can be effectively initialized by means of a digital filter. Computation time is reduced by using an optimal filter. The construction of optimal filters involves the solution of a nonlinear minimization problem using an iterative procedure. In this paper a simple filter based on the Dolph–Chebyshev window, which has properties similar to those of an optimal filter, is described. It is shown to be optimal for an appropriate choice of parameters. It has an explicit analytical expression and is easily implemented. Its effectiveness is demonstrated by application to Richardson’s forecast: the initial pressure tendency is reduced from 145 hPa per 6 h to −0.9 hPa per 6 h. Use of the filter is not restricted to initialization; it may also be applied as a weak constraint in four-dimensional data assimilation.

1. Introduction

To eliminate spurious high-frequency oscillations, the initial data for numerical weather prediction models must be modified to reduce gravity wave components to a realistic level. This process is called initialization. Of the many methods of initialization that have been developed, one of the simplest is based on digital filtering (Lynch 1990). In Lynch and Huang (1992, hereafter LH92) an adiabatic initialization is performed by carrying out two short model integrations, one forward and one backward from the initial time. For each model variable at each grid point and level, this produces a sequence of values centered on the initial time. Each sequence is processed with a simple low-pass filter, and the initialized data comprises the resulting values. In LH92 a filter based on the Fourier transform of an ideal frequency response function, modified by a Lanczos window (defined below) was used. It was found that a filter span of 6 h was required to achieve adequate suppression of spurious oscillations.

The generalization of the filtering procedure to account for diabatic effects was demonstrated in Huang and Lynch (1993, hereafter HL93). The idea is to integrate the model adiabatically backward for half the span and use the terminal values so obtained as initial data for a forward diabatic integration over the total span. The sequences of values produced by the forward integration are centered on the initial time, and they may be low-pass filtered to produce the initialized data. This paper also showed how a more efficient initialization is possible using an optimal filter: the total filter span was reduced from 6 to 3 h.

The theoretical background for the design of optimal filters is the Chebyshev alternation theorem (Oppenheim and Schafer 1989). The construction of the optimal filter involves the solution of a nonlinear minimization problem using an iterative procedure called the Remez exchange algorithm (Esch 1990). In this note we describe a simple filter based on the Dolph–Chebyshev window, which gives results similar to those achieved with the optimal filter but which is much simpler to implement.

2. The Dolph–Chebyshev window

The function to be described is constructed using the well-known Chebyshev polynomials and was first used by Dolph (1946) to solve the problem of designing a radio antenna having optimal directional characteristics (Kraus 1988). The Chebyshev polynomials are defined by the equations

 
formula

The following properties are easily derived from the definition: Tn(x) is an nth-order polynomial in x, even or odd accordingly as n is even or odd; Tn(x) has n zeros in the open interval (−1, +1) and n + 1 extrema in the closed interval [−1, +1]; Tn(x) oscillates between −1 and +1 for x in [−1, +1]; Tn(x) > 1 if x > 1; for large x, Tn(x) ≈ 2n−1xn.

Now consider the function defined in the frequency domain by

 
formula

where x0 > 1. Let θs be such that x0 cos(θs/2) = 1. As θ varies from 0 to θs, W(θ) falls from 1 to r = 1/T2M(x0). For θsθπ, W(θ) oscillates in the range ±r. Clearly, W(θ) is symmetric about the origin. Thus, the form of W(θ) is that required of the response function of a low-pass filter: a maximum at θ = 0 and small values as θ → ±π. The remarkable thing about W(θ) is that it has a finite Fourier transform. By means of the definition of Tn(x) and basic trigonometric identities, W(θ) can be written as a finite expansion

 
formula

The coefficients {wn} may be evaluated from the inverse transform

 
formula

where |n| ≤ M, N = 2M + 1, and θm = 2πm/N (Antoniou 1993). Since W(θ) is real and even, wn are also real and wn = wn. The weights {wn: −Mn ≤ +M} define the Dolph–Chebyshev or, for short, Dolph window. This window may be used to modify the coefficients of a low-pass filter, as was done with the Lanczos window in LH92, to reduce Gibbs oscillations. Alternatively, the window may be used directly as a low-pass filter, as will be described below.

3. Design of low-pass filter

There are several ways to specify the Dolph window. The order N = 2M + 1 and ripple ratio r may be chosen, where r is defined as the maximum amplitude in the stop band [θs, π]. Then the width of the main lobe [−θs, θs] can be computed from

 
formula

The resulting window has the minimum main-lobe width (i.e., minimum θs) for the given ripple ratio and order N. The amplitude response was calculated for filters of several orders N = 2M + 1 with r = 0.1. The results are plotted in Fig. 1; the response (dB) is defined as

 
δ = 20 log10|W(θ)|,

so the minimum attenuation in the stop band is 20 dB (the energy of components with frequencies in the interval [θs, π] is reduced to 1% or less). Note from (4) that x0 and θs depend on M. In Fig. 1 we see that the width of the pass band decreases as M increases; thus, the order may be chosen to obtain the required frequency selectivity.

Fig. 1.

Frequency response (dB) for Dolph filter with ripple ratio r = 0.1 for filter orders N = 2M + 1 with M = 2, 4, 8, and 16.

Fig. 1.

Frequency response (dB) for Dolph filter with ripple ratio r = 0.1 for filter orders N = 2M + 1 with M = 2, 4, 8, and 16.

An alternative procedure is to specify the filter order N = 2M + 1 and stop-band edge θs. Then x0 and r are obtained from

 
formula

The window thus defined has minimum ripple ratio for given main-lobe width and filter order. Let us suppose components with periods less than 3 h are to be eliminated (τs = 3 h) and the time step is Δt = 0.5 h. Then θs = 2πΔt/τs ≈ 1.05. The responses for filters of order 5, 7, and 9 (or M = 2, 3, and 4) are plotted in Fig. 2. We see that damping in the stop band increases with filter order and that a filter of order N = 7, or span T = 2MΔt = 3 h, attenuates high-frequency components by more than 20 dB.

Fig. 2.

Frequency response (dB) for Dolph filter with stop-band edge τs = 3 h (θs ≈ 1.05) for filter orders N = 2M + 1 with M = 2, 3, and 4.

Fig. 2.

Frequency response (dB) for Dolph filter with stop-band edge τs = 3 h (θs ≈ 1.05) for filter orders N = 2M + 1 with M = 2, 3, and 4.

Finally, we may specify the ripple ratio r and stop-band edge θs and determine what order filter is required to achieve these. Solving (5) for M and eliminating x0, we find

 
formula

For θsπ the denominator is close to θs/2. It follows immediately that the minimum required time span T = 2MΔt is given approximately by

 
formula

For τs = 3 h and r = 0.1, this yields Tmin ≈ 2.86 h.

4. Comparison with other window functions

For initialization, we wish to keep the filter span as short as possible to minimize computation. However, if the filter is to be used as a constraint in four-dimensional data assimilation, the time span is set to the period over which data is to be assimilated. As the span increases and, with it, the order N, closer approximation to the ideal square-wave frequency response should be possible. For the Fourier expansion, this is so; but the amplitude of the Gibbs oscillations does not diminish with increasing order of truncation, so windowing is still necessary. For the optimal filter discussed below, excellent approximation to the ideal is attainable for higher order: it is possible to limit the approximation error in the pass band 0 ≤ θθp. As the Dolph function is monotonic in the range 0 ≤ θθs it is not possible to guarantee a response in the pass band whose flatness increases sufficiently quickly with increasing order. Thus, the Dolph function may be unsuitable for direct use as a filter; however, it can be used in the same way as other windows, in combination with the truncated Fourier transform of the response function for an ideal low-pass filter, to control the Gibbs oscillations and achieve a high accuracy of approximation to the ideal.

The coefficients for an ideal low-pass filter with frequency cutoff θc are

 
formula

Application of a uniform window of length N = 2M + 1 corresponds to setting hn = 0 for |n| > M:

 
formula

Other windows wn are applied by multiplying hn pointwise by the window value, wnhn; this corresponds to convolution in the frequency domain. The Dolph–Chebyshev function can be used in this way. The highest frequency component present in H(θ), the response function corresponding to hn, is cosMθ. To remove it, a window with main-lobe width wider than the period of this component is required. The main-lobe width of the Dolph window is 2θs; thus, θs = 2π/M should give reasonable damping of Gibbs oscillations. We apply three Dolph windows with differing values of θs to the filter given by (7) with span T = 24 h, time step Δt = 0.5 h, and cutoff period τc = 6 h (so that M = 24 and θc ≈ 0.5); see Fig. 3. Clearly, the central value θs = 2π/M (or τs = T/2) gives adequate attenuation of high-frequency oscillations. The smaller value of θs is insufficiently damping while the larger value widens the transition band to an unacceptable degree.

Fig. 3.

Frequency response for low-pass filter with parameter values T = 24 h, Δt = 0.5 h, and τc = 6 h (M = 24 and θc ≈ 0.5) modified by a Dolph window with stop-band edge θs ∈ {π/M, 2π/M, 4π/M} or τs ∈{T, T/2, T/4}.

Fig. 3.

Frequency response for low-pass filter with parameter values T = 24 h, Δt = 0.5 h, and τc = 6 h (M = 24 and θc ≈ 0.5) modified by a Dolph window with stop-band edge θs ∈ {π/M, 2π/M, 4π/M} or τs ∈{T, T/2, T/4}.

A large number of windows are defined and compared in Harris (1978), where their advantages and demerits are discussed. The frequency responses for several windows are compared in Fig. 4. The Lanczos window

 
formula

was used in LH92. Another frequently used window, due to Hamming, is defined by

 
formula

with α = 0.54. In the cases depicted in Fig. 4, the parameter values are T = 24 h, Δt = 0.5 h, and τc = 6 h (so M = 24 and θc ≈ 0.5) and, for the Dolph window, τs = T/2. The Dolph and Hamming windows appear very similar in effect. Closer examination shows that the Dolph window gives better attenuation of high frequencies: the minimum damping for θ > θs is 63 dB, compared to 50 dB for the Hamming window.

Fig. 4.

Frequency response for low-pass filter with parameter values T = 24 h, Δt = 0.5 h, and τc = 6 h (M = 24 and θc ≈ 0.5) modified by uniform, Lanczos, Hamming, and Dolph windows.

Fig. 4.

Frequency response for low-pass filter with parameter values T = 24 h, Δt = 0.5 h, and τc = 6 h (M = 24 and θc ≈ 0.5) modified by uniform, Lanczos, Hamming, and Dolph windows.

5. Comparison with optimal filter

An optimal filter has the smallest maximum approximation errors in the pass and stop bands for a prescribed transition band. Once the filter order N, pass-band edge θp, and stop-band edge θs are given, the optimal filter coefficients may be calculated by an iterative numerical procedure. The examples below were generated using the code in McClellan et al. (1973). In Fig. 5 we compare the optimal filter response for four values of the pass-band edge, τp = 4 h, 6 h, 8 h, and 10 h with the Dolph filter. The fixed parameters are T = 3 h, Δt = 0.5 h, and τs = 3 h. Thus M = 3, θs ≈ 1.05, and θp varies from about 0.75 to about 0.3. The response of the optimal filter approaches that of the Dolph filter (with the same values of M and θs) as τp increases (or θp decreases): for τp = 10 h, the two curves are indistinguishable on the plot.

Fig. 5.

Frequency response (dB) for Dolph filter with M = 3 and τs = 3 h, and for optimal filters for τp = 4 h, 6 h, 8 h, and 10 h. Other parameters: T = 3 h, Δt = 0.5 h, and τs = 3 h.

Fig. 5.

Frequency response (dB) for Dolph filter with M = 3 and τs = 3 h, and for optimal filters for τp = 4 h, 6 h, 8 h, and 10 h. Other parameters: T = 3 h, Δt = 0.5 h, and τs = 3 h.

In HL93 the filter parameters used for the optimal filter were T = 3 h, Δt = 360 s, τs = 3 h, and τp = 15h (so M = 15, θs ≈ 0.2, and θp ≈ 0.04). The response of this filter (see Fig. 12 of HL93) was compared to the Dolph filter with the same values of M and θs: the results (not shown here) were, for practical purposes, identical.

To further illustrate the close similarity between the Dolph and optimal filters, Table 1 presents the filter coefficients {hn: 0 ≤ nM} for two filters. In each case the total span is T = 3 h and the time step Δt = 300 s, so that M = 18 and N = 37. The stop-band edge is τs = 3 h for each filter and the pass-band edge for the optimal filter is τp = 15 h. It can be seen from Table 1 that the coefficients agree to three significant figures: for practical purposes the two filters may be considered to be essentially equivalent.

Table 1.

Filter coefficients for a Dolph filter and for an optimal filter with τp = 15 h. In both cases the stop-band edge is τs = 3 h, the span T = 3 h and Δt = 300 s, so M = 18 and N = 37.

Filter coefficients for a Dolph filter and for an optimal filter with τp = 15 h. In both cases the stop-band edge is τs = 3 h, the span T = 3 h and Δt = 300 s, so M = 18 and N = 37.
Filter coefficients for a Dolph filter and for an optimal filter with τp = 15 h. In both cases the stop-band edge is τs = 3 h, the span T = 3 h and Δt = 300 s, so M = 18 and N = 37.

The optimal filter is more general than the Dolph filter: it can be designed to have multiple pass and stop bands and may have ripples in the pass bands. The Dolph window cannot replicate this behavior, as it is monotone in the interval [0, θs]. But for the parameter values of interest here the Dolph filter gives comparable results. Since the optimal filter is, by construction, the best possible solution to minimizing the maximum deviation from the ideal in the pass and stop bands, the Dolph filter shares this property provided the equivalence holds. In the appendix it is proved that the Dolph window is, in fact, an optimal filter whose pass-band edge, θp, is the solution of the equation W(θ) = 1 − r. Note the essential distinction: for the general optimal filter, θp can be freely chosen; for the Dolph window, it is determined by the other parameters. The algorithm for the optimal filter is complex, involving about one thousand lines of code; calculation of the Dolph filter coefficients is simplicity: the Chebyshev polynomials are easily generated from the recurrence relation, and the coefficients follow immediately from (3).

6. Application to initialization

An adiabatic initialization using a filter with weights {wn: −Mn ≤ + M} is performed by carrying out two short model integrations of length MΔt, one forward and one backward from the initial time. For each model variable x at each grid point and level, this provides a set of values {xn: −Mn ≤ + M}. The initialized values are then defined as

 
formula

This sum may be calculated cumulatively as the integrations proceed; full details are provided in LH92. Generalization to include diabatic effects is discussed in HL93.

In view of the practical indistinguishability of the Dolph filter and the optimal filter for the parameter values chosen in HL93 (see Table 1), we may expect that the results obtained by initializing with a Dolph filter would be virtually identical to those reported in that paper. Another application will be described here. Richardson (1922) calculated the pressure tendency using observations valid at 0700 UTC 20 May 1910. Richardson’s data tables have been extended using original sources, and a model based on his formulation of the primitive equations has been written (Lynch 1994). For the unmodified data, the initial pressure tendency at a central point calculated using the model was 145 hPa per 6 h, in agreement with Richardson’s value. When an initialization was performed using a Lanczos windowed filter with time step Δt = 300 s, cutoff period τc = 6 h, and span T = 6 h (as in LH92), the tendency was drastically reduced to a value of −2.3 hPa per 6h. The same data was initialized using a Dolph filter with a 3-h span and stop-band edge τs = 3 h; the filter coefficients are shown in Table 1. The initial pressure tendency was, in this case, further reduced to a value of −0.9 hPa per 6 h. Richardson reported observations for the date and time in question showing that the barometer was almost steady in the region of the central point. Thus, the value produced with the Dolph filter is the more realistic result.

The initial tendency at a central point is a useful indicator of local gravity wave activity. A more global measure is provided by the quantity N1, defined by

 
formula

which is related to the absolute pressure tendency averaged over the forecast area. The evolution of N1 for three forecasts is depicted in Fig. 6. In all cases the time step was Δt = 300 s. For uninitialized data, the solid line shows that N1 starts at a relatively high value of about 4 × 10−6 s−1 and falls to around 10−6 s−1 within 3 h. The two initialized runs, one using the Lanczos filter and one the Dolph filter, have much lower starting values for N1, and this parameter remains small throughout the forecasts. Although the span of the Dolph filter is only half that of the Lanczos filter, there is essentially no difference in noise levels between the two initialization methods.

Fig. 6.

Evolution of the noise parameter N1 (10−6 s−1) for the first 3 h of forecasts from uninitialized data (solid), and data initialized using the Lanczos filter (dashed) and the Dolph filter (dotted).

Fig. 6.

Evolution of the noise parameter N1 (10−6 s−1) for the first 3 h of forecasts from uninitialized data (solid), and data initialized using the Lanczos filter (dashed) and the Dolph filter (dotted).

7. Conclusions

The Dolph–Chebyshev window has properties similar to those of an optimal filter: it has been shown to be optimal for an appropriate choice of parameters. It has an explicit analytical expression, making it especially easy to implement. Its effectiveness has been demonstrated by application to the forecast first made by Richardson. The initial pressure tendency was reduced from an unrealistic to a reasonable level. The use of the filter is not restricted to initialization; it may also be applied as a weak constraint in four-dimensional data assimilation. Preliminary four-dimensional variational assimilation experiments with the HIRLAM model have been made using the Dolph window as a weak constraint on high-frequency components. The procedure was successful in controlling gravity waves. A comparison with the normal mode method is planned; this work will be reported elsewhere.

Acknowledgments

It is a pleasure to thank the following for helpful discussions and comments on an earlier draft: Jean-François Geleyn, Dominique Giard, Nils Gustafsson, Xiang-Yu Huang, Jan Erik Haugen, Erland Källén, Aidan McDonald, and two anonymous reviewers.

REFERENCES

REFERENCES
Antoniou, A., 1993: Digital Filters: Analysis, Design, and Applications. 2d ed. McGraw-Hill, 689 pp.
Dolph, C. L., 1946: A current distribution for broadside arrays which optimizes the relationship between beam width and side-lobe level. Proc. IRE, 34, 335–348.
Esch, R., 1990: Functional approximation. Handbook of Applied Mathematics, C. E. Pearson, Ed., Van Nostrand Reinhold, 928–987.
Harris, F. J., 1978: On the use of windows for harmonic analysis with the discrete Fourier transform. Proc. IEEE, 66, 51–83.
Huang, X.-Y., and P. Lynch, 1993: Diabatic digital filtering initialization: Application to the HIRLAM model. Mon. Wea. Rev., 121, 589–603.
Kraus, J. D., 1988: Antennas. 2d ed. McGraw-Hill, 892 pp.
Lynch, P., 1990: Initialization using a digital filter. Research Activities in Atmospheric and Ocean Modeling, CAS/JSC Working Group on Numerical Experimentation. Rep. 14, WMO, 1.5–1.6. [Available from Case Postale 2300, CH-1211 Geneva 2, Switzerland.]
.
——, 1994: Richardson’s marvellous forecast. Proc. Int. Symp. on the Life Cycles of Extratropical Cyclones. S. Grønås and M. A. Shapiro, Eds., 38–48. [Available from Geophysics Inst., University of Bergen, Allegaten 70 Bergen N-5007, Norway.]
.
——, and X.-Y. Huang, 1992: Initialization of the HIRLAM model using a digital filter. Mon. Wea. Rev., 120, 1019–1034.
McClellan, J. H., T. W. Parks, and L. R. Rabiner, 1973: A computer program for designing optimum FIR linear phase digital filters. IEEE Trans. Aud. Electro., 21, 506–526.
Oppenheim, A. V., and R. W. Schafer, 1989: Discrete-Time Signal Processing. Prentice-Hall International, Inc., 879 pp.
Richardson, L. F., 1922: Weather Prediction by Numerical Process. Cambridge University Press, 236 pp.

APPENDIX

Optimality of the Dolph Filter

Let J(θ) be the frequency response of an ideal low-pass filter with cutoff frequency θc. An optimal filter is defined by specifying pass and stop bands with edges θp and θs such that θp < θc < θs and minimizing the maximum deviation of the frequency response H(θ) from the ideal in these bands. Let ε(θ) = H(θ) − J(θ) and define the maximum deviation as

 
δ = max{|ε(θ)|: 0 ≤ θθp or θsθπ}.

The extreme points are those for which ε(θ) = ±δ. The function H(θ) can be written as a polynomial of order M:

 
formula

where x = cosθ. The Chebyshev alternation theorem (Oppenheim and Schafer 1989, 468–469) states that if H(θ) is a polynomial of order M that minimizes δ, there must be M + 2 extreme points with alternately positive and negative deviations. Moreover, H(θ) is unique and both θp and θs are extreme points.

The Dolph window W(θ) = rT2M [x0 cos(θ/2)] has M zeros in θs < θ < π. There are M + 1 extreme points in θsθπ (these include θs and π) for which W(θ) = ±r. If we define θp such that W(θp) = 1 − r, the point θp is a further extreme point, bringing the total to M + 2. Thus, W(θ) fulfils the conditions of the alternation theorem and must be the unique optimal solution with the pass-band edge given by θp.

Footnotes

Corresponding author address: Peter Lynch, Met Éireann, Glasnevin Hill, Dublin 9, Ireland.