# A fast inverse polynomial reconstruction method based on conformal Fourier transformation.

1. INTRODUCTIONTraditional representation of a smooth and periodic function in the Fourier basis has been widely applied in the spectral method [1], the pseudospectral method [2,3], and other computational physics [4] and signal processing research domains [5]. But Gibbs phenomenon occurs when the Fourier representation method is applied to a non-periodic or discontinuous function [6]. Neither accuracy nor efficiency can be achieved anymore because the finite Fourier representation cannot converge at the discontinuities or the boundaries of the function. Consequently, Gibbs phenomenon greatly hinders Fourier-representation based applications. A class of re-expansion methods have been developed to eliminate the Gibbs phenomenon in reconstructing discontinuous functions, with the locations of the discontinuities known a priori. They re-express the function in a non-periodic basis set, such as the Gegenbauer polynomial basis set [6, 7], the Pade Fourier rational basis function [8] and the Freund polynomial basis set [9], to obtain the exponential convergence over the entire interval of the function, including the discontinuities or boundaries. However, the performance of these methods depends on the specific basis set and imposes several stringent requirements on the basis parameters [10, 11].

The Inverse Polynomial Reconstruction Method (IPRM) presents a basis-independent reconstruction frame to overcome the Gibbs phenomenon [12,15,17]. The IPRM reconstructs a function as a finite sum of polynomials, based on the equation between the Fourier series of the reconstructed function and that of the original function. The IPRM provides more flexibility and faster convergence than other re-expansion methods, but Jung and Shizgal found that the increment in the polynomial degree can dramatically deteriorate the conditioning of the transformation matrix and the reconstruction may fail to converge [13, 14].

The Generalized IPRM method (G-IPRM) has been proposed in [16] to avoid the ill-conditioning of the transformation matrix. The condition number of the transformation matrix in the G-IPRM is kept close to one by increasing the number of entries in the transformation matrix, specifically by increasing the number of Fourier series to be the square of the degree of polynomials. Thus, the well-conditioned transformation matrix makes the G-IPRM more robust, and significantly generalizes the application domains of the traditional IPRM, but with increased memory cost and computation complexity. The time complexity and the memory cost of the G-IPRM are both O([M.sup.3.sub.G]), where [M.sub.G] is the degree of polynomials, if the iterative Least Square QR (LSQR) method [18] is used to invert for the reconstruction coefficients. Such heavy computation complexity and memory cost are undesirable in many applications, e.g. the reconstruction of large-scale induced electric or current density in discontinuous media in computational electromagnetics [19-24] or in radar signal research domain [25, 26]. Therefore, it is necessary to reduce the computation time and the memory complexity in the G-IPRM.

In this paper, a fast IPRM method is presented to reduce the computation complexity and the memory cost of the G-IPRM. The framework of the fast IPRM is modified by reconstructing the function in the discretized elements, so that the Conformal Fourier Transform (CFT) [28] and the Chirp Z-Transform (CZT) [29, 30] algorithms can be combined with the LSQR to reduce the computation complexity. The memory cost of the fast IPRM is reduced due to the transformation matrix being discretized, while the robustness of the reconstruction is preserved through configuring the relationships among the number of the discretized elements L, the degree of polynomials M and the number of Fourier series N. The computation complexity and memory cost of the fast IPRM are O(MNlog2L) and O(MN), respectively. Numerical results demonstrate the robustness and the efficiency of the fast IPRM.

This paper is organized as follows: The fast IPRM method is presented in Section 2. Numerical results, including reconstruction of a meromorphic function and the induced current density in discontinuous material distributions, are shown in Section 3. The conclusions are drawn in Section 4.

2. THE FAST IPRM METHOD.

Given the locations of discontinuities, the essence of the IPRM is solving an inverse problem from the already-known finite Fourier series of a discontinuous function and that of the polynomial basis functions. The objective of the proposed fast IPRM is to solve this inverse problem more efficiently through a discretization technique and two fast algorithms. There are three major steps in the fast IPRM: (1) modify the IPRM reconstruction framework by reconstructing the function in the discretized elements; (2) accelerate the reconstruction by combining two fast algorithms-CFT and CZT with the classic LSQR in the modified IPRM framework; (3) configure the discretization and reconstruction parameters to guarantee the robustness and efficiency enhancement of the fast IPRM.

In this section, first the fast IPRM will be presented for reconstructing a single-region discontinuous function;, and then it will be extended to a multiple-region discontinuous function.

2.1. The Fast IPRM Method for a Single-region Discontinuous Function

As a discontinuous function can be divided into multiple regions of piecewise continuous functions, we first consider the fast IPRM method for only one region of such a function.

2.1.1. Modification of the IPRM Framework by Discretization

Assume that f(x) is continuous in the region [OMEGA] = [a, b], b > a and zero outside. F = [[F.sub.-N/2],...,[F.sub.N/2-1]] is the vector of the first N Fourier coefficients of f(x) and the superscript T denotes matrix transpose.

The region [OMEGA] is equally segmented into L uniformly-spaced elements [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where [[OMEGA].sub.l] = [[a.sub.l],[b.sub.l]], and [[OMEGA].sub.i] [intersection] [[OMEGA].sub.j] = 0 if i [not equal to] j.

The boundaries of the lth element are

[a.sub.l] = a + (l - 1) q, [b.sub.l] = a + lq (1)

where q is the length of each element, q = [b.sub.l] - [a.sub.l] = b-a / L.

The function f(x) is locally reconstructed in these uniformly-discretized elements

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

where [[??].sub.l](x) is the locally-reconstructed function in the 1th element [[OMEGA].sub.l]. [??](x) is the reconstructed function in the region [OMEGA], and [direct sum] denotes that [??](x) is a collection of [[??].sub.i](x). [[xi].sup.l](x) = 2x-([b.sub.l]+[a.sub.l]) / q is the linear mapping from [[OMEGA].sub.l] = [[a.sub.l], [b.sub.l]] to [-1,1]. M is the degree of the polynomial, [C.sub.m](-) is the mth order polynomial, and [g.sup.l] = [[g.sup.l.sub.0] [g.sup.l.sub.1]...[g.sup.l.sub.M-1]] is the vector of coefficients for reconstructing [f.sub.l](x).

The reconstruction coefficients can be obtained by equating the Fourier series of the reconstructed function [??](x) and that of the original function f(x)

W x g = F (3)

where g is the vector of reconstruction coefficients constituted with that of the local reconstructions g = [[[g.sup.1] [g.sup.2]...[g.sup.L]].sup.T], and W is the transformation matrix

W = [[w.sup.1]...[w.sup.l]...[w.sup.L]] (4)

where [w.sup.l] is the sub-matrix of Fourier series of polynomials for the local reconstruction [??](x).

The entry of the sub-matrix [w.sup.l] at the n + N/2 + 1th row and the m + 1th column is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and [f.sub.0] = 1/(b - a) is the frequency resolution. n = -N/2,..., N/2-1, m = 0,1,..., M - 1, l = 1,..., L, and i is the imaginary unit.

The modified IPRM framework for reconstructing single-region discontinuous function is achieved by Equations (1)-(5). The reconstruction framework of the G-IPRM is a special case with L = 1. From (4) and (5), it can be seen that, due to the uniform discretization and the local reconstruction, the transformation matrix W is uniformly discretized into L sub-matrices, which are differentiated only by the exponential scaling term independent of the polynomial order m. Given the Fourier series F, the aim of the fast IPRM is to efficiently solve the reconstruction coefficients g in (3) by taking advantage of the uniformly-discretized transformation matrix W in the modified IPRM framework.

2.1.2. Acceleration of the Reconstruction by the CFT and the CZT

When the iterative LSQR is utilized to solve the Equation (3), the computation complexity of the IPRM is determined by the number of the multiplication operations to calculate the following two matrix-vector products

y = W x u (6)

z = [W.sup.[dagger]] x v (7)

where the superscript [dagger] denotes the complex conjugate transpose, u and v are already known vectors, u = [[u.sup.1]...[u.sup.l]...[[u.sup.L]].sup.T], [u.sup.l] = [[u.sup.l.sub.0] [u.sup.l.sub.1]...[u.sup.l.sub.M-1]], v = [[v.sub.-N/2]...[V.sub.N/2-1]]; y and z are respectively N x 1 and LM x 1 dimension unknown vectors.

The computation complexities of (6) and (7) by direct multiplication are both O (LMN), which is very expensive. Fortunately, in the modified framework, taking advantage of the uniformly discretized transformation matrix W, the computation complexity can be reduced by employing fast algorithms to evaluate products in (6) and (7).

Substituting (5) into (6), the (n + N/2 + 1)th entry of the vector y is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

where n = -N/2,..., N/2 - 1, [[gamma].sup.l](x) is the (M - 1) order polynomial function in the lth element given by [[gamma].sup.l](x) = [M-1.summation over (m=0)] [u.sup.l.sub.m] x [C.sub.m] ([[xi].sup.l](x)), x [member of] [[a.sub.l], [b.sub.l]].

The product vector y in (8) is the summation of Fourier series of the L piecewise continuous functions [[gamma].sup.l](x), [[gamma].sup.2](x),..., [[gamma].sup.L](x). These functions are in the same polynomial degree M and have the same interval space q, thus they can be considered as the uniformly-spaced elements of a discontinuous polynomial function [gamma](x), [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in [a, b] and zero outside.

Instead of calculating the Fourier series of the L piecewise continuous functions [[gamma].sup.l](x), [[gamma].sup.2](x),..., [[gamma].sup.L](x) independently, the product vector y can be efficiently and accurately achieved by calculating the Fourier series of [gamma](x) using the CFT algorithm [28], with the computation complexity of O (MN [log.sub.2] 2L). The CFT algorithm, in which uniform discretization and high-order polynomial interpolation techniques are utilized, is a fast and accurate method to compute the Fourier series of a discontinuous function. Since the discretization is conformal to the discontinuities, the CFT algorithm can accelerate the evaluation of y. Details of the CFT algorithm can be found in [28].

The other product vector to be calculated z in (7) is composed of L sub-vectors z = [[[z.sup.1]...[z.sup.l]...[z.sup.L]].sup.T] and [z.sup.l] = [[z.sup.l.sub.0] [z.sup.l.sub.1]...[z.sup.l.sub.M-1]], l = 1, 2,...,L. Substituting (5) into (7), the (m + 1)th entry of [z.sup.l] is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

where [d.sub.m](n) = [v.sub.n] x [([w.sup.1.sub.n,m]).sup.*] and the superscript x denotes the complex conjugate.

Due to the transformation matrix W being uniformly discretized, for each constant integer m, [z.sup.1.sub.m], [z.sup.2.sub.m],..., [z.sup.L.sub.m] are the first L discrete Fourier series of the N-point discrete function [d.sub.m](x), with the frequency resolution [f.sub.0]q = 1/L and L [much less than] N. Therefore, the product vector z can be efficiently evaluated by performing the Chirp Z-Transform (CZT) [29, 30] on the discrete function dm(n), for each m, m = 0,1,..., M - 1, hence the computation complexity to evaluate z is O (MN [log.sub.2] 2L).

By combining with the CZT and the CFT algorithms, the computation time in one iteration of the LSQR for evaluating the coefficient vector g is reduced from O (LMN) to O (MN [log.sub.2] 2L). Furthermore, it is not necessary to store the entire transformation matrix W, and only the sub-matrix [w.sup.1] and several vectors are required, leading to a memory requirement of O(MN) in the fast IPRM algorithm.

2.1.3. Discretization and Reconstruction Parameters

To guarantee the robustness as well as the efficiency enhancement of the fast IPRM, the relationship among the discretization and reconstruction parameters L, M and N are discussed.

The fast IPRM can possess the same robustness as the G-IPRM by keeping the condition number of the transformation matrix close to one. From (5), it can be seen that, the different sub-matrices of the transformation matrix W are independent of each other. Therefore, the condition number of transformation matrix W is equivalent to that of the matrix A

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. When [C.sub.m](x) is the normalized Legendre polynomial basis set, the entry in matrix A is [H.sub.m]([pi]n/L) = [(-i).sup.m] [(m + 1/2).sup.1/2] [j.sub.m]([pi]n/L), where [j.sub.m](x) is the mth order spherical Bessel function of the first kind [27]. According to Theorem 4.1 and Theorem 4.2 in [16], the condition number of the matrix A is decided by N/L. If N/L is the square of the degree M, i.e., if

N = L[M.sup.2] (11)

then the condition number of the transformation matrix W is kept close to one. This not only preserves the robustness of the fast IPRM as that of the G-IPRM, but also ensures that the numbers of iterations of the two methods to evaluate g are close to each other once the termination criterion of the LSQR [epsilon] is determined [16]. Thus, the reduction of the computation complexity in the fast IPRM is decided by that in one iteration of the LSQR.

To guarantee the efficiency improvement of the fast IPRM over the G-IPRM, the relationship among the reconstruction parameters of the fast IPRM L, M and that of the G-IPRM [M.sub.G], is now considered. As known from the convergence theory for polynomial expansion [1] and Theorem 5.1 in [16], the IPRM converges faster in a smaller interval. Because of the discretization and local reconstruction of the fast IPRM, given a certain accuracy requirement, M < [M.sub.G]. Generally, a larger element number L accompanies with a lower polynomial degree M. For a given reconstruction accuracy requirement, L and M can be selected to satisfy

[L.sup.2/3] M < [M.sub.G] (12)

From (12), the conditions L[M.sup.3] < [L.sup.2][M.sup.3] < [M.sup.3.sub.G] and L[M.sup.3] [log.sub.2] 2L < [L.sup.2][M.sup.3] < [M.sup.3.sub.G] are satisfied. Therefore, with the same robustness and the same reconstruction accuracy requirement, the fast IPRM, whose computation complexity and memory cost are O (L[M.sup.3] [log.sub.2] L) and O (L[M.sup.3]), respectively, can achieve much higher computation and memory cost efficiency than the G-IPRM, whose one-iteration computation complexity and memory cost are both O [M.sup.3.sub.G].

2.2. The Fast IPRM method for a Multiple-region Discontinuous Function

In this subsection, we will extend the fast IPRM to reconstruct a multiple-region discontinuous function from its Fourier coefficients.

Assume that f(x) is a multiple-region discontinuous function in the interval [OMEGA] = [[d.sub.0],[d.sub.s]], with S > 1 and zero outside. The discontinuities of f(x) locate at [[d.sub.0], [d.sub.1],..., [d.sub.s]], and [d.sub.0] < [d.sub.1] < ... < [d.sub.s]. The first N Fourier series of f(x) is F = [[F.sub.-N/2]...[F.sub.N/2-1]][.sup.T].

To obtain the modified IPRM framework for reconstructing f(x), firstly the interval of the function [OMEGA] is divided into S regions within which the function is continuous. Then by the same procedures as described in Subsection 2.1, each region is segmented into [L.sub.1], [L.sub.2],..., [L.sub.S] uniformly elements, where the function is locally reconstructed as the polynomials with degree of [M.sub.1], [M.sub.2],..., [M.sub.S], respectively.

The transformation matrix W, with the size of N x [S.summation over (s=1)] [L.sub.s][M.sub.s], is composed of [S.summation over (s=1)] [L.sub.s] sub-matrices of Fourier series of polynomial basis for local reconstruction in each element

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

where [w.sup.l.sub.s] is the submatrix of Fourier series, with the size of N x [M.sub.s], for local reconstruction of the lth element in the sth region, s = 1,..., S, l = 1,..., L. The entry of the sub-matrix wis at the n + N/2 + 1th row and the m + 1th column is

[w.sup.l.sub.s,n,m] = exp [-i2[pi]n (l - 1)[f.sub.0][q.sub.s]] x [w.sup.1.sub.s,n,m] (14)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. [q.sub.s] = ([d.sub.s] - [d.sub.s-1])/[L.sub.s] is element length in the sth region, and [f.sub.0] = 1/([d.sub.S] - [d.sub.0]) is the frequency resolution.

The reconstruction of the multiple-region discontinuous function is also an inverse problem to evaluate the coefficients g based on Equation (3). As known from (13) and (14), the transformation matrix W is uniformly discretized in each region. Therefore, similar to the single-region case, the CFT and the CZT can be combined with the LSQR to accelerate the evaluation of the reconstruction coefficients, and the computation complexity is O (N [S.summation over (s=1)] [M.sub.s] [log.sub.2] 2[L.sub.s]). Meanwhile, the memory cost is O (N [S.summation over (s=1)] [M.sub.s]) since only sub-matrices [w.sup.1.sub.1], [w.sup.1.sub.2], ..., [w.sup.1.sub.S] and several vectors are required to be stored.

For the reconstruction of the multiple-region discontinuous function, the robustness and the efficiency of the fast IPRM are further guaranteed by selecting the reconstruction parameters to meet the following conditions

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

where [M.sub.G,s] is the polynomial degree for reconstructing the sth region by the G-IPRM, with the same reconstruction accuracy as the fast IPRM, and [N.sub.G] is the number of Fourier series in the G-IPRM, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

3. NUMERICAL RESULTS

In this section, three examples are simulated to demonstrate the performance of the fast IPRM. Example in Subsection 3.1 takes the same meromorphic function as in [16] to compare the fast IPRM and the G-IPRM proposed in [16]. Subsection 3.2 compares the reconstruction results of a single-region induced current density by the fast IPRM and the G-IPRM, respectively. Subsection 3.3 shows the reconstruction results by the fast IPRM and the G-IPRM for an induced current density whose support domain includes multiple piecewise continuous regions.

To compare the efficiency of the fast IPRM and the G-IPRM, in these three examples, the memory costs of the two methods are estimated by the number of the entries in the transformation matrix that need to be stored, and their computation complexities are measured by the number of the complex multiplication operations to obtain the reconstruction coefficients g.

3.1. Reconstruction of a Single-region Meromorphic Function

The analytic expression of the function is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This function, as shown in Fig. 1(a), is a single-region discontinuous function with discontinuities located at x = [+ or -]1.

In the fast IPRM, for different degree of polynomials M, to meet condition (12), the number of the discretized elements is configured to be L = 2. For the fixed L, according to (11), the number of the Fourier series is N = 2[M.sup.2].

As shown in Fig. 1(b), similar to the G-IPRM, the condition number of the transformation matrix W in the fast IPRM is kept close to one with the increasing of polynomial degree M. This not only ensures the robustness of the fast IPRM, but also makes the number of iterations in the fast IPRM almost the same as that in the G-IPRM when the termination criterion of the LSQR algorithm [epsilon] is decided, as shown in Fig. 1(c). Hence the difference of computation complexity between these two methods is determined by the computation complexity in one iteration.

The accuracy of the reconstruction is measured by the relative maximum error (RME) [[parallel] f (x) - [??] (x) [parallel].sub.infinity]/[[parallel]f (x)[parallel].sub.infinity]. From Fig. 1(d), it can be seen that, the fast IPRM achieve the same reconstruction accuracy with a lower polynomial degree than the G-IPRM.

The memory cost and the computation complexity of the two methods are shown in Fig. 1(e) and Fig. 1(f), respectively. It can be seen that, although the computation complexity of the fast IPRM is not obviously improved, the memory cost of the fast IPRM is greatly reduced, e.g., for a relative maximum error as [10.sup.-5], the memory cost in the fast IPRM is about 1/12 of that in the G-IPRM.

[FIGURE 1 OMITTED]

Therefore, for this single-region meromorphic function, the fast IPRM not only inherits the robustness of the G-IPRM but also possesses higher memory efficiency than the G-IPRM.

3.2. Reconstruction of a Single-region Induced Current Density

In this subsection, an induced electric current density due to a plane wave at normal incidence to a three-layer medium with frequency [f.sub.c] = 5 GHz in a vacuum background is reconstructed. The layer interfaces are located at 2.2 m and 4.4 m, respectively, and the relative magnetic permeability for all three layers is 1. The relative permittivity is [[epsilon].sub.r] = 1, 3 and 1 respectively. The conductivity for the center layer is [sigma] = 0.003 S/m. The wave number in the central layer is 63.5085. As shown in Fig. 2(a), this wave distribution is a single-region discontinuous function with discontinuities located at x = 2.2 m and x = 4.4 m, respectively.

To guarantee the efficiency improvement of the fast IPRM, L and M are selected to satisfy (12). As shown in Fig. 2(h) multiple values of L can meet formula (12) and reduce the computation complexity compared to the G-IPRM. In the simulation, the value of L is taken as the half of the wave numbers in the center layer, i.e., L = 32, and for the fixed L, the number of the Fourier series varies with the degree of polynomial, N = 32[M.sup.2].

[FIGURE 2 OMITTED]

Figure 2(b) shows that, just as in the G-IPRM, the condition number of the transformation matrix W in the fast IPRM is kept close to one with the increasing of polynomial degree. Fig. 2(c) demonstrates that the numbers of iterations in the fast IPRM and the G-IPRM are close to each other, with a fixed termination criterion [epsilon] in the LSQR. Thus the difference of computation complexity between these two methods is determined by that in one iteration of the LSQR.

As shown in Fig. 2(d), the G-IPRM method converges much slower and requires higher degree polynomials, whereas the fast IPRM with elements number L = 32, needs a much lower polynomial degree to obtain the same reconstruction accuracy. Fig. 2(g) shows that the RME of the fast IPRM is reduced with the increase of L for a fixed value of M. Meanwhile, a larger element number L accompanies a lower polynomial degree M for a given accuracy requirement.

The memory cost and the computation complexity are significantly reduced in the fast IPRM compared with the G-IPRM, as shown in Figs. 2(e) and (f). With the relative maximum error [10.sup.-5], the computation complexity and the memory cost in fast IPRM are only about 1 /4 and 1 /90 of that in the G-IPRM, respectively.

3.3. Reconstruction of a Multiple-region Induced Current Density

A multiple-region induced electric current density due to a plane wave at normal incidence to a five-layer medium with frequency [f.sub.c] = 5 GHz in a vacuum background, is reconstructed in this subsection. The layer interfaces are located at 0.8 m, 1 m, 5.4 m and 5.8 m, respectively. The relative permittivity is [[epsilon].sub.r] = 1, 5, 3, 4 and 1, respectively. The conductivity for the center three layers is [sigma] = 0.005, 0.003 and 0.004 S/m, respectively, and the wave numbers in the three center layers are 7.4536, 127.02, 13.3333, respectively. As shown in Fig. 3(a), the wave distribution is a three-region discontinuous function with discontinuities located at 0.8 m, 1m, 5.4 m and 5.8 m, respectively.

In the fast IPRM, to meet the condition in (16), we configure the element number as the half of the wave numbers in the central three layers, [L.sub.1] = 4, [L.sub.2] = 64 and [L.sub.3] = 7. The degrees of polynomials in the three regions are the same [M.sub.1] = [M.sub.2] = [M.sub.3] = M, and according to (15), the number of the Fourier series is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

As in the previous two numerical examples, the robustness and the efficiency of the fast IPRM are illustrated in Fig. 3(b)-Fig. 3(f). Fig. 3(b) shows that, in the fast IPRM, the condition number of transformation matrix W is kept close to one with the increasing of polynomial degree, while that in the G-IPRM is slightly larger. Fig. 3(c) demonstrates that the numbers of iterations in the fast IPRM and the G-IPRM are close to each other. As shown in Fig. 3(d), the polynomial degree in the fast IPRM is much lower than that in the G-IPRM. From Fig. 3(e) and Fig. 3(f), it can be found that, the memory cost and the computation complexity are significantly reduced in the fast IPRM compared with the G-IPRM, e.g. for a maximum relative error as 10-5, the computation complexity and the memory cost in fast IPRM are about 1/5 and 1/94 of that in the G-IPRM, respectively.

[FIGURE 3 OMITTED]

3.4. Discussion

The above numerical examples demonstrate that the fast IPRM method works efficiently for eliminating the Gibbs phenomenon of both single-region and multiple-region discontinuous functions. It is also seen that, provided that the number of the discretized elements L satisfy (12) or (16), the efficiency improvement of the fast IPRM over the G-IPRM is more distinct with the increment of L. For the multiple-region case, the improvement of the efficiency in the fast IPRM over the G-IPRM is decided by the region, which has the largest number of the discretized elements.

4. CONCLUSION

In this paper, we propose a fast IPRM method to efficiently overcome Gibbs phenomenon in the reconstruction of discontinuous or non-periodic functions. The Conformal Fourier Transform (CFT) and the Chirp Z-Transform (CZT) algorithms are combined in the modified IPRM framework to decrease both computation time and memory complexity. For a discontinuous function with S piecewise-continuous regions, the memory cost and the computation complexity of the fast IPRM method are O (N [S.summation over (s=1)] [M.sub.s]) and O (N [S.summation over (s=1)] [M.sub.s] [log.sub.2] 2[L.sub.s]), respectively. Numerical results show that the fast IPRM method works well for both single-region functions and multiple-region functions. Therefore, the fast IPRM method can be useful in pseudospectral time-domain methods and in volume integral equation solvers for discontinuous material distributions or other areas where Gibbs phenomenon needs to be efficiently eliminated.

ACKNOWLEDGMENT

ZL, JYY are supported by the Natural Science Foundation of China (NSFC) under Grants 61102143. ZL is grateful to the support of China Scholar Council (CSC) for the chance of researching in Duke University.

Received 20 September 2011, Accepted 3 November 2011, Scheduled 16 November 2011

REFERENCES

[1.] Boyd, J. P., Chebyshev and Fourier Spectral Methods, Dover publications, 2000.

[2.] Fornberg, B., A Practical Guide to Pseudospectral Methods, Cambridge University Press, 1996.

[3.] Liu, Q. H., "The PSTD algorithm: A time-domain method requiring only two cells per wavelength," Microwave and Optical Technology Letters, Vol. 15, 158-165, 1997.

[4.] Giordano, N. J. and H. Nakanishi, Computational Physics, Prentice Hall Publishing, 1997.

[5.] Gonzalez, R. C. and R. E. Woods, Digital Image Processing, Addison-Wesley Publishing, 1992.

[6.] Gottlieb, D. and C. W. Shu, "On the Gibbs phenomenon and its resolution," SIAM Review, Vol. 39, 644-668, 1997.

[7.] Gottlieb, D., C. W. Shu, S. Alex, and H. Vandeven, "On the Gibbs phenomenon I: Recovering exponential accuracy from the fourier partial sum of a nonperiodic analytic function," Journal of Computational and Applied Mathematics, Vol. 43, 81-98, 1992.

[8.] Driscoll, T. A. and B. Fornberg, "A pade-based algorithm for overcoming Gibbs phenomenon," Numerical Algorithms, Vol. 26, 77-92, 2001.

[9.] Gelb, A. and J. Tanner, "Robust reprojection methods for the resolution of the Gibbs phenomenon," Applied and Computational Harmonic Analysis, Vol. 20, 3-25, 2006.

[10.] Boyd, J. P., "Trouble with Gegenbauer reconstruction for defeating Gibbs' phenomenon: Runge phenomenon in the diagonal limit of Gegenbauer polynomial approximations," Journal of Computational Physics, Vol. 204, 253-264, 2005.

[11.] Min, M., T. Lee, P. F. Fischer, and S. K. Gray, "Fourier spectral simulations and Gegenbauer reconstructions for electromagnetic waves in the presence of a metal nanoparticle," Journal of Computational Physics, Vol. 213, 730-747, 2006.

[12.] Shizgal, B. D. and J. Jung, "Towards the resolution of the Gibbs phenomena," Journal ofComputational and Applied Mathematics, Vol. 161, 41-65, 2003.

[13.] Jung, J. and B. D. Shizgal, "Generalization of the inverse polynomial reconstruction method in the resolution of the Gibbs phenomenon," Journal of Computational and Applied Mathematics, Vol. 172, 131-151, 2004.

[14.] Jung, J. and B. D. Shizgal, "On the numerical convergence with the inverse polynomial reconstruction method for the resolution of the Gibbs phenomenon," Journal of Computational Physics, Vol. 224, 477-488, 2007.

[15.] Jung, J. and B. D. Shizgal, "Inverse polynomial reconstruction of two-dimensional fourier images," Journal of Science Computing, Vol. 25, 367-399, 2005.

[16.] Hrycak, T. and K. Grochenig, "Pseudospectral fourier reconstruction with the modified inverse polynomial reconstruction method," Journal of Computational Physics, Vol. 229, 933-946, 2010.

[17.] Jung, J., "A hybrid method for the resolution of the Gibbs phenomenon," Lecture Notes in Computational Science and Engineering, Vol. 76, 219-227, 2011.

[18.] Paige, C. C. and M. A. Saunders, "LSQR: An algorithm for sparse linear equations and sparse least squares," ACM Transactions on Mathematical Software, Vol. 8, 43-71, 1982.

[19.] Carvalho, S. and L. S. Mendes, "Scattering of EM waves by inhomogeneous dielectrics with the use of the method of moments and 3-D solenoidal basis functions," Microwave and Optical Technology Letters, Vol. 23, 42-46, 1999.

[20.] Li, M.-K. and W. C. Chew, "Applying divergence-free condition in solving the volume integral equation," Progress In Electromagnetics Research, Vol. 57, 311-333, 2006.

[21.] Fan, Z., R.-S. Chen, H. Chen, and D.-Z. Ding, "Weak form nonuniform fast Fourier transform method for solving volume integral equations," Progress In Electromagnetics Research, Vol. 89, 275-289, 2009.

[22.] Hu, L., L.-W. Li, and T. S. Yeo, "Analysis of scattering by large inhomogeneous bi-anisotropic objects using AIM," Progress In Electromagnetics Research, Vol. 99, 21-36, 2009.

[23.] Taboada, J. M., M. G. Araujo, J. M. Bertolo, L. Landesa, F. Obelleiro, and J. L. Rodriguez, "MLFMA-FFT parallel algorithm for the solution of large-scale problems in electromagnetics," Progress In Electromagnetics Research, Vol. 105, 15-30, 2010.

[24.] Ergul, O., T. Malas, and L. Gurel, "Solutions of large-scale electromagnetics problems using an iterative inner-outer scheme with ordinary and approximate multilevel fast multipole algorithms," Progress In Electromagnetics Research, Vol. 106, 203-223, 2010.

[25.] Huang, Y., Y. Liu, Q. H. Liu, and J. Zhang, "Improved 3-D GPR detection by NUFFT combined with MPD method," Progress In Electromagnetics Research, Vol. 103, 185-199, 2010.

[26.] Zhu, X., Z. Zhao, W. Yang, Y. Zhang, Z.-P. Nie, and Q. H. Liu, "Iterative time-reversal mirror method for imaging the buried object beneath rough ground surface," Progress In Electromagnetics Research, Vol. 117, 19-33, 2011.

[27.] Abramowitz, M. and I. A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Government Printing Office Publishing, 1972.

[28.] Zhu, C. H., Q. H. Liu, Y. Shen, and L. J. Liu, "A high accuracy conformal method for evaluating the discontinuous fourier transform," Progress In Electromagnetics Research, Vol. 109, 425-440, 2010.

[29.] Rabiner, L. R., R. W. Schafer, and C. M. Rader, "The chirp Z-transform algorithm and its application," IEEE Transaction on Audio Electroacoust, Vol. 17, 86-92, 1969.

[30.] Franceschetti, G., R. Lanari, and E. S. Marzouk, "A new two-dimensional squint mode SAR processor," IEEE Transactions on Aerospace and Electronic Systems, Vol. 32, 854-863, 1996.

Z. Liu (1), *, Q. H. Liu (2), C. H. Zhu (3), and J. Y. Yang (1)

(1) School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China

(2) Department of Electrical and Computer Engineering, Duke University, Durham, NC 27708, USA

(3) Department of Control Science and Engineering, Harbin Institute of Technology, 92 West Street, Harbin 150001, China

* Corresponding author: Zhe Liu (zl52@duke.edu).

Printer friendly Cite/link Email Feedback | |

Author: | Liu, Z.; Liu, Q.H.; Zhu, C.H.; Yang, J.Y. |
---|---|

Publication: | Progress In Electromagnetics Research |

Article Type: | Report |

Geographic Code: | 9CHIN |

Date: | Dec 1, 2011 |

Words: | 5627 |

Previous Article: | A hybrid method based on differential evolution and continuous ant colony optimization and its application on wideband antenna design. |

Next Article: | Mapping the SBR and TW-ILDCs to heterogeneous CPU-GPU architecture for fast computation of electromagnetic scattering. |

Topics: |