Nyquist–Shannon sampling theorem

From Wikipedia, the free encyclopedia

https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

The Nyquist–Shannon sampling theorem is an essential principle for digital signal processing linking the frequency range of a signal and the sample rate required to avoid a type of distortion called aliasing. The theorem states that the sample rate must be at least twice the bandwidth of the signal to avoid aliasing. In practice, it is used to select band-limiting filters to keep aliasing below an acceptable amount when an analog signal is sampled or when sample rates are changed within a digital signal processing function.

Example of magnitude of the Fourier transform of a bandlimited function

The Nyquist–Shannon sampling theorem is a theorem in the field of signal processing which serves as a fundamental bridge between continuous-time signals and discrete-time signals. It establishes a sufficient condition for a sample rate that permits a discrete sequence of samples to capture all the information from a continuous-time signal of finite bandwidth.

Strictly speaking, the theorem only applies to a class of mathematical functions having a Fourier transform that is zero outside of a finite region of frequencies. Intuitively we expect that when one reduces a continuous function to a discrete sequence and interpolates back to a continuous function, the fidelity of the result depends on the density (or sample rate) of the original samples. The sampling theorem introduces the concept of a sample rate that is sufficient for perfect fidelity for the class of functions that are band-limited to a given bandwidth, such that no actual information is lost in the sampling process. It expresses the sufficient sample rate in terms of the bandwidth for the class of functions. The theorem also leads to a formula for perfectly reconstructing the original continuous-time function from the samples.

Perfect reconstruction may still be possible when the sample-rate criterion is not satisfied, provided other constraints on the signal are known (see § Sampling of non-baseband signals below and compressed sensing). In some cases (when the sample-rate criterion is not satisfied), utilizing additional constraints allows for approximate reconstructions. The fidelity of these reconstructions can be verified and quantified utilizing Bochner's theorem.

The name Nyquist–Shannon sampling theorem honours Harry Nyquist and Claude Shannon, but the theorem was also previously discovered by E. T. Whittaker (published in 1915), and Shannon cited Whittaker's paper in his work. The theorem is thus also known by the names Whittaker–Shannon sampling theorem, Whittaker–Shannon, and Whittaker–Nyquist–Shannon, and may also be referred to as the cardinal theorem of interpolation.

Introduction

Sampling is a process of converting a signal (for example, a function of continuous time or space) into a sequence of values (a function of discrete time or space). Shannon's version of the theorem states:

Theorem — If a function ${\displaystyle x(t)}$ contains no frequencies higher than B hertz, then it can be completely determined from its ordinates at a sequence of points spaced less than ${\displaystyle 1/(2B)}$ seconds apart.

A sufficient sample-rate is therefore anything larger than ${\displaystyle 2B}$ samples per second. Equivalently, for a given sample rate ${\displaystyle f_s}$, perfect reconstruction is guaranteed possible for a bandlimit ${\displaystyle B<f_s/2}$.

When the bandlimit is too high (or there is no bandlimit), the reconstruction exhibits imperfections known as aliasing. Modern statements of the theorem are sometimes careful to explicitly state that ${\displaystyle x(t)}$ must contain no sinusoidal component at exactly frequency ${\displaystyle B,}$ or that ${\displaystyle B}$ must be strictly less than ½ the sample rate. The threshold ${\displaystyle 2B}$ is called the Nyquist rate and is an attribute of the continuous-time input ${\displaystyle x(t)}$ to be sampled. The sample rate must exceed the Nyquist rate for the samples to suffice to represent ${\displaystyle x(t).}$  The threshold ${\displaystyle f_s/2}$ is called the Nyquist frequency and is an attribute of the sampling equipment. All meaningful frequency components of the properly sampled ${\displaystyle x(t)}$ exist below the Nyquist frequency. The condition described by these inequalities is called the Nyquist criterion, or sometimes the Raabe condition. The theorem is also applicable to functions of other domains, such as space, in the case of a digitized image. The only change, in the case of other domains, is the units of measure attributed to ${\displaystyle t,}$ ${\displaystyle f_s,}$ and ${\displaystyle B.}$

The normalized sinc function: sin(πx) / (πx) ... showing the central peak at x = 0, and zero-crossings at the other integer values of x.

The symbol ${\displaystyle T\triangleq 1/f_s}$ is customarily used to represent the interval between samples and is called the sample period or sampling interval. The samples of function ${\displaystyle x(t)}$ are commonly denoted by ${\displaystyle x[n]\triangleq x(nT)}$ (alternatively ${\displaystyle x_n}$ in older signal processing literature), for all integer values of ${\displaystyle n.}$  Another convenient definition is ${\displaystyle x[n]\triangleq T\cdot x(nT),}$ which preserves the energy of the signal as ${\displaystyle T}$ varies.

A mathematically ideal way to interpolate the sequence involves the use of sinc functions. Each sample in the sequence is replaced by a sinc function, centered on the time axis at the original location of the sample ${\displaystyle nT,}$ with the amplitude of the sinc function scaled to the sample value, ${\displaystyle x[n].}$ Subsequently, the sinc functions are summed into a continuous function. A mathematically equivalent method uses the Dirac comb and proceeds by convolving one sinc function with a series of Dirac delta pulses, weighted by the sample values. Neither method is numerically practical. Instead, some type of approximation of the sinc functions, finite in length, is used. The imperfections attributable to the approximation are known as interpolation error.

Practical digital-to-analog converters produce neither scaled and delayed sinc functions, nor ideal Dirac pulses. Instead they produce a piecewise-constant sequence of scaled and delayed rectangular pulses (the zero-order hold), usually followed by a lowpass filter (called an "anti-imaging filter") to remove spurious high-frequency replicas (images) of the original baseband signal.

Aliasing

Main article: Aliasing

The samples of two sine waves can be identical when at least one of them is at a frequency above half the sample rate.

When ${\displaystyle x(t)}$ is a function with a Fourier transform ${\displaystyle X(f)}$:

${\displaystyle X(f)\triangleq {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}x(t)e^{-i2\pi ft}\mathrm{d}t,}$

the Poisson summation formula indicates that the samples, ${\displaystyle x(nT),}$ of ${\displaystyle x(t)}$ are sufficient to create a periodic summation of ${\displaystyle X(f).}$ The result is:

${\displaystyle X_s(f)\triangleq \sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}X\left(f-k{f}_s\right)=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}T\cdot x(nT)e^{-i2\pi nTf},}$

(Eq.1)

${\displaystyle X(f)}$ (top blue) and ${\displaystyle X_A(f)}$ (bottom blue) are continuous Fourier transforms of two different functions, ${\displaystyle x(t)}$ and ${\displaystyle x_A(t)}$ (not shown). When the functions are sampled at rate ${\displaystyle f_s}$, the images (green) are added to the original transforms (blue) when one examines the discrete-time Fourier transforms (DTFT) of the sequences. In this hypothetical example, the DTFTs are identical, which means the sampled sequences are identical, even though the original continuous pre-sampled functions are not. If these were audio signals, ${\displaystyle x(t)}$ and ${\displaystyle x_A(t)}$ might not sound the same. But their samples (taken at rate ${\displaystyle f_s}$) are identical and would lead to identical reproduced sounds; thus ${\displaystyle x_A(t)}$ is an alias of ${\displaystyle x(t)}$ at this sample rate.

which is a periodic function and its equivalent representation as a Fourier series, whose coefficients are ${\displaystyle T\cdot x(nT)}$. This function is also known as the discrete-time Fourier transform (DTFT) of the sample sequence.

As depicted, copies of ${\displaystyle X(f)}$ are shifted by multiples of the sampling rate ${\displaystyle f_s}$ and combined by addition. For a band-limited function ${\displaystyle (X(f)=0,\text{ for all }|f|\ge B)}$ and sufficiently large ${\displaystyle f_s,}$ it is possible for the copies to remain distinct from each other. But if the Nyquist criterion is not satisfied, adjacent copies overlap, and it is not possible in general to discern an unambiguous ${\displaystyle X(f).}$ Any frequency component above ${\displaystyle f_s/2}$ is indistinguishable from a lower-frequency component, called an alias, associated with one of the copies. In such cases, the customary interpolation techniques produce the alias, rather than the original component. When the sample-rate is pre-determined by other considerations (such as an industry standard), ${\displaystyle x(t)}$ is usually filtered to reduce its high frequencies to acceptable levels before it is sampled. The type of filter required is a lowpass filter, and in this application it is called an anti-aliasing filter.

Spectrum, ${\displaystyle X_s(f)}$, of a properly sampled bandlimited signal (blue) and the adjacent DTFT images (green) that do not overlap. A brick-wall low-pass filter, ${\displaystyle H(f)}$, removes the images, leaves the original spectrum, ${\displaystyle X(f)}$, and recovers the original signal from its samples.

The figure on the left shows a function (in gray/black) being sampled and reconstructed (in gold) at steadily increasing sample-densities, while the figure on the right shows the frequency spectrum of the gray/black function, which does not change. The highest frequency in the spectrum is half the width of the entire spectrum. The width of the steadily-increasing pink shading is equal to the sample-rate. When it encompasses the entire frequency spectrum it is twice as large as the highest frequency, and that is when the reconstructed waveform matches the sampled one.

Derivation as a special case of Poisson summation

When there is no overlap of the copies (also known as "images") of ${\displaystyle X(f)}$, the ${\displaystyle k=0}$ term of Eq.1 can be recovered by the product:

$${\displaystyle X(f)=H(f)\cdot X_s(f),}$$

where:

$${\displaystyle H(f)\triangleq \{\begin{array}{ll}1& |f|<B\\ 0& |f|>f_s-B.\end{array}}$$

The sampling theorem is proved since ${\displaystyle X(f)}$ uniquely determines ${\displaystyle x(t)}$.

All that remains is to derive the formula for reconstruction. ${\displaystyle H(f)}$ need not be precisely defined in the region ${\displaystyle [B,f_s-B]}$ because ${\displaystyle X_s(f)}$ is zero in that region. However, the worst case is when ${\displaystyle B=f_s/2,}$ the Nyquist frequency. A function that is sufficient for that and all less severe cases is:

$${\displaystyle H(f)=\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\left(\frac{f}{f_s}\right)=\{\begin{array}{ll}1& |f|<\frac{f_s}{2}\\ 0& |f|>\frac{f_s}{2},\end{array}}$$

where ${\displaystyle \mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}}$ is the rectangular function. Therefore:

${\displaystyle X(f)=\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\left(\frac{f}{f_s}\right)\cdot X_s(f)}$

${\displaystyle =\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}(Tf)\cdot \sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}T\cdot x(nT)e^{-i2\pi nTf}}$      (from  Eq.1, above).

${\displaystyle =\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}x(nT)\cdot \underset{\mathcal{F}\left\{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\frac{t-nT}{T}\right)\right\}}{\underbrace{T\cdot \mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}(Tf)\cdot e^{-i2\pi nTf}}}.}$     

The inverse transform of both sides produces the Whittaker–Shannon interpolation formula:

${\displaystyle x(t)=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}x(nT)\cdot \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\frac{t-nT}{T}\right),}$

which shows how the samples, ${\displaystyle x(nT)}$, can be combined to reconstruct ${\displaystyle x(t)}$.

• Larger-than-necessary values of ${\displaystyle f_s}$ (smaller values of ${\displaystyle T}$), called oversampling, have no effect on the outcome of the reconstruction and have the benefit of leaving room for a transition band in which ${\displaystyle H(f)}$ is free to take intermediate values. Undersampling, which causes aliasing, is not in general a reversible operation.
• Theoretically, the interpolation formula can be implemented as a low-pass filter, whose impulse response is ${\displaystyle \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}(t/T)}$ and whose input is ${\displaystyle \sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}x(nT)\cdot \delta (t-nT),}$ which is a Dirac comb function modulated by the signal samples. Practical digital-to-analog converters (DAC) implement an approximation like the zero-order hold. In that case, oversampling can reduce the approximation error.

Shannon's original proof

Poisson shows that the Fourier series in Eq.1 produces the periodic summation of ${\displaystyle X(f)}$, regardless of ${\displaystyle f_s}$ and ${\displaystyle B}$. Shannon, however, only derives the series coefficients for the case ${\displaystyle f_s=2B}$. Virtually quoting Shannon's original paper:

Let ${\displaystyle X(\omega )}$ be the spectrum of ${\displaystyle x(t).}$  Then

${\displaystyle x(t)=\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}X(\omega )e^{i\omega t}\mathrm{d}\omega =\frac{1}{2\pi }{\int }_{-2\pi B}^{2\pi B}X(\omega )e^{i\omega t}\mathrm{d}\omega ,}$

because ${\displaystyle X(\omega )}$ is assumed to be zero outside the band ${\displaystyle \left|\frac{\omega }{2\pi }\right|<B.}$  If we let ${\displaystyle t=\frac{n}{2B},}$ where ${\displaystyle n}$ is any positive or negative integer, we obtain:

${\displaystyle x\left(\frac{n}{2B}\right)=\frac{1}{2\pi }{\int }_{-2\pi B}^{2\pi B}X(\omega )e^{i\omega \frac{n}{2B}}\mathrm{d}\omega .}$

(Eq.2)

On the left are values of ${\displaystyle x(t)}$ at the sampling points. The integral on the right will be recognized as essentially the ${\displaystyle n^{th}}$ coefficient in a Fourier-series expansion of the function ${\displaystyle X(\omega ),}$ taking the interval ${\displaystyle -B}$ to ${\displaystyle B}$ as a fundamental period. This means that the values of the samples ${\displaystyle x(n/2B)}$ determine the Fourier coefficients in the series expansion of ${\displaystyle X(\omega ).}$  Thus they determine ${\displaystyle X(\omega ),}$ since ${\displaystyle X(\omega )}$ is zero for frequencies greater than ${\displaystyle B,}$ and for lower frequencies ${\displaystyle X(\omega )}$ is determined if its Fourier coefficients are determined. But ${\displaystyle X(\omega )}$ determines the original function ${\displaystyle x(t)}$ completely, since a function is determined if its spectrum is known. Therefore the original samples determine the function ${\displaystyle x(t)}$ completely.

Shannon's proof of the theorem is complete at that point, but he goes on to discuss reconstruction via sinc functions, what we now call the Whittaker–Shannon interpolation formula as discussed above. He does not derive or prove the properties of the sinc function, as the Fourier pair relationship between the rect (the rectangular function) and sinc functions was well known by that time.

Let ${\displaystyle x_n}$ be the ${\displaystyle n^{th}}$ sample. Then the function ${\displaystyle x(t)}$ is represented by:

${\displaystyle x(t)=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{x}_n\frac{\sin (\pi (2Bt-n))}{\pi (2Bt-n)}.}$

As in the other proof, the existence of the Fourier transform of the original signal is assumed, so the proof does not say whether the sampling theorem extends to bandlimited stationary random processes.

Notes

1. 
   

1. Multiplying both sides of Eq.2 by ${\displaystyle T=1/2B}$ produces, on the left, the scaled sample values ${\displaystyle (T\cdot x(nT))}$ in Poisson's formula (Eq.1), and, on the right, the actual formula for Fourier expansion coefficients.

Application to multivariable signals and images

Main article: Multidimensional sampling

Subsampled image showing a Moiré pattern

Properly sampled image

The sampling theorem is usually formulated for functions of a single variable. Consequently, the theorem is directly applicable to time-dependent signals and is normally formulated in that context. However, the sampling theorem can be extended in a straightforward way to functions of arbitrarily many variables. Grayscale images, for example, are often represented as two-dimensional arrays (or matrices) of real numbers representing the relative intensities of pixels (picture elements) located at the intersections of row and column sample locations. As a result, images require two independent variables, or indices, to specify each pixel uniquely—one for the row, and one for the column.

Color images typically consist of a composite of three separate grayscale images, one to represent each of the three primary colors—red, green, and blue, or RGB for short. Other colorspaces using 3-vectors for colors include HSV, CIELAB, XYZ, etc. Some colorspaces such as cyan, magenta, yellow, and black (CMYK) may represent color by four dimensions. All of these are treated as vector-valued functions over a two-dimensional sampled domain.

Similar to one-dimensional discrete-time signals, images can also suffer from aliasing if the sampling resolution, or pixel density, is inadequate. For example, a digital photograph of a striped shirt with high frequencies (in other words, the distance between the stripes is small), can cause aliasing of the shirt when it is sampled by the camera's image sensor. The aliasing appears as a moiré pattern. The "solution" to higher sampling in the spatial domain for this case would be to move closer to the shirt, use a higher resolution sensor, or to optically blur the image before acquiring it with the sensor using an optical low-pass filter.

Another example is shown here in the brick patterns. The top image shows the effects when the sampling theorem's condition is not satisfied. When software rescales an image (the same process that creates the thumbnail shown in the lower image) it, in effect, runs the image through a low-pass filter first and then downsamples the image to result in a smaller image that does not exhibit the moiré pattern. The top image is what happens when the image is downsampled without low-pass filtering: aliasing results.

The sampling theorem applies to camera systems, where the scene and lens constitute an analog spatial signal source, and the image sensor is a spatial sampling device. Each of these components is characterized by a modulation transfer function (MTF), representing the precise resolution (spatial bandwidth) available in that component. Effects of aliasing or blurring can occur when the lens MTF and sensor MTF are mismatched. When the optical image which is sampled by the sensor device contains higher spatial frequencies than the sensor, the under sampling acts as a low-pass filter to reduce or eliminate aliasing. When the area of the sampling spot (the size of the pixel sensor) is not large enough to provide sufficient spatial anti-aliasing, a separate anti-aliasing filter (optical low-pass filter) may be included in a camera system to reduce the MTF of the optical image. Instead of requiring an optical filter, the graphics processing unit of smartphone cameras performs digital signal processing to remove aliasing with a digital filter. Digital filters also apply sharpening to amplify the contrast from the lens at high spatial frequencies, which otherwise falls off rapidly at diffraction limits.

The sampling theorem also applies to post-processing digital images, such as to up or down sampling. Effects of aliasing, blurring, and sharpening may be adjusted with digital filtering implemented in software, which necessarily follows the theoretical principles.

A family of sinusoids at the critical frequency, all having the same sample sequences of alternating +1 and –1. That is, they all are aliases of each other, even though their frequency is not above half the sample rate.

Critical frequency

To illustrate the necessity of ${\displaystyle f_s>2B,}$ consider the family of sinusoids generated by different values of ${\displaystyle \theta }$ in this formula:

${\displaystyle x(t)=\frac{\cos (2\pi Bt+\theta )}{\cos (\theta )}=\cos (2\pi Bt)-\sin (2\pi Bt)\tan (\theta ),-\pi /2<\theta <\pi /2.}$

With ${\displaystyle f_s=2B}$ or equivalently ${\displaystyle T=1/2B,}$ the samples are given by:

${\displaystyle x(nT)=\cos (\pi n)-\underset{0}{\underbrace{\sin (\pi n)}}\tan (\theta )=(-1{)}^n}$

regardless of the value of ${\displaystyle \theta .}$ That sort of ambiguity is the reason for the strict inequality of the sampling theorem's condition.

Sampling of non-baseband signals

As discussed by Shannon:

A similar result is true if the band does not start at zero frequency but at some higher value, and can be proved by a linear translation (corresponding physically to single-sideband modulation) of the zero-frequency case. In this case the elementary pulse is obtained from ${\displaystyle \sin (x)/x}$ by single-side-band modulation.

That is, a sufficient no-loss condition for sampling signals that do not have baseband components exists that involves the width of the non-zero frequency interval as opposed to its highest frequency component. See sampling for more details and examples.

For example, in order to sample FM radio signals in the frequency range of 100–102 MHz, it is not necessary to sample at 204 MHz (twice the upper frequency), but rather it is sufficient to sample at 4 MHz (twice the width of the frequency interval).

A bandpass condition is that ${\displaystyle X(f)=0,}$ for all nonnegative ${\displaystyle f}$ outside the open band of frequencies:

${\displaystyle \left(\frac{N}{2}{f}_{\mathrm{s}},\frac{N+1}{2}{f}_{\mathrm{s}}\right),}$

for some nonnegative integer ${\displaystyle N}$. This formulation includes the normal baseband condition as the case ${\displaystyle N=0.}$

The corresponding interpolation function is the impulse response of an ideal brick-wall bandpass filter (as opposed to the ideal brick-wall lowpass filter used above) with cutoffs at the upper and lower edges of the specified band, which is the difference between a pair of lowpass impulse responses:

$${\displaystyle (N+1)\mathrm{sinc}\left(\frac{(N+1)t}{T}\right)-N\mathrm{sinc}\left(\frac{Nt}{T}\right).}$$

Other generalizations, for example to signals occupying multiple non-contiguous bands, are possible as well. Even the most generalized form of the sampling theorem does not have a provably true converse. That is, one cannot conclude that information is necessarily lost just because the conditions of the sampling theorem are not satisfied; from an engineering perspective, however, it is generally safe to assume that if the sampling theorem is not satisfied then information will most likely be lost.

Nonuniform sampling

The sampling theory of Shannon can be generalized for the case of nonuniform sampling, that is, samples not taken equally spaced in time. The Shannon sampling theory for non-uniform sampling states that a band-limited signal can be perfectly reconstructed from its samples if the average sampling rate satisfies the Nyquist condition. Therefore, although uniformly spaced samples may result in easier reconstruction algorithms, it is not a necessary condition for perfect reconstruction.

The general theory for non-baseband and nonuniform samples was developed in 1967 by Henry Landau. He proved that the average sampling rate (uniform or otherwise) must be twice the occupied bandwidth of the signal, assuming it is a priori known what portion of the spectrum was occupied.

In the late 1990s, this work was partially extended to cover signals for which the amount of occupied bandwidth is known but the actual occupied portion of the spectrum is unknown. In the 2000s, a complete theory was developed (see the section Sampling below the Nyquist rate under additional restrictions below) using compressed sensing. In particular, the theory, using signal processing language, is described in a 2009 paper by Mishali and Eldar. They show, among other things, that if the frequency locations are unknown, then it is necessary to sample at least at twice the Nyquist criteria; in other words, you must pay at least a factor of 2 for not knowing the location of the spectrum. Note that minimum sampling requirements do not necessarily guarantee stability.

Sampling below the Nyquist rate under additional restrictions

Main article: Undersampling

The Nyquist–Shannon sampling theorem provides a sufficient condition for the sampling and reconstruction of a band-limited signal. When reconstruction is done via the Whittaker–Shannon interpolation formula, the Nyquist criterion is also a necessary condition to avoid aliasing, in the sense that if samples are taken at a slower rate than twice the band limit, then there are some signals that will not be correctly reconstructed. However, if further restrictions are imposed on the signal, then the Nyquist criterion may no longer be a necessary condition.

A non-trivial example of exploiting extra assumptions about the signal is given by the recent field of compressed sensing, which allows for full reconstruction with a sub-Nyquist sampling rate. Specifically, this applies to signals that are sparse (or compressible) in some domain. As an example, compressed sensing deals with signals that may have a low overall bandwidth (say, the effective bandwidth ${\displaystyle EB}$) but the frequency locations are unknown, rather than all together in a single band, so that the passband technique does not apply. In other words, the frequency spectrum is sparse. Traditionally, the necessary sampling rate is thus ${\displaystyle 2B.}$ Using compressed sensing techniques, the signal could be perfectly reconstructed if it is sampled at a rate slightly lower than ${\displaystyle 2EB.}$ With this approach, reconstruction is no longer given by a formula, but instead by the solution to a linear optimization program.

Another example where sub-Nyquist sampling is optimal arises under the additional constraint that the samples are quantized in an optimal manner, as in a combined system of sampling and optimal lossy compression.^[9] This setting is relevant in cases where the joint effect of sampling and quantization is to be considered, and can provide a lower bound for the minimal reconstruction error that can be attained in sampling and quantizing a random signal. For stationary Gaussian random signals, this lower bound is usually attained at a sub-Nyquist sampling rate, indicating that sub-Nyquist sampling is optimal for this signal model under optimal quantization.

Historical background

The sampling theorem was implied by the work of Harry Nyquist in 1928, in which he showed that up to ${\displaystyle 2B}$ independent pulse samples could be sent through a system of bandwidth ${\displaystyle B}$; but he did not explicitly consider the problem of sampling and reconstruction of continuous signals. About the same time, Karl Küpfmüller showed a similar result and discussed the sinc-function impulse response of a band-limiting filter, via its integral, the step-response sine integral; this bandlimiting and reconstruction filter that is so central to the sampling theorem is sometimes referred to as a Küpfmüller filter (but seldom so in English).

The sampling theorem, essentially a dual of Nyquist's result, was proved by Claude E. Shannon. The mathematician E. T. Whittaker published similar results in 1915, J. M. Whittaker in 1935, and Gabor in 1946 ("Theory of communication").

In 1948 and 1949, Claude E. Shannon published the two revolutionary articles in which he founded the information theory. In Shannon 1948 the sampling theorem is formulated as "Theorem 13": Let ${\displaystyle f(t)}$ contain no frequencies over W. Then

$${\displaystyle f(t)=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{X}_n\frac{\sin \pi (2Wt-n)}{\pi (2Wt-n)},}$$

where ${\displaystyle X_n=f\left(\frac{n}{2W}\right).}$

It was not until these articles were published that the theorem known as "Shannon's sampling theorem" became common property among communication engineers, although Shannon himself writes that this is a fact which is common knowledge in the communication art. A few lines further on, however, he adds: "but in spite of its evident importance, [it] seems not to have appeared explicitly in the literature of communication theory".

Other discoverers

Others who have independently discovered or played roles in the development of the sampling theorem have been discussed in several historical articles, for example, by Jerri and by Lüke. For example, Lüke points out that H. Raabe, an assistant to Küpfmüller, proved the theorem in his 1939 Ph.D. dissertation; the term Raabe condition came to be associated with the criterion for unambiguous representation (sampling rate greater than twice the bandwidth). Meijering mentions several other discoverers and names in a paragraph and pair of footnotes:

As pointed out by Higgins, the sampling theorem should really be considered in two parts, as done above: the first stating the fact that a bandlimited function is completely determined by its samples, the second describing how to reconstruct the function using its samples. Both parts of the sampling theorem were given in a somewhat different form by J. M. Whittaker and before him also by Ogura. They were probably not aware of the fact that the first part of the theorem had been stated as early as 1897 by Borel. As we have seen, Borel also used around that time what became known as the cardinal series. However, he appears not to have made the link. In later years it became known that the sampling theorem had been presented before Shannon to the Russian communication community by Kotel'nikov. In more implicit, verbal form, it had also been described in the German literature by Raabe. Several authors have mentioned that Someya introduced the theorem in the Japanese literature parallel to Shannon. In the English literature, Weston introduced it independently of Shannon around the same time.

1. 
   

Several authors, following Black, have claimed that this first part of the sampling theorem was stated even earlier by Cauchy, in a paper published in 1841. However, the paper of Cauchy does not contain such a statement, as has been pointed out by Higgins.

2. As a consequence of the discovery of the several independent introductions of the sampling theorem, people started to refer to the theorem by including the names of the aforementioned authors, resulting in such catchphrases as "the Whittaker–Kotel'nikov–Shannon (WKS) sampling theorem" or even "the Whittaker–Kotel'nikov–Raabe–Shannon–Someya sampling theorem". To avoid confusion, perhaps the best thing to do is to refer to it as the sampling theorem, "rather than trying to find a title that does justice to all claimants".

— Eric Meijering, "A Chronology of Interpolation From Ancient Astronomy to Modern Signal and Image Processing" (citations omitted)

In Russian literature it is known as the Kotelnikov's theorem, named after Vladimir Kotelnikov, who discovered it in 1933.

Why Nyquist?

Exactly how, when, or why Harry Nyquist had his name attached to the sampling theorem remains obscure. The term Nyquist Sampling Theorem (capitalized thus) appeared as early as 1959 in a book from his former employer, Bell Labs, and appeared again in 1963, and not capitalized in 1965. It had been called the Shannon Sampling Theorem as early as 1954, but also just the sampling theorem by several other books in the early 1950s.

In 1958, Blackman and Tukey cited Nyquist's 1928 article as a reference for the sampling theorem of information theory, even though that article does not treat sampling and reconstruction of continuous signals as others did. Their glossary of terms includes these entries:

Sampling theorem (of information theory)

Nyquist's result that equi-spaced data, with two or more points per cycle of highest frequency, allows reconstruction of band-limited functions. (See Cardinal theorem.)

Cardinal theorem (of interpolation theory)

A precise statement of the conditions under which values given at a doubly infinite set of equally spaced points can be interpolated to yield a continuous band-limited function with the aid of the function

$${\displaystyle \frac{\sin (x-x_i)}{x-x_i}.}$$

Exactly what "Nyquist's result" they are referring to remains mysterious.

When Shannon stated and proved the sampling theorem in his 1949 article, according to Meijering, "he referred to the critical sampling interval ${\displaystyle T=\frac{1}{2W}}$ as the Nyquist interval corresponding to the band ${\displaystyle W,}$ in recognition of Nyquist's discovery of the fundamental importance of this interval in connection with telegraphy". This explains Nyquist's name on the critical interval, but not on the theorem.

Similarly, Nyquist's name was attached to Nyquist rate in 1953 by Harold S. Black:

If the essential frequency range is limited to ${\displaystyle B}$ cycles per second, ${\displaystyle 2B}$ was given by Nyquist as the maximum number of code elements per second that could be unambiguously resolved, assuming the peak interference is less half a quantum step. This rate is generally referred to as signaling at the Nyquist rate and ${\displaystyle \frac{1}{2B}}$ has been termed a Nyquist interval.

— Harold Black, Modulation Theory (bold added for emphasis; italics as in the original)

According to the Oxford English Dictionary, this may be the origin of the term Nyquist rate. In Black's usage, it is not a sampling rate, but a signaling rate.

Van der Pol oscillator

From Wikipedia, the free encyclopedia

https://en.wikipedia.org/wiki/Van_der_Pol_oscillator 

In the study of dynamical systems, the van der Pol oscillator (named for Dutch physicist Balthasar van der Pol) is a non-conservative, oscillating system with non-linear damping. It evolves in time according to the second-order differential equation

$${\displaystyle \frac{d^{2}x}{d{t}^2}-\mu (1-x^2)\frac{dx}{dt}+x=0,}$$

where x is the position coordinate—which is a function of the time t—and μ is a scalar parameter indicating the nonlinearity and the strength of the damping.

History

The Van der Pol oscillator was originally proposed by the Dutch electrical engineer and physicist Balthasar van der Pol while he was working at Philips. Van der Pol found stable oscillations, which he subsequently called relaxation-oscillations and are now known as a type of limit cycle, in electrical circuits employing vacuum tubes. When these circuits are driven near the limit cycle, they become entrained, i.e. the driving signal pulls the current along with it. Van der Pol and his colleague, van der Mark, reported in the September 1927 issue of Nature that at certain drive frequencies an irregular noise was heard, which was later found to be the result of deterministic chaos.

The Van der Pol equation has a long history of being used in both the physical and biological sciences. For instance, in biology, Fitzhugh and Nagumo extended the equation in a planar field as a model for action potentials of neurons. The equation has also been utilised in seismology to model the two plates in a geological fault, and in studies of phonation to model the right and left vocal fold oscillators.

Two-dimensional form

Liénard's theorem can be used to prove that the system has a limit cycle. Applying the Liénard transformation ${\displaystyle y=x-x^3/3-\dot{x}/\mu }$, where the dot indicates the time derivative, the Van der Pol oscillator can be written in its two-dimensional form:

${\displaystyle \dot{x}=\mu \left(x-\frac{1}{3}{x}^3-y\right)}$

${\displaystyle \dot{y}=\frac{1}{\mu }x}$.

Another commonly used form based on the transformation ${\displaystyle y=\dot{x}}$ leads to:

${\displaystyle \dot{x}=y}$

${\displaystyle \dot{y}=\mu (1-x^2)y-x}$.

Results for the unforced oscillator

Relaxation oscillation in the Van der Pol oscillator without external forcing. The nonlinear damping parameter is equal to μ = 5.

• When μ = 0, i.e. there is no damping function, the equation becomes
  $${\displaystyle \frac{d^{2}x}{d{t}^2}+x=0.}$$

  
  This is a form of the simple harmonic oscillator, and there is always conservation of energy.

• When μ > 0, all initial conditions converge to a globally unique limit cycle. Near the origin ${\displaystyle x=\frac{dx}{dt}=0,}$ the system is unstable, and far from the origin, the system is damped.
• The Van der Pol oscillator does not have an exact, analytic solution. However, such a solution does exist for the limit cycle if f(x) in the Lienard equation is a constant piece-wise function.
• The period at small μ has serial expansion
  $${\displaystyle T=\frac{2\pi }{1-{\mu }^2/16+17{\mu }^4/3072+O({\mu }^6)}.}$$

  
  See Poincaré–Lindstedt method for a derivation to order 2. See chapter 10 of  for a derivation up to order 3, and  for a numerical derivation up to order 164.
• For large μ, the behavior of the oscillator has a slow buildup, fast release cycle (a cycle of building up the tension and releasing the tension, thus a relaxation oscillation). This is most easily seen in the form
  $${\displaystyle \dot{x}=\mu \left(x-\frac{1}{3}{x}^3-y\right),\dot{y}=\frac{1}{\mu }x.}$$

  
  In this form, the oscillator completes one cycle as follows:
  • Slowly ascending the right branch of the cubic curve ${\displaystyle y=x-\frac{x^3}{3},}$ from (2, –2/3) to (1, 2/3).
  • Rapidly moving to the left branch of the cubic curve, from (1, 2/3) to (–2, 2/3).
  • Repeat the two steps on the left branch.
• The leading term in the period of the cycle is due to the slow ascending and descending, which can be computed as:
  $${\displaystyle T=2\int dt=2\int \mu \frac{dy}{x}=2\mu {\int }_2^1\frac{dy}{dx}\frac{dx}{x}=(3-2\ln 2)\mu }$$

  
  Higher orders of the period of the cycle is
  $${\displaystyle T=(3-2\ln 2)\mu +3\alpha {\mu }^{-1/3}-\frac{23}{{\mu }^{-1}}\ln \mu +O(\ln {\mu }^{-1}).}$$

  
  where α ≈ 2.338 is the smallest root of Ai(–α) = 0, where Ai is the Airy function.(Section 9.7) (contains a derivation, but has a misprint of 3α to 2α.)
• The amplitude of the cycle is 
  $${\displaystyle 2+\frac{\alpha }{3}{\mu }^{-4/3}-\frac{16}{27}{\mu }^{-2}\ln \mu +O({\mu }^{-2})}$$

  

Hopf bifurcation

As μ moves from less than zero to more than zero, the spiral sink at origin becomes a spiral source, and a limit cycle appears "out of the blue" with radius two. This is because the transition is not generic: when ε = 0, both the differential equation becomes linear, and the origin becomes a circular node.

Knowing that in a Hopf bifurcation, the limit cycle should have size ${\displaystyle \propto {\epsilon }^{1/2},}$ we may attempt to convert this to a Hopf bifurcation by using the change of variables ${\displaystyle u={\epsilon }^{1/2}x,}$ which gives

$${\displaystyle \ddot{u}+u+u^2\dot{u}-\epsilon \dot{u}=0}$$

This indeed is a Hopf bifurcation.

Hamiltonian for Van der Pol oscillator

Randomly chosen initial conditions are attracted to a stable orbit.

One can also write a time-independent Hamiltonian formalism for the Van der Pol oscillator by augmenting it to a four-dimensional autonomous dynamical system using an auxiliary second-order nonlinear differential equation as follows:

${\displaystyle \ddot{x}-\mu (1-x^2)\dot{x}+x=0,}$

${\displaystyle \ddot{y}+\mu (1-x^2)\dot{y}+y=0.}$

Note that the dynamics of the original Van der Pol oscillator is not affected due to the one-way coupling between the time-evolutions of x and y variables. A Hamiltonian H for this system of equations can be shown to be

${\displaystyle H(x,y,p_x,p_y)=p_x{p}_y+xy-\mu (1-x^2)y{p}_y,}$

where ${\displaystyle p_x=\dot{y}+\mu (1-x^2)y}$ and ${\displaystyle p_y=\dot{x}}$ are the conjugate momenta corresponding to x and y, respectively. This may, in principle, lead to quantization of the Van der Pol oscillator. Such a Hamiltonian also connects the geometric phase of the limit cycle system having time dependent parameters with the Hannay angle of the corresponding Hamiltonian system.

Quantum oscillator

The quantum van der Pol oscillator, which is the quantum mechanics version of the classical van der Pol oscillator, has been proposed using a Lindblad equation to study its quantum dynamics and quantum synchronization. Note the above Hamiltonian approach with an auxiliary second-order equation produces unbounded phase-space trajectories and hence cannot be used to quantize the van der Pol oscillator. In the limit of weak nonlinearity (i.e. μ→0) the van der Pol oscillator reduces to the Stuart–Landau equation. The Stuart–Landau equation in fact describes an entire class of limit-cycle oscillators in the weakly-nonlinear limit. The form of the classical Stuart–Landau equation is much simpler, and perhaps not surprisingly, can be quantized by a Lindblad equation which is also simpler than the Lindblad equation for the van der Pol oscillator. The quantum Stuart–Landau model has played an important role in the study of quantum synchronisation (where it has often been called a van der Pol oscillator although it cannot be uniquely associated with the van der Pol oscillator). The relationship between the classical Stuart–Landau model (μ→0) and more general limit-cycle oscillators (arbitrary μ) has also been demonstrated numerically in the corresponding quantum models.

Forced Van der Pol oscillator

Chaotic behaviour in the Van der Pol oscillator with sinusoidal forcing. The nonlinear damping parameter is equal to μ = 8.53, while the forcing has amplitude A = 1.2 and angular frequency ω = 2π/10.

The forced, or driven, Van der Pol oscillator takes the 'original' function and adds a driving function Asin(ωt) to give a differential equation of the form:

${\displaystyle \frac{d^{2}x}{d{t}^2}-\mu (1-x^2)\frac{dx}{dt}+x-A\sin (\omega t)=0,}$

where A is the amplitude, or displacement, of the wave function and ω is its angular velocity.

Popular culture

Electrical circuit involving a triode, resulting in a forced Van der Pol oscillator. The circuit contains: a triode, a resistor R, a capacitor C, a coupled inductor-set with self inductance L and mutual inductance M. In the serial RLC circuit there is a current i, and towards the triode anode ("plate") a current i_a, while there is a voltage u_g on the triode control grid. The Van der Pol oscillator is forced by an AC voltage source E_s.

Author James Gleick described a vacuum tube Van der Pol oscillator in his book from 1987 Chaos: Making a New Science. According to a New York Times article, Gleick received a modern electronic Van der Pol oscillator from a reader in 1988.

APMonitor

From Wikipedia, the free encyclopedia

https://en.wikipedia.org/wiki/APMonitor 

 

APMonitor

Developer(s)	APMonitor
Stable release	v1.0.1 / January 31, 2022
Repository	https://github.com/APMonitor/
Operating system	Cross-platform
Type	Technical computing
License	Proprietary, BSD
Website	APMonitor product page

Advanced process monitor (APMonitor) is a modeling language for differential algebraic (DAE) equations. It is a free web-service or local server for solving representations of physical systems in the form of implicit DAE models. APMonitor is suited for large-scale problems and solves linear programming, integer programming, nonlinear programming, nonlinear mixed integer programming, dynamic simulation, moving horizon estimation, and nonlinear model predictive control. APMonitor does not solve the problems directly, but calls nonlinear programming solvers such as APOPT, BPOPT, IPOPT, MINOS, and SNOPT. The APMonitor API provides exact first and second derivatives of continuous functions to the solvers through automatic differentiation and in sparse matrix form.

Programming language integration

Julia, MATLAB, Python are mathematical programming languages that have APMonitor integration through web-service APIs. The GEKKO Optimization Suite is a recent extension of APMonitor with complete Python integration. The interfaces are built-in optimization toolboxes or modules to both load and process solutions of optimization problems. APMonitor is an object-oriented modeling language and optimization suite that relies on programming languages to load, run, and retrieve solutions. APMonitor models and data are compiled at run-time and translated into objects that are solved by an optimization engine such as APOPT or IPOPT. The optimization engine is not specified by APMonitor, allowing several different optimization engines to be switched out. The simulation or optimization mode is also configurable to reconfigure the model for dynamic simulation, nonlinear model predictive control, moving horizon estimation or general problems in mathematical optimization.

As a first step in solving the problem, a mathematical model is expressed in terms of variables and equations such as the Hock & Schittkowski Benchmark Problem #71used to test the performance of nonlinear programming solvers. This particular optimization problem has an objective function ${\displaystyle \underset{x\in \mathbb{R}}{min}{x}_1{x}_4(x_1+x_2+x_3)+x_3}$ and subject to the inequality constraint ${\displaystyle x_1{x}_2{x}_3{x}_4\ge 25}$ and equality constraint ${\displaystyle {x_1}^2+{x_2}^2+{x_3}^2+{x_4}^2=40}$. The four variables must be between a lower bound of 1 and an upper bound of 5. The initial guess values are ${\displaystyle x_1=1,x_2=5,x_3=5,x_4=1}$. This mathematical model is translated into the APMonitor modeling language in the following text file.

! file saved as hs71.apm
Variables
  x1 = 1, >=1, <=5
  x2 = 5, >=1, <=5
  x3 = 5, >=1, <=5
  x4 = 1, >=1, <=5
End Variables

Equations
  minimize x1*x4*(x1+x2+x3) + x3

  x1*x2*x3*x4 > 25
  x1^2 + x2^2 + x3^2 + x4^2 = 40
End Equations

The problem is then solved in Python by first installing the APMonitor package with pip install APMonitor or from the following Python code.

# Install APMonitor
import pip
pip.main(['install','APMonitor'])

Installing a Python is only required once for any module. Once the APMonitor package is installed, it is imported and the apm_solve function solves the optimization problem. The solution is returned to the programming language for further processing and analysis.

# Python example for solving an optimization problem
from APMonitor.apm import *

# Solve optimization problem
sol = apm_solve("hs71", 3)

# Access solution
x1 = sol["x1"]
x2 = sol["x2"]

Similar interfaces are available for MATLAB and Julia with minor differences from the above syntax. Extending the capability of a modeling language is important because significant pre- or post-processing of data or solutions is often required when solving complex optimization, dynamic simulation, estimation, or control problems.

High Index DAEs

The highest order of a derivative that is necessary to return a DAE to ODE form is called the differentiation index. A standard way for dealing with high-index DAEs is to differentiate the equations to put them in index-1 DAE or ODE form (see Pantelides algorithm). However, this approach can cause a number of undesirable numerical issues such as instability. While the syntax is similar to other modeling languages such as gProms, APMonitor solves DAEs of any index without rearrangement or differentiation. As an example, an index-3 DAE is shown below for the pendulum motion equations and lower index rearrangements can return this system of equations to ODE form (see Index 0 to 3 Pendulum example).

Pendulum motion (index-3 DAE form)

Model pendulum
  Parameters
    m = 1
    g = 9.81
    s = 1
  End Parameters

  Variables
    x = 0
    y = -s
    v = 1
    w = 0
    lam = m*(1+s*g)/2*s^2
  End Variables

  Equations
    x^2 + y^2 = s^2
    $x = v
    $y = w
    m*$v = -2*x*lam
    m*$w = -m*g - 2*y*lam
  End Equations
End Model

Applications in APMonitor Modeling Language

Many physical systems are naturally expressed by differential algebraic equation. Some of these include:

• cell cultures
• chemical reactors
• cogeneration (power and heat)
• distillation columns
• drilling automation
• essential oil steam distillation
• friction stir welding
• hydrate formation in deep-sea pipelines
• infectious disease spread
• oscillators
• severe slugging control
• solar thermal energy production
• solid oxide fuel cells
• space shuttle launch simulation
• Unmanned Aerial Vehicles (UAVs)

Models for a direct current (DC) motor and blood glucose response of an insulin dependent patient are listed below. They are representative of differential and algebraic equations encountered in many branches of science and engineering.

Direct current (DC) motor

Parameters
  ! motor parameters (dc motor)
  v   = 36        ! input voltage to the motor (volts)
  rm  = 0.1       ! motor resistance (ohms)
  lm  = 0.01      ! motor inductance (henrys)
  kb  = 6.5e-4    ! back emf constant (volt·s/rad)
  kt  = 0.1       ! torque constant (N·m/a)
  jm  = 1.0e-4    ! rotor inertia (kg m²)
  bm  = 1.0e-5    ! mechanical damping (linear model of friction: bm * dth)

  ! load parameters
  jl = 1000*jm    ! load inertia (1000 times the rotor)
  bl = 1.0e-3     ! load damping (friction)
  k = 1.0e2       ! spring constant for motor shaft to load
  b = 0.1         ! spring damping for motor shaft to load
End Parameters

Variables
  i     = 0       ! motor electric current (amperes)
  dth_m = 0       ! rotor angular velocity sometimes called omega (radians/sec)
  th_m  = 0       ! rotor angle, theta (radians)
  dth_l = 0       ! wheel angular velocity (rad/s)
  th_l  = 0       ! wheel angle (radians)
End Variables

Equations
  lm*$i - v = -rm*i -    kb *$th_m
  jm*$dth_m =  kt*i - (bm+b)*$th_m - k*th_m +     b *$th_l + k*th_l
  jl*$dth_l =             b *$th_m + k*th_m - (b+bl)*$th_l - k*th_l
  dth_m = $th_m
  dth_l = $th_l 
End Equations

Blood glucose response of an insulin dependent patient

! Model source:
! A. Roy and R.S. Parker. “Dynamic Modeling of Free Fatty 
!   Acids, Glucose, and Insulin: An Extended Minimal Model,”
!   Diabetes Technology and Therapeutics 8(6), 617-626, 2006.
Parameters
  p1 = 0.068       ! 1/min
  p2 = 0.037       ! 1/min
  p3 = 0.000012    ! 1/min
  p4 = 1.3         ! mL/(min·µU)
  p5 = 0.000568    ! 1/mL
  p6 = 0.00006     ! 1/(min·µmol)
  p7 = 0.03        ! 1/min
  p8 = 4.5         ! mL/(min·µU)
  k1 = 0.02        ! 1/min
  k2 = 0.03        ! 1/min
  pF2 = 0.17       ! 1/min
  pF3 = 0.00001    ! 1/min
  n = 0.142        ! 1/min
  VolG = 117       ! dL
  VolF = 11.7      ! L
  ! basal parameters for Type-I diabetic
  Ib = 0           ! Insulin (µU/mL)
  Xb = 0           ! Remote insulin (µU/mL)
  Gb = 98          ! Blood Glucose (mg/dL)
  Yb = 0           ! Insulin for Lipogenesis (µU/mL)
  Fb = 380         ! Plasma Free Fatty Acid (µmol/L)
  Zb = 380         ! Remote Free Fatty Acid (µmol/L)
  ! insulin infusion rate
  u1 = 3           ! µU/min
  ! glucose uptake rate
  u2 = 300         ! mg/min
  ! external lipid infusion
  u3 = 0           ! mg/min
End parameters

Intermediates
  p9 = 0.00021 * exp(-0.0055*G)  ! dL/(min*mg)
End Intermediates

Variables
  I = Ib
  X = Xb
  G = Gb
  Y = Yb
  F = Fb
  Z = Zb
End variables

Equations
  ! Insulin dynamics
  $I = -n*I  + p5*u1
  ! Remote insulin compartment dynamics
  $X = -p2*X + p3*I
  ! Glucose dynamics
  $G = -p1*G - p4*X*G + p6*G*Z + p1*Gb - p6*Gb*Zb + u2/VolG
  ! Insulin dynamics for lipogenesis
  $Y = -pF2*Y + pF3*I
  ! Plasma-free fatty acid (FFA) dynamics
  $F = -p7*(F-Fb) - p8*Y*F + p9 * (F*G-Fb*Gb) + u3/VolF
  ! Remote FFA dynamics
  $Z = -k2*(Z-Zb) + k1*(F-Fb)
End Equations