Our natural world behaves mostly analogically, at least at a macroscopic level, therefore to interface it to digital computers one has to sample analog signals. The universally known sampling theorem, credited to Nyquist and Shannon, but the story is more articulated  , states that to reconstruct correctly the information carried by a bandlimited signal, the sampling frequency must be at least twice the highest frequency of the signal : . In practice we use because there could be ambiguity in reconstructing the component associated to and, depending on the non ideal shape of real-world bandlimited spectra, also folding of the upper part of the spectrum. In this form the theorem is always valid but sometimes it is stated as where is the bandwidth of the signal. The last formulation implicitly assumes that the lowest frequency of the signal is zero, , otherwise it is not generally true, as we will see. If the sampling frequency is lower than , i.e. , each of the frequencies above will be aliased, i.e. superimposed or confused, in particular, with a corresponding frequency in the range (page 1 of Fig.1) according to the relation or with integer (Fig.1).
Fig.1 Folding around the Nyquist frequency and its multiples .
You may think of the diagram of Fig.1 as pages of length that fold over each other, in particular over the first, alternately like an “accordion-pleated”  strip of paper. Because this phenomenon the spectrum of a signal, when sampled, will be aliased or replicated over all the pages. When the spectrum of the signal is contained entirely in one of the pages, the spectral aliases of the sampled signal will not overlap. If the original spectrum is not on the first page, one of the aliases will be positioned on the first page with the result of having converted down the frequencies of the original spectrum without modification of the bins power. The order of the bins of the original spectrum will be preserved if the original spectrum of the signal is contained in the odd pages and inverted for the even pages. All this occurs because sampling in time domain is a multiplication of the signal by a comb of unitary pulses, which in frequency domain becomes a convolution of the Fourier transformed unitary pulses with the spectrum of the signal. A nice detailed explanation, both mathematical and visual, is given in . The complete spectrum of a sampled bandlimited signal is constituted of replicas of the original spectrum symmetrically disposed around multiples, positive and negative, of the sampling frequency, as illustrated in Fig.2. For example, let the original spectrum (the diagram on the upper side) be composed of two pure tones at 1 and 12 (arbitrary frequency units) and of a continuous spectrum ranging from 6 to 7. Sampling the signal at 5 we obtain the coloured diagram on the bottom, where we can see the symmetrical replicas of the original spectrum centred at multiples of and also note that the pure tones, originally positioned at 5 and 12, and then isolated, are now superimposed to the borders of the continuous spectrum. Of course, when we sample a signal with continuous spectrum for a finite duration , the spectrum of the sampled signal will be constituted of discrete components (bins), whose frequency resolution is , which are not visible in Fig.2.
Fig.2 Replicas of the spectrum of a sampled bandlimited or bandpass signal.
The “accordion-pleated” paper model leads to a straightforward mathematical formulation.
Let be the bandwidth of the signal to be sampled.
which means that the original spectrum of the signal must be contained entirely in one of the pages of Fig.1. Actually the spectrum could be segmented in different pages. In that case we have to state the above conditions for each segment and others have to be verified so that the segments will not fold on each other when the signal is sampled.
The connections among the page number, and are: and where is the floor function, i.e. the largest integer not greater than “arg”.
Isolating from the above inequalities we have
and eliminating we obtain from which
These are the fundamental formulae for undersampling (bandpass sampling, subsampling, downsampling, sub-Nyquist, super-Nyquist, harmonic sampling) even if I did not use any of these terms in my original report  because, at that time, I was not concerned about any specific terminology for this kind of operation.
For example, let it be and . Applying (2) we have , i.e. , and then from (1) we obtain all the allowable sampling frequencies: : , : and, of course, : . As anticipated, the order of the bins of the aliased spectrum of the bandpass signal is reversed or not depending on the position of the original spectrum of the signal to respect to the chosen : if the corresponding is odd the order is preserved, if it is even the order is reversed.
Note that for doing a correct undersampling you cannot use all the sampling frequencies , in fact in the above example it would be which is clearly wrong. If you are interested only in the lowest bound of the sampling frequencies, substituting (2) into (1) the left term furnishes
which is the same as that reported in . An expression equivalent to (1), but in the time domain, is given in . We did not use the equal sign in (1) and (2) to avoid a possible folding and ambiguities of the frequencies at the borders of the spectrum, but if the power in the bins outside the open range is zero or negligible we may write and , i.e. we have to consider the shape of the real-world bandpass spectrum and choose so that to avoid ambiguities and minimize the folding of the spectrum on itself. It is possible to give to (3) a different form. Let it be and call it “band index”. Consider . The range of is , furthermore it is and then . The range is and its diagram is shown on Fig.3. Note that for not .
Fig.3 Lowest bound of the sampling frequency normalized to the bandwidth versus the band index.
In a recent article  were reported the two formulae and to compute an allowable sampling frequency for undersampling a bandpass signal, being the carrier, and the bandwidth of the signal. The procedure to compute is: as a first approximation put , insert this value in then use this rounded-down integer value of to calculate the true . I think that this method, as an illustration of the undersampling concepts, is useless and even misleading at least for two reasons: the first because it is not simpler than the more general approach given by the inequalities (1) and (2), the second and worse, because it gives only a single sampling frequency instead of all the permitted frequencies, and the computed frequency is not even the lowest bound, but only that particular sampling frequency for which the spectral aliases of the sampled signal are centred on the pages of Fig.1. The above formulae express this last property in a foggy way and even as an algorithm they are twisted, compared to the sunny logic of (1) e (2). In fact take and , being and , as in the preceding example, we will have , instead the lowest bound for the sampling frequency is . If the only data at our disposal are and it is easy to switch to the general method taking and and carrying on the computation as I suggested. Anyway, the formulae and can be easily deduced from . By definition , therefore is always satisfied, because implicitly contained in (1): note that you cannot use every for undersampling, as already shown, because you have to satisfy the other constraint. To deduce , assume , i.e. the arithmetic mean of the two bounds of (1). Doing the following sequential manipulations we have: .
Because , and substituting deduced from (1), we finally obtain . The substitution produces, as expected, an lower than the arithmetic mean assumed in advance. Eventually it is and then with . A more perspicuous deduction can be done from the key conditions  considering that the aim of the above formulae for is to centre the spectral aliases of the sampled signal on the pages of Fig.1. It has to be from which .
Even in practical applications it is important to be able to calculate all the permitted sampling frequencies, because you may have some constraints that force you to choose a particular range of sampling frequencies, so that it is better to rely on the general method for this computation.
My interest in signal processing started in the mid of 1975 when I began doing my thesis in Physics  which consisted in the design and in the realization of a SODAR system for use in atmospheric boundary layer studies. For the hardware I basically followed the work done by E.J.Owens , adding some original solutions, anyway I was the first in Italy to design and build a SODAR system that really worked, and even today many people use my scientific and technical ideas and solutions, some of which are described in   , even if not all of them recognize it. SODAR is the acronym of SOund Detecting And Ranging and it is an acoustic RADAR which emits in the atmosphere short acoustic bursts which are scattered by turbulence. The repetition rate of the bursts is, normally, between 3 and 6 seconds depending on the spatial range to be explored. The emitter is usually a power loudspeaker, placed at the focus of a microwave parabolic dish generally enclosed in an acoustic shield to attenuate the environmental noise. There are also SODARs whose antennas are constituted of an array of loudspeakers. To respect to the frequency range emitted there are, basically, two kind of instruments: low and high frequency SODARs. Typically, the low frequency range is 1000-3000 Hz, the high frequency 5000-7000 Hz. In Fig.4 and Fig.5 are shown the antennas of recent versions of SODARs I contributed to design: the big truncated cone-shaped antennas of Fig. 4 are for low-frequency range, the three horn-reflectors of Fig.5 are for high-frequency range.
Fig.4 Antennas of the low-frequency SODAR system
Fig.5 Antennas of the high-frequency SODAR system
During the 1976, and for many years after, the first version of the SODAR, and its upgrades, I designed and built personally, were extensively used in measurement campaigns and there emerged the need of an efficient sampling of the signal and the necessity of a real-time processing of the data. The first need came also from the fact that we had old computers with limited A/D and poor storage units, the second because we needed the wind profiles immediately for certain applications in the air pollution monitoring. The SODAR is capable of producing a cumbersome amount of data even for today standards, especially if you want to store the raw data for advanced future analysis and because you have to digitize the signal continuously for many days, and sometimes for months. So that I had to reduce the rate and the amount of sampled data without losing the information we were interested on. The solution proceeded by successive approximations. My first approach was hardware and I realized, in 1980, an audio heterodyne that translated down the spectrum of the echo. It was also tried the decimation of the sampled data, comparing the spectra before and after and observing empirically that, in certain conditions, the result was only a down translation of the frequency bins without modification of the bins power. Then at the beginning of 1981 I ran into , p.230, and imagined that the “accordion-pleated” paper model had a useful mathematical formulation in terms of the fundamental formulae (1) and (2) for undersampling shown above. Only much more later I read  and  and realized that, at least (1) was already known, even if the topic was understated and treated differently (it was never named undersampling but “bandpass sampling theorem”) and partially and without proof in the quoted references, instead I think that my proof is simple and smart. In  the fundamental formula is stated differently and in the time domain instead of frequency domain, as I did. Furthermore no formula for is given. In  the “bandpass sampling theorem” is listed among the problems left to the reader and the formula shown refers only to the critical sampling frequency (3), but one of the terms may suggest, to an attentive reader, the way to compute . At that time, to my knowledge, people working on SODAR systems did not use the undersampling technique to digitize the signal, and even FFT was not so popular. Hence I think I was the first to introduce the undersampling in this area, and in a very simple form well suited for practical use. My fault was not to publicize enough my results with the consequence that a few people have tried to catch the merit for them even people I informed of personally. But, even if my report of 1983 was late, having I achieved the results in 1981 and even before, the papers of the others are all of two and more years later and, in a number of cases, I know why: at the beginning they did not believe in my results!
The emitted signal of the SODAR is a sinusoidal
burst typically of 100 milliseconds every 6 seconds. The receiver channels open
after the emission of the bursts. Basically, the signal received is strongly
modulated in amplitude at a low frequency, broaded in frequency by turbulence
and shifted in frequency because the Doppler effect, furthermore it is embedded
in a variable amount of environmental noise. With the SODAR it is possible to
visualize the atmospheric turbulence and measure the wind profile remotely up
Fig.6 Received power spectrum with very good S/N ratio
Our aim is to measure the vertical profiles of the speed and direction of the wind and the intensity of turbulence for each scan or for a number of averaged scans and at a number of height levels. Before digitizing we have to filter the signal with a bandpass filter with bandwidth , as shown in Fig.6, to avoid aliasing too much out of band noise. Furthermore, I lock the sampling frequency to the transmitted frequency to compensate for possible relative drifts that influence directly the precision of the wind measurement. Moreover, when the frequency bands of the three antennas are separated in such a manner to not interfere each other and to satisfy, in the whole, the requirements for undersampling, the three received signals are mixed together  before sampling, reducing at one third the required sampling rate. In a particular arrangement the whole bandwidth of the low frequency SODAR was so that the numbers are the same as in the previous example, except for the measuring units, which are uninfluential in this context. This is the moment when undersampling enters the scene: the ratio determines, in the example shown, the maximum factor of reduction on the number of sampled data which is a great saving in memory storage and a strong reduction (considering also the reduction in the sampling rate by mixing, when possible, of the three received signals) of the sampling rate that speed up the data analysis, permitting the real-time processing with cheaper instrumentation. After having undersampled the signal of each scan, we partition them in a number of segments and for each segment we apply the FFT to extract the spectrum shown in Fig.6.
In the backscattering mode, i.e. receiver coaxial or coincident with the transmitter, the axial component of the wind is given by
where is the sound velocity, the transmitted frequency and the received frequency. This is a first order approximation but it is pretty good in our atmosphere because . Typically the received frequency is defined as the first normalized moment of the spectrum
in which is the interval shown in Fig.6, is the power of the bin and is the resolution of the spectrum. In a typical triaxial arrangement the antennas are positioned as in Fig.7
Fig.7 Orientation of the antennas of a triaxial SODAR system
Having computed the three axial components , , of the wind with the formula (4), we need, for meteorological use, to transform them in cartesian components , , directed respectively along the axes of a orthogonal frame of reference. By definition is and where and are respectively the unit vectors of the antennas axes and of the cartesian axes. It is
in which we use the summation convention for the repeated indices in the products. It is hence with . Because we already know from (4), knowing also the geometric coefficients (cosine directors) we will be able to solve (5) to respect to the three unknowns . For the particular orientations depicted in Fig.7, we have
Solving to respect to the we finally obtain
Eventually, what began with undersampling SODAR
signals, has given its fruits. We are able to save, plot, print wind and
intensity profiles of the atmospheric boundary layer, up to about
Fig.8 Vertical profile of the intensity of the vertical backscattered signal versus time.
Fig.9 Vertical profile of the vertical wind component versus time.
 Abdul J. Jerry, “The Shannon Sampling Theorem – Its Various Extensions and Applications: A Tutorial Review”, Proc. of IEEE, vol.65, no. 11, November 1977.
 Julius S. Bendat, Allan G. Piersol, “Random Data: Analysis and Measurement Procedures”, Wiley-Interscience, 1971, (p.230).
 E. Oran Brigham, “The Fast Fourier Transform”, Prentice-Hall, Inc., 1974, (p.87).
 “Reference Data For Radio Engineers, Fifth Edition”, Howard W. Sams & Co., Inc., ITT, 1970, (p.21-14).
 Bonnie Baker, “Turning Nyquist upside down by undersampling”, EDN 12 May 2005.
 Angelo Ricotta, “Development of an
acoustic radar and its applications to the planetary boundary layer dynamics
studies”, Thesis in Physics,
 E. J. Owens, NOAA MARK VII Acoustic Echo
Sounder, NOAA Tech. Mem.,
 P. L. Butzer, J. R. Higgins, R. L. Stens, “Sampling theory of signal analysis”, Development of Mathematics 1950-2000, Editor J. P. Pier, Birkhäuser, 2000.
 A. Ricotta, M. Berico, “Triaxial Sodar”, Technical Report LPS 80-6, LPS-CNR, Frascati, 1980 (in Italian).