Chirp spectrum

Often (as in radar applications) a(t) is a slowly varying function of time and the phase θ(t) is oscillatory and varies rapidly, over the range of integration. With such waveforms, the stationary phase approximation can be used to investigate the spectrum.[6]^: 34[8][9][10] The method relies on the fact that the major contributions to the Fourier integral come from the region where the rate of change of phase is minimal, i.e. when

${\displaystyle {\frac {d}{dt}}[(\omega _{0}-\omega )t+\theta (t)]=0\qquad or\qquad (\omega -\omega _{0})-\theta '(t)=0}$

Unless θ(t) is a constant, the point in time t_s at which the phase is stationary will vary according to the instantaneous frequency ω_s.
Expressing the difference between (ω_s-ω₀).t and θ(t) as a Taylor series about the time t_s, but discarding all but the first three terms (of which the second term is zero, here), the Fourier integral can be written, approximately, as

${\displaystyle S(\omega )\approxeq a(t_{s})\int _{t_{s}-\delta }^{t_{s}+\delta }\exp \left[-j\left\{(\omega _{s}-\omega _{0})\cdot t-\theta (t_{s})-{\frac {\theta ''(t_{s})}{2}}\cdot (t-t_{s})^{2}\right\}\right]\cdot dt}$

In this equation t_s represents a constant time point, so terms depending on t_s alone may be taken outside the integral. The expression simplifies to[6]^: 39[10]
${\displaystyle S(\omega )\approxeq {\sqrt {2\pi }}\cdot {\frac {a(t_{s})}{\sqrt {|\theta ''(t)|}}}\cdot \exp \left[j\left\{(\omega _{0}-\omega _{s})t+\theta (t_{s})+{\frac {\pi }{4}}\right\}\right]}$
so
${\displaystyle |S(\omega _{t})|^{2}\approxeq 2\pi \cdot {\frac {a^{2}(t)}{|\theta ''(t)|}}}$

where ω_t is used to indicate the dependence of the frequency variable on t.
This is a very useful expression linking, as it does, the spectrum profile to the amplitude and phase characteristics of the chirp.

To carry out the inverse process, i.e. to find the time domain function s(t) given frequency domain data, the inverse Fourier transform is derived.

${\displaystyle s(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }|S(\omega )|\cdot \exp[j(\Phi (\omega )+\omega \cdot t)]\cdot d\omega }$

where Φ(x) is the phase function of the spectrum. The stationary phase points for this integrand are located at

${\displaystyle \Phi '(\omega )=-t}$

and the corollary relationship, equivalent to that derived for the spectrum, can be obtained by the stationary phase method, and is

${\displaystyle a^{2}(t_{\omega })\approxeq {\frac {1}{2\pi }}\cdot {\frac {|S(\omega )|^{2}}{|\Phi ''(\omega )|}}}$

In effect, stationary phase analysis gives the following (approximate) Fourier pair relationships:[6]^: 43
${\displaystyle a(t)\cdot \exp[j\theta (t)]\approxeq {\frac {1}{2\pi }}\cdot \int _{-\infty }^{\infty }S(\omega )|\cdot \exp[j\{\Phi (\omega )+\omega \cdot t\}]\cdot d\omega }$
and
${\displaystyle |S(\omega )|\cdot \exp[j\Phi (\omega )]\approxeq \int _{-\infty }^{\infty }a(t)\cdot \exp[-j\{\omega t-\theta (t)\}]dt}$

Consequently, approximate expressions for a(t) and θ(t) can be obtained when the spectrum, including its phase function Φ(ω) is given and, similarly, approximate expressions for |S(ω| and Φ(ω) can be obtained when the signal characteristics are given. Several examples of the procedure are given in the literature[6]^: 43[8][10]

Although the relationships are only approximate, their accuracy improves as the time-bandwidth product increases. In cases where the signal envelope and spectrum modulus are defined by smoothly varying Gaussian function then a T.ΔF product as low as 15 will give acceptable results, but if both a(t) and |S(ω)| are defined by rectangular functions, then the product T.ΔF needs to be much greater, typically over 100.[6]^: 49

Examples

Typically, in the radar case, a(t) is a constant over the duration of the signal and, for convenience, is assumed here to be unity. So the phase and amplitude characteristics, in the frequency domain, are related by

${\displaystyle \Phi ''(\omega )=\pm {\frac {1}{2\pi }}\cdot |S(\omega )|^{2}}$

There are two solutions for Φ(ω), which are complex conjugates of each other. The two filters with these characteristics can be used as the transmitter and receiver filters of a radar system and are interchangeable.
The group delay characteristic D(ω), (where D(ω)=-dΦ/dω), is

${\displaystyle D(\omega )=-\Phi '(\omega )=-\int _{0}^{\infty }\Phi ''(\omega )\cdot d\omega +K}$
so
${\displaystyle D(\omega )=\pm {\frac {1}{2\pi }}\cdot \int _{0}^{\infty }|S(\omega )|^{2}\cdot d\omega +K}$

So in the case of a rectangular time envelope, the dispersive delay characteristic is given by the integral of the square of the envelope.[10] If the positive sign is taken, then the group delay increases with increasing frequency and vice versa. The result is only approximate, but is more accurate for large values of the time bandwidth product.
Consider, as an example, the case of a spectrum that is uniform over the range -ω_max/2 to ω_max/2, then

${\displaystyle |S(\omega )|^{2}=A\qquad for\qquad |\omega |<{\frac {\omega _{max}}{2}}}$
so
${\displaystyle D(\omega )={\frac {1}{2\pi }}\cdot \int A\cdot d\omega +K={\frac {A}{2\pi }}\cdot \omega +K}$

Put D(-ω_max/2) = 0 and D(ω_max/2) = T, where T is the pulse duration, then K = T/2 and A = (2πT)/ω_max
so, finally
${\displaystyle D(\omega )=T\cdot \left[{\frac {1}{2}}+{\frac {\omega }{\omega _{max}}}\right]}$

As expected, a flat topped frequency spectrum corresponds to a linear frequency sweep.

The linear chirp is just one special case which, in any case, can be calculated more precisely by the methods of the earlier section. The particular usefulness of the stationary phase method lies in its ability to provide results when the frequency sweep is non linear. In such cases the spectral response can be shaped to meet some desired design criteria, for example, low side-lobes when a chirp is compressed. One such family of spectral functions that has been studied[6]^: 51 is given by

${\displaystyle |S(\omega )|^{2}=A_{n}\cdot cos^{n}\left({\frac {\pi \omega }{\omega _{max}}}\right)\qquad where\qquad |\omega |<{\frac {\omega _{max}}{2}}\qquad and\ n\ is\ an\ integer}$

It is possible to find the group delay characteristics of these functions in a similar manner to that carried out above and the results for n = 1 to 4 have been calculated.[6]^: 51
Although these cosine functions are amenable to mathematical manipulation, they are rarely chosen to define the spectral characteristics of a chirp, in practice, because when compressed they give broad main pulses with high side-lobe levels. A better characteristic (among many)[11] is the Hamming function, given by

${\displaystyle |S(\omega )|^{2}=A\cdot \left[0.54+0.46\cdot cos\left({\frac {2\pi \omega }{\omega _{max}}}\right)\right]=A\cdot \left[0.08+0.92\cdot cos^{2}\left({\frac {\pi \omega }{\omega _{max}}}\right)\right]}$

A plot of this characteristic is shown, plotted over the range -ω_max/2 to ω_max/2.

Applying the equations given above, the group delay characteristic which achieves this spectral shape can be obtained. It is

${\displaystyle D_{H}(\omega )=T\cdot \left[{\frac {1}{2}}+{\frac {\omega }{\omega _{max}}}+{\frac {1.7037}{4\pi }}\cdot sin\left({\frac {2\pi \omega }{\omega _{max}}}\right)\right]}$

Now because the principle of stationary phase shows that there is a direct relationship between elapsed time and the instantaneous signal delay then, for the Hamming window, t/T can be related to ω/ω_max by

${\displaystyle {\frac {t}{T}}={\frac {1}{2}}+{\frac {\omega }{\omega _{max}}}+{\frac {1.7037}{4\pi }}\cdot sin\left({\frac {2\pi \omega }{\omega _{max}}}\right)}$

This characteristic which is time as a function of frequency is shown here. Inverting the plot gives the more usual (and more useful) plot of frequency as a function of time, which is also shown.

Other spectral shapes can be investigated in the same way and the results, although approximate, are surprisingly accurate, especially when the time bandwidth product of the pulse is high.

The stationary phase method does not predict or deal with Fresnell ripples, so it is unable to offer any means by which these ripples can be minimized. As an example, the figure below shows a chirp spectrum with T.ΔF =250 obtained for a non-linear chirp aiming to match the Hamming window, using the methods described above. The figure shows that the spectral profile matches the Hamming characteristic quite well, but Fresnell ripples, not predicted by the method, are very much in evidence.

Whenever a Fourier integral cannot be evaluated by analytical means, an approximate solution is usually possible by numerical analysis. Such a procedure requires the function to be sampled, usually at equi-spaced intervals in time.

One consequence of sampling is that the resultant spectrum is periodic in the frequency domain. In addition to the (desired) baseband spectrum, additional versions of the spectrum occur, centered on multiples of the sampling frequency. To ensure that there is no overlapping of frequency data (i.e. no aliasing) the Nyquist sampling theorem must be satisfied. In practice, a sampling rate substantially higher than that dictated by the sampling theorem is advisable[12]^: 11

A straightforward way to approximate an integral, such as a Fourier integral, is to use the standard 'rectangle rule' for numerical integration. The method assumes the signal value taken at a sample instant remains constant for one sampling interval, until the next sample is taken. This procedure is sometimes referred to as a 'box-car generator', or a zero order sample and hold.[13]^: 114[14]^: 34 If the time interval between samples is W, then s_n = s(nW), and the desired integral is obtained, approximately, by summing the rectangular areas.

The result so obtained is the convolution of a rectangular pulse with step size W with the impulses located at the sampling instants with weights equal to the sample values.[12]^: 12 In consequence, the spectrum of interest will have superimposed upon it the frequency response of the sample and hold,[13]^: 135[14]^: 36 and the spectrum of the sampled signal Ss is given by:[12]^: 12

${\displaystyle Ss(\omega ))=W{\frac {\sin(\omega W/2)}{\omega W/2}}\cdot \left[\sum _{n=-\infty }^{\infty }s_{n}\cdot \exp(-jn\omega W)\right]}$

The first part of the expression, i.e. the 'sin(x)/x' part, is the frequency response of the sample and hold. Its amplitude decreases with frequency and it falls to 63% of its peak value at half the sampling frequency and it is zero at multiples of that frequency (since f_s =1/W).

The second term in the equation is called the Fourier transform of the discrete signal s_n.[12]^: 12[15] It is a continuous function over all ω and involves an infinite number of summations. In practice the summation process can be truncated to a finite number of samples, N, possibly because the waveform is periodic or zero outside the range of samples. Furthermore, because the same spectrum is endlessly repeated, it is possible to confine interest to spectral data within the range -ω_s/2 to +ω_s/2.

As an example, an exponential chirp (with its top frequency well below the Nyquist limit) is sampled at 256 points, as shown.

The sampled spectrum, Ss(ω) of this waveform, calculated using the equation given above, is shown. To simplify the plot, only the results at positive frequencies have been displayed. The influence of the frequency spectrum of the zero order hold circuit is clearly seen in the diagram.

The baseband portion of the spectrum is shown in more detail in the next figure and the response shows a distinct slope, being significantly lower at the higher frequencies.

Although the characteristic of the zero order hold has a small influence on this result, the slope is mainly due to the properties of the chirp. The waveform sweeps relatively quickly over the high frequencies and spends more time sweeping the low frequencies, consequently there is less energy content at the high frequencies with more at the lower ones. (A linear chirp, on the other hand, has a nominally flat spectrum because its frequencies are swept at the same rate, as shown in some earlier plots).

If we limit interest in the output spectrum to a finite number of discrete data points (= N), at frequencies ω_m given by
${\displaystyle \omega _{m}={\frac {2\pi m}{NW}}\qquad for\qquad m=0,1,2,...N-1}$

then the formula for calculating the discrete Fourier transform is

${\displaystyle Ss_{m}\quad =Ss\left(j{\frac {2\pi m}{NW}}\right)\quad =\quad \sum _{n=0}^{N-1}s_{n}\cdot exp(-j\left({\frac {2\pi mn}{N}}\right)}$

The calculations can be carried out by means of a straightforward computer algorithm,[12]^: 21 but this is not very efficient in computer usage. Consequently, more efficient algorithms have been developed, especially Fast Fourier Transforms (FFT). Computer programs which implement the FFT are widely available in the literature[12]^: 54[15]^{: 119, 412}[16] and in proprietary CAD programs such as Mathcad, MATLAB, and Mathematica.
In the following example a linear chirp with time bandwidth product of 25 is sampled at 128 points (i.e. N = 128). In the figure samples of the real part of the waveform are shown - note that these are samples in the time domain. The FFT process assumes the waveform is cyclic, so these 128 data points can be considered to be part of an endlessly repeating sequence in time.

Linear chirp with TB=25 and N=128

By calculating the N-point FFT of this data, the discrete spectrum of the sequence is obtained. The magnitude of this spectrum is shown in the attached figure, where these data points are samples in frequency. The data is cyclic so, in the plot, the zero frequency point is at n = 0 and also at n = 128 (i.e. both points are the same frequency). The point n = 64 corresponds to +fs/2 (and also to -fs/2).

Spectrum of Linear Chirp, TB=25, N=128

To display the spectrum in more detail (but not necessarily with more resolution[17]), the time sequence can be extended by zero padding.[15]^: 80–85[18][19] For example, extending the 128 point time sequence with zeros to give N = 4096 results in that part of the spectrum originally presented in 16 samples, now being presented in 512 samples, as shown.

The Fresnel ripples on a chirp spectrum are very obtrusive, especially when time-bandwidth products are low (under 50, say) and their presence leads to high time sidelobe levels when chirps are subject to pulse compression as in radar and sonar systems. They arise because of the sudden discontinuities in the chirp waveform at the commencement and termination of the pulse.

Although there are a number of procedures that can be applied to reduce the ripple levels, they are not all equally effective. Furthermore, some of the methods require amplitude shaping, or amplitude modulation, of the chirp pulse and this makes those methods unsuitable when, for example, the chirp pulses are to be transmitted by a power amplifier operating in a near-limiting condition. For such systems only the methods using frequency (or phase) pre-distortion are appropriate.

If the transitions at the start and end of the chirp are made less sudden (or more 'rounded'), then a reduction in ripple amplitude is achieved.[6]^: 213[20][21] The durations of the two transition regions need only be a small fraction of the pulse duration, and suggested values are between 2/ΔF and 3/ΔF [20] but, as expected, when the time-bandwidth product of the pulse is small, longer transition periods are needed. The actual profiles of these rise and fall regions of a pulse do not seem to be critical and may be provided, for example, by band limiting filters in analogue implementations and a linear slope in digital ones.

Two examples show the spectra of linear chirps with finite rise-times. The first is for a chirp with time-bandwidth of 250, where the rise and fall times are 4% of the total pulse duration and the second is for a chirp with time-bandwidth of 25, where the rise and fall times are 10% of the total. These two spectra show a marked reduction in ripple amplitude compared to the spectra of unmodified linear chirps shown earlier.

A analogous technique can be applied to the frequency characteristic of the chirp waveform by adding linear FM distortion segments (quadratic phase modulation distortion) to the frequency characteristic of the chirp, as shown. The method is effective because amplitude and phase distortions having functional similarity can produce similar effects when the distortion factors are small.[20][22]

This added distortion is referred to as 'predistortion'. Suggested values for these distortion regions, to achieve good results, were initially determined to be:
${\displaystyle \Delta f_{p}=0.75\Delta F\qquad and\qquad \delta =1/\Delta F}$

In 1988, later work[23] proposed slightly different values, namely:
${\displaystyle \Delta f_{p}=0.73\Delta F\qquad and\qquad \delta =0.86/\Delta F}$

The outcomes can be further improved by optimizing values for each particular situation.

Below are two plots showing the effects of frequency pre-correction, for time bandwidths of 25 and 250. These can be compared to the results in the earlier sections.

The ripple reduction achieved by frequency pre-correction, although significant, is seen to be less successful than that achieved by the amplitude modulation methods of the previous section. However, it has been suggested[21] that by implementing cubic (rather than quadratic) phase pre-correction, comparable results can be achieved.

This method uses an inverse Fourier transform in order to derive a waveform which has a spectrum with the phase characteristic of a chosen chirp but a new amplitude profile which is rectangular and ripple free. The method is very effective but, unfortunately, the waveform that is so derived has a semi-infinite time duration. If, for convenience, the newly derived waveform is truncated to a practical length, then some ripple is reintroduced onto the spectrum.

As an example, a linear chirp waveform with a time bandwidth of 25 is shown together with its spectrum magnitude (shown by a full line) which, as demonstrated earlier, has a large ripple component. It is possible to find, by means of an inverse FFT, a chirp waveform which, in the frequency domain, has the same phase characteristic as before, but with the rectangular magnitude characteristic shown by the dashed line on the plot. The chirp waveform resulting from this process has a very long time duration, but when it is truncated to say, a length 2T, then the spectrum acquires some ripple once more, as shown.

There are many applications in which a spectrum with a rectangular magnitude profile is not ideal. For example, when a chirp waveform is compressed by means of its matched filter, then the resultant waveform approximates to the sinc function and, consequently, has annoyingly high sidelobes. Often, to improve the characteristics of the pulse and lower the sidelobe levels, its spectrum is modified, typically to a bell-shaped profile. Similar problems arise in digital signal processing where the spectral shaping is provided by a window function, a process sometimes called apodization. In the case of an antenna array, similar profiling by "weighting functions" is used to reduce the spatial sidelobes of the radiation pattern.

Although spectral shaping of a chirp could be applied in the frequency domain, better results are obtained if the shaping is carried out in the time domain.[24][25] Examples of this process are shown for linear chirps with time-bandwidth products of 250 and 25. They have been shaped by a 3-term Blackman-Harris window[11] given by
${\displaystyle |U(\omega )|^{2}=0.42323-0.49755\cos \left({\frac {2\pi \omega }{\omega _{max}}}\right)+0.07922\cos \left({\frac {4\pi \omega }{\omega _{max}}}\right)}$
The spectra, now bell-shaped, are seen to be free of ripples.

Non-linear chirps can be devised that have a bell shaped spectrum, such as the Blackman-Harris window just discussed, and consequently will exhibit reduced ripple compared to the linear chirp. By means of the stationary phase method described earlier, an approximate relationship between time and frequency can be obtained:

${\displaystyle {\frac {t}{T}}={\frac {1}{2}}+{\frac {\omega }{\omega _{max}}}+0.1871\sin \left({\frac {2\pi \omega }{\omega _{max}}}\right)+0.014895\sin \left({\frac {4\pi \omega }{\omega _{max}}}\right)}$
Rearranging the equation, one can plot frequency against time, which looks roughly as shown to the right.

Non-linear chirps with spectral profiles of a Blackman-Harris window show ripple reduction, but can have disappointing performance due to amplitude profiles with fast rise and fall times.

As examples, below are two plots of the spectral magnitudes of such chirps, with time-bandwidth products of 250 and 25, respectively. These chirps have reduced energy content in their outer frequency regions, but still show less than ideal performance.