DeepAI

# Multichannel Interpolation for Periodic Signals via FFT, Error Analysis and Image Scaling

This paper describes a new method for the multichannel interpolation of a discrete signal. It is shown that a bandlimited periodic signal f can be exactly reconstructed from finite samples of g_k (1≤ k≤ M) which are the responses of M linear systems with input f. The proposed interpolation can also be applied to approximate non-bandlimited periodic signals. Quantitative error analysis is provided to ensure its effectiveness for approximating non-bandlimited periodic signals and its Hilbert transform. Importantly, involving the fast Fourier transform (FFT) computational technique, bring high computational efficiency and reliability of the proposed algorithm. The superior performance of the proposed algorithm is demonstrated by several simulations. Additionally, the proposed interpolation is applied to the image scaling problem. The empirical studies show that the proposed method performs better than the other conventional image scaling methods.

• 3 publications
• 12 publications
04/29/2019

### Experimental Demonstration of Data Transmission Based on the Exact Inverse Periodic Nonlinear Fourier Transform

We design a two-dimensional signal constellation based on the exact peri...
07/13/2020

### Fast approximation by periodic kernel-based lattice-point interpolation with application in uncertainty quantification

This paper deals with the kernel-based approximation of a multivariate p...
01/29/2021

### A Fast Template Periodogram for Detecting Non-sinusoidal Fixed-shape Signals in Irregularly Sampled Time Series

Astrophysical time series often contain periodic signals. The large and ...
11/11/2021

### Processing of large sets of stochastic signals: filtering based on piecewise interpolation technique

Suppose K__Y and K__X are large sets of observed and reference signals, ...
01/30/2023

### Time Encoding Sampling of Bandpass Signals

This paper investigates the problem of sampling and reconstructing bandp...
04/22/2022

### Speaking-Rate-Controllable HiFi-GAN Using Feature Interpolation

This paper presents a speaking-rate-controllable HiFi-GAN neural vocoder...
06/06/2019

### ZeLiC and ZeChipC: Time Series Interpolation Methods for Lebesgue or Event-based Sampling

Lebesgue sampling is based on collecting information depending on the va...

## 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

 gm(t)=1√2π∫σ−σF(ω)Hm(ω)eitωdω,  m=1,…,M

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 real-valued 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 quaternion-valued 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 pseudo-orthogonal 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 non-uniform 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:

1. 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 .

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

3. 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 non-bandlimited periodic signals by the proposed MIP.

4. 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

 ∥f∥p:=(12π∫T|f(t)|pdt)1p<∞.

For , it can be written as

 f(t)=∑n∈Za(n)eint (1)

with , where the Fourier series is convergent to in norm. It is well known that is a Hilbert space with the inner product

 (f,h):=12π∫Tf(t)¯¯¯¯¯¯¯¯¯h(t)dt,   ∀f,h∈L2(T).

If and their Fourier coefficients are respectively given by

 f∼a(n),    h∼b(n).

By Hölder inequality, we have and the general version of Parseval’s identity is of the form:

 (f,h)=∑n∈Za(n)¯¯¯¯¯¯¯¯¯¯b(n).

Moreover, the convolution theorem gives

 (f∗h)(t):=12π∫Tf(s)h(t−s)ds=∑n∈Za(n)b(n)eint. (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

 Hf(t)=∑n∈Z(−isgn(n))a(n)eint=12πp.v.∫Tf(s)cot(t−s2)ds,

where sgn is the signum function taking values , or for , or respectively. By definition, we have

 H2f(t)=−f(t)+a(0). (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:

 f(t)=∑n∈INa(n)eint,   IN={n:N1≤n≤N2}. (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

 Ik={n:N1+(k−1)L≤n≤N1+kL−1},  Jk=M+k⋃l=k+1Il.

Then we have and .

For , let

 hm(t) =∑n∈Zbm(n)eint, (5) gm(t) =(f∗hm)(t)=12π∫Tf(s)hm(t−s)ds.

It follows from (2) that

 gm(t)=∑n∈Zcm(n)eint,

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

 An=[a(n),a(n+L),a(n+2L),⋯,a(n+(M−1)L)], En(t)=[eint,ei(n+L)t,ei(n+2L)t,⋯,ei(n+(M−1)L)t]T.

It follows that

 f(t)=∑n∈INa(n)eint=∑n=I1M−1∑k=0a(n+kL)ei(n+kL)t=∑n=I1AnEn(t). (6)

Similar considerations applying to , we have

 gm(t)=∑n∈I1Cm,nEn(t),

where

 Cm,n=[c(n),c(n+L),c(n+2L),⋯,c(n+(M−1)L)].

Owing to the periodicity of , we easily obtain

 En(2πpL)=ein2πpL[1,1,⋯,1]

for every . This leads to a simple expression for samples of , that is

 gm(2πpL)=∑n∈I1Cm,nEn(2πpL)=∑n∈I1ein2πpLM−1∑k=0cm(n+kL)=∑n∈I1dm(n)ein2πpL, (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,

 1L[gm(t0),gm(t1),⋯,gm(tL−1)]=˜DmF−1LU−1L, (8)

where is the -th order DFT matrix

 FL=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ω0ω0ω0⋯ω0ω0ω1ω2⋯ωL−1ω0ω2ω4⋯ω2(L−1)⋮⋮⋮⋱⋮ω0ωL−1ω2(L−1)⋯ω(L−1)2⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (9)

with and is a diagonal matrix

 UL=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ω0ωN1{\huge 0}ω2N1{\huge 0}⋱ω(L−1)N1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (10)

Proof. From Eq. (7), we can rewrite as

 ⎡⎢ ⎢ ⎢ ⎢ ⎢⎣dm(N1)dm(N1+1)⋮dm(N1+L−1)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦T⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣eiN12πpLei(1+N1)2πpL⋮ei(L−1+N1)2πpL⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

Note that , and , it follows that

 1L[eiN12πpL,ei(1+N1)2πpL,⋯,ei(L−1+N1)2πpL]T=Product of F−1L % and (p+1)-th column of U−1L.

Hence we immediately obtain Eq. (8) which completes the proof.

Following the definition of , for , we further set by matrix

 Hn=[bk(n+jL−L)]jk, (11)

which means the element of is . Suppose that is invertible for every and we denote the inverse matrix as

 H−1n=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣q11(n)q12(n)⋯q1M(n)q21(n)q22(n)⋯q2M(n)⋮⋮⋮qM1(n)qM2(n)⋯qMM(n)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦.

Next we use the elements of to construct the interpolation functions. Let

 rm(n)={qmk(n+L−kL),if n∈Ik, k=1,2,⋯,M,0if n∉IN, (12)

and define

 ym(t)=∑n∈INrm(n)eint. (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

Let and be defined by (13). The shifted functions of can be expressed by

 [ym(t−t0),ym(t−t1),⋯,ym(t−tL−1)]=˜Vm(t)FLUL (14)

where are given respectively by (9), (10) and is defined as

 ˜Vm(t)=[vm,N1(t),vm,N1+1(t),⋯,vm,L+N1−1(t)]

with , for .

Proof. By the definition of and rearranging the summation terms of Eq. (13), we get

 ym(t−tp)=∑n∈INrm(n)ein(t−2πpL)=M∑k=1∑n∈I1rm(n+kL−L)ei(n+kL−L)(t−2πpL)=∑n∈I1(M∑k=1rm(n+kL−L)ei(n+kL−L)t)e−in2πpL.

From the piecewise defined in Eq. (12), it is easy to verify that for all and . It follows that can be rewritten as

 ym(t−tp)=∑n∈I1vm,n(t)e−in2πpL. (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

 Dn =[d1(n),d2(n),d3(n),⋯,dM(n)], Vn(t) =[v1,n(t),v2,n(t),v3,n(t),⋯,vM,n(t)].

By straightforward computations, we get

 AnHn =Dn, (16) Hn[Vn(t)]T =En(t). (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

Let and , be given above. If is invertible for all , then

 f(t)=1LM∑m=1L−1∑p=0gm(tp)ym(t−tp) (18)

where and is given by Eq. (13).

Proof. Multiplying Eq. (17) by and using Eq. (16), the expression (4) can be rewritten as

 f(t) =∑n∈I1Dn[Vn(t)]T=∑n∈I1M∑m=1dm(n)vm,n(t) =M∑m=1∑n∈I1dm(n)vm,n(t)=M∑m=1˜Dm[˜Vm(t)]T =M∑m=1˜DmF−1LU−1LULFL[˜Vm(t)]T (19) =1LM∑m=1L−1∑p=0gm(tp)ym(t−tp)

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,

 f(t)=12N+12N∑p=0f(tp)DN(t−tp) (20)

where

 DN(t)=⎧⎪⎨⎪⎩sin(12+N)tsin12t,t≠2kπ,2N+1t=2kπ;

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

 f′(t)=∑n∈INina(n)eint,  f′′(t)=∑n∈IN−n2a(n)eint.

Let , then should be defined by

 Hn=⎡⎢⎣1in−n21i(n+L)−(n+L)21i(n+2L)−(n+2L)2⎤⎥⎦.

By straightforward computations,

 H−1n=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣2L2+3Ln+n22L2−n(2L+n)L2n(L+n)2L2i(3L+2n)2L2−2i(L+n)L2i(L+2n)2L2−12L21L2−12L2⎤⎥ ⎥ ⎥ ⎥ ⎥⎦.

In order to satisfy the assumption that is divisible by , we have to assume that . Suppose that , then . Consequently, we obtain as follows:

 v1,n(t)= 2L2+3Ln+n22L2eint−n(2L+n)L2ei(n+L)t +n(L+n)2L2ei(n+2L)t, v2,n(t)= i(3L+2n)2L2eint−2i(L+n)L2ei(n+L)t +i(L+2n)2L2ei(n+2L)t, v3,n(t)= −12L2eint+1L2ei(n+L)t−12L2ei(n+2L)t.

By using and letting in Eq. (15), we have that equals

 sin3((N0+12)t)(N02+N0+1−(N0+1)N0cost)(2N0+1)2sin3(t2)

and

 y2(t) =∑n∈I1v2,n(t)=sin(t)sin3((N0+12)t)(2N0+1)2sin3(t2), y3(t) =∑n∈I1v3,n(t)=2sin3((N0+12)t)(2N0+1)2sin(t2).

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 anti-involution property (3), MIP can be applied to compute Hilbert transform as well.

###### Example 3

Let and . Note that does not satisfy the conditions of Theorem 1 because it is not invertible at . Nevertheless, By Remark 1, this inconsistence can be avoided if . Since for , then equals

 ∑0<|n|≤Nisgn(n)eint=−2csc(t2)sin(Nt2)sin(12(1+N)t)

and the sampling formula

 f(t)=−22N+12N∑p=0Hf(tp)csc(t−tp2)            ×sin(N(t−tp)2)sin(1+N2(t−tp)) (21)

holds if and .

Note that satisfies the preconditions for establishment of Eq. (21) inherently. By substituting with in (21) and using (3), we conclude that the Hilbert transform of can be computed by

 Hf(t)=22N+12N∑p=0(f(tp)−a(0))csc(t−tp2)            ×sin(N(t−tp)2)sin(1+N2(t−tp)). (22)

Unlike (21), the formula (22) is valid for the case of . In fact, if , we have

 a(0)=12π∫Tf(t)dt=12N+12N∑p=0f(tp) (23)

by direct computation.

###### Remark 2

Given a real-valued -point discrete signal , . By replacing with and letting in the right hand side of (22), we get a discrete signal . The DFT for is

 Z(k)=⎧⎨⎩X(0),k=0,2X(k),1≤k≤N,0,N+1≤k≤2N,

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

 fA(t)=N∑n=0˜aneint

where for and for . Applying Theorem 1 to , we have the following. If , then

 f(t)+iHf(t)=11+NN∑p=0(f(tp)+iHf(tp))y(t−tp),

where and with

 yr(t)=csc(t2)cos(Nt2)sin(12(1+N)t)

and

 yi(t)=csc(t2)sin(Nt2)sin(12(1+N)t).

Moreover, if is real-valued, then

 f(t) =11+NN∑p=0f(tp)yr(t−tp)−Hf(tp)yi(t−tp), (24) Hf(t) =11+NN∑p=0f(tp)yi(t−tp)+Hf(tp)yr(t−tp). (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 real-valued 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

 Hn=[1−isgn(n)1−isgn(n+L)].

It is clear that

 H−1n=[1212−i2i2]  for  −N≤n≤−1,  H−10=[10