1 Introduction
In practical applications, one often has a number of data points, obtained by sampling or experimentation, which represent the values of a continuous signal. Interpolation by simple functions (e.g., trigonometric functions or rational functions), as one of the important mathematical techniques used in physics and engineering sciences, deals with the problem of reconstructing or approximating the continuous signals from a series of discrete points. An ideal interpolation usually tells us that a continuous signal can be exactly and uniquely reconstructed from the discrete points if it satisfies some suitable conditions (e.g., bandlimited). Such an interpolation is also referred to as a sampling theorem in signal processing [1]
. In general, sampling formulas interpolate the given data even if the specific conditions for perfect reconstruction cannot be met. The error analysis for interpolation or sampling formulas is of great importance because the continuous signals that need to be recovered do not satisfy the conditions for ideal interpolation in most circumstances. Due to the wide range applications of interpolation, finding new interpolation or sampling formulas with error estimations as well as developing their fast algorithms for implementation have received considerable attentions in recently years (see e.g.,
[2, 3, 4, 5, 6, 7]).Most of classical sampling formulas have centered around reconstructing a signal from its own samples [1, 6]. In fact, reconstructing a signal from data other than the samples of the original signal is possible. Papoulis [8] first proposed generalized sampling expansion (GSE) for bandlimited signals defined on real line . The GSE indicates that a bandlimited signal can be reconstructed from the samples of output signals of linear systems. That is, one can reconstruct from the samples of the output signals
where is Fourier transform of and are the system functions and . The idea of GSE has been extended in different directions. Cheung [9] introduced GSE for realvalued multidimensional signals associated with Fourier transform, while Wei, Ran and Li [10, 11] presented the GSE with generalized integral transformation, such as fractional Fourier transform and linear canonical transform. In [12], the authors studied GSE for quaternionvalued signals associated with quaternion Fourier transform. Importantly, the applications based on GSE as well as its generalizations have been conducted by researchers [3, 13].
A finite duration signal is commonly assumed to be periodic in practice. Accordingly, there are certain studies concerning the interpolation theory of periodic signals. A sampling theorem for periodic signals was first introduced by Goldman [14]. The author in [15] proposed pseudoorthogonal bases for reconstructing periodic signals. In a series of papers [16, 17, 18], researchers extensively discussed the sinc interpolation of discrete periodic signals and they also derived several equivalent interpolation formulas in distinct forms. As an extension of periodic sinc interpolation, decomposing a periodic signal in a basis of shifted and scaled versions of a generating function was studied in [19]. Moreover, they further presented an error analysis for their approximation method. Recently, the nonuniform sampling theorems for periodic signals were also presented [20, 21].
In this paper, a multichannel interpolation for periodic signals is studied. We derive a general interpolation formula that depends on the transfer functions of the linear systems. The formula bears a resemblance to the classical GSE defined on real line. Nevertheless, the recipe of derivation is different from the traditional one, and moreover, the proposed formula is given by a finite summation, as opposed to a infinite summation in the traditional case. The contributions of this paper are summarized as follows:

We propose a multichannel interpolation for periodic signals (MIP). The proposed MIP is capable of generating various useful interpolation formulas by selecting suitable parameters according to the types and the amount of collected data. In addition, it would restore the original signal and some integral transformations (such as Hilbert transform) of .

A fast algorithm based on FFT is presented which brings high computational efficiency and reliability of the proposed algorithm of MIP.

To the authors’ knowledge, there seems no error analysis for classical GSE in approximating general finite energy signal on real line in the literature. By contrast, we analyze the errors that arise in approximating nonbandlimited periodic signals by the proposed MIP.

The proposed MIP is applied to image scaling problem. Its main advantage makes good use of multifaceted information of image, so that the scaled image retains lots of information of the original unscaled image. The experimental simulations show the superiority of our method over the other three popularly used conventional methods: nearest neighbour, bilinear and bicubic.
The rest of the paper is organized as follows: Section 2 recalls some preliminaries of Fourier series. Section 3 formulates MIP and presents some examples to illustrate how to use MIP flexibly. The error analysis is drawn in Section 4. In Section 5, the effectiveness of the proposed MIP for approximating signals is demonstrated by several numerical examples and the application of MIP to image scaling is also addressed. Finally, conclusions are made in Section 6.
2 Preliminaries
This part recalls some preparatory knowledge of Fourier series (see e.g. [22]). Without loss of generality, we restrict attention to the periodic signals, or alternatively, we will consider the signals defined on unit circle .
Let be the totality of functions such that
For , it can be written as
(1) 
with , where the Fourier series is convergent to in norm. It is well known that is a Hilbert space with the inner product
If and their Fourier coefficients are respectively given by
By Hölder inequality, we have and the general version of Parseval’s identity is of the form:
Moreover, the convolution theorem gives
(2) 
The analytic signal defined by the linear combination of original signal and its (circular) Hilbert transform is regarded as an useful representation from which the phase, energy and frequency may be estimated [23, 24]. The circular Hilbert transform [25, 26] for periodic signal is given by
where sgn is the signum function taking values , or for , or respectively. By definition, we have
(3) 
We will first concentrate on reconstruction problem of finite order trigonometric polynomials. To maintain consistent terminology with the classical case, in what follows, a finite order trigonometric polynomial is called a bandlimited signal. Specifically, a periodic signal is said to be bandlimited if its sequence of Fourier coefficients possesses finite nonzero elements. Let , in the sequel we denote by the totality of bandlimited signals with the following form:
(4) 
3 Formulation of MIP
Throughout the paper, let , and assume that is divisible by , namely, . We cut the set of integers into pieces for convenience. Let us set
Then we have and .
For , let
(5)  
It follows from (2) that
where . We particularly mention that the series (5) may not be convergent in general. Nevertheless, is well defined when . In this case, may be regarded as a distribution.
The proposed MIP is to reconstruct from the samples of . To achieve this, it is naturally to expect that the simultaneous sampling of signals will reduces sampling rate by . we shall note that the signal expressed as Eq. (4
) may be represented by a shorter length of summation as long as we introduce the following vectors.
For , let
It follows that
(6) 
Similar considerations applying to , we have
where
Owing to the periodicity of , we easily obtain
for every . This leads to a simple expression for samples of , that is
(7) 
where . The factor gives us a hint that we may be able to express the samples of in terms of the DFT matrix or its inverse matrix.
Lemma 1
Let and There is a matrix representation for samples of in terms of the inverse DFT matrix. That is,
(8) 
where is the th order DFT matrix
(9) 
with and is a diagonal matrix
(10) 
Proof. From Eq. (7), we can rewrite as
Note that , and , it follows that
Hence we immediately obtain Eq. (8) which completes the proof.
Following the definition of , for , we further set by matrix
(11) 
which means the element of is . Suppose that is invertible for every and we denote the inverse matrix as
Next we use the elements of to construct the interpolation functions. Let
(12) 
and define
(13) 
As with most conventional bandlimited interpolation methods, the interpolation functions of MIP are generated from some fixed functions by translations. The following result provides a matrix representation of shifted functions generated by .
Lemma 2
Proof. By the definition of and rearranging the summation terms of Eq. (13), we get
From the piecewise defined in Eq. (12), it is easy to verify that for all and . It follows that can be rewritten as
(15) 
In view of the resemblance of Eq. (15) and Eq. (7) and by similar arguments to Lemma 1, we immediately have Eq. (14) which completes the proof.
Until now, in addition to representing the sampled values of by , the interpolation functions (, ) have been constructed from the elements of . we shall now make a connection between and . To achieve this, let
By straightforward computations, we get
(16)  
(17) 
Having introduced above notations, we are in position to show our main result which is referred to multichannel interpolation for periodic signals (MIP).
Theorem 1
Proof. Multiplying Eq. (17) by and using Eq. (16), the expression (4) can be rewritten as
(19)  
The last equality is a consequence of Lemma 1 and 2. The proof is complete.
Remark 1
In general, has to be invertible for every because the definition of depends on () which are the elements of . However, if there is an index set such that for every , then is not necessary to be invertible for . This is due to the fact that the zero terms of (6) may be removed from the summation. In this case, we just need to set () for .
Now we use some examples to show how Theorem 1 can be flexibly used to derive various sampling formulas for periodic signals. For simplicity we shall restrict to the case where .
Example 1
We first present the most basic example for our results. Let , , . By Theorem 1, we conclude that can be recovered from its samples. That is,
(20) 
where
is the th order Dirichlet kernel. This formula is referred to trigonometric interpolation for data (). It should be stressed that the interpolation formulas given in [16, 17, 18] are mathematically equivalent to Eq. (20). Note that the expression (20) is not numerically stable at which the function takes form. Based on the DFT, the author in [17] obtained a numerically stable formulation which is somehow equivalent to interpolation by the FFT [27]. In fact, the DFT or FFT based expressions for trigonometric interpolation can be subsumed in Eq. (19).
Example 2
Theorem 1 permits us to express from its samples and samples of its derivatives. Note that
Let , then should be defined by
By straightforward computations,
In order to satisfy the assumption that is divisible by , we have to assume that . Suppose that , then . Consequently, we obtain as follows:
By using and letting in Eq. (15), we have that equals
and
Note that for (), the determinant of defined by (11) is a Vandermonde determinant. Therefore is invertible for all . This would imply that we can further derive a sampling formula for recovering from its samples as well as the samples of its first derivatives.
Having considered interpolation for the samples of original signal and its derivatives. We now continue to show how the samples of Hilbert transform of can be used to reconstruct itself. Interestingly, based on antiinvolution property (3), MIP can be applied to compute Hilbert transform as well.
Example 3
Remark 2
Given a realvalued point discrete signal , . By replacing with and letting in the right hand side of (22), we get a discrete signal . The DFT for is
where is the DFT of . That means that is the discrete Hilbert transform of [28]. This fact reveals the relationship between discrete Hilbert transform and continuous circular Hilbert transform. For the discrete signal of even number, similar arguments can be made, we omit the details.
Example 4
The analytic signal associated with is defined by , it can be rewritten as
where for and for . Applying Theorem 1 to , we have the following. If , then
where and with
and
Moreover, if is realvalued, then
(24)  
(25) 
Example 5
In Example 4, we express (or ) by the samples of and . However, the formulas (24) and (25) are only valid for the realvalued signals. By Theorem 1, it is natural to expect that we can deduce a new expression of (resp. ) in terms of the samples of and by letting . Rewrite as with . Let and . Then , and
It is clear that