NumSBT: A subroutine for calculating spherical Bessel transforms numerically

A previous subroutine, LSFBTR, for computing numerical spherical Bessel (Hankel) transforms is updated with several improvements and modifications. The procedure is applicable if the input radial function and the output transform are defined on logarithmic meshes and if the input function satisfies reasonable smoothness conditions. Important aspects of the procedure are that it is simply implemented with two successive applications of the fast Fourier transform, and it yields accurate results at very large values of the transform variable. Applications to the evaluation of overlap integrals and the Coulomb potential of multipolar charge distributions are described.

Program summary

Program title: NumSBTCatalogue identifier: AANZ_v2_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AANZ_v2_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 476No. of bytes in distributed program, including test data, etc.: 4451Distribution format: tar.gzProgramming language: Fortran 90Computer: GenericOperating system: LinuxClassification: 4.6Catalogue identifier of previous version: AANZ_v1_0Journal reference of previous version: Comput. Phys. Comm. 30 (1983) 93Does the new version supersede the previous version?: NoNature of problem: This program is a subroutine which, for a function defined numerically on a logarithmic mesh in the radial coordinate, generates the spherical Bessel, or Hankel, transform on a logarithmic mesh in the transform variable. Accurate results for large values of the transform variable are obtained, that would not be otherwise obtainable.Solution method: The program applies a procedure proposed by the author [1] that treats the problem as a convolution. The calculation then requires two applications of the fast Fourier transform method.Reasons for new version: The method of computing the transform at small values of the transform variable has been substantially changed and the whole procedure simplified. In addition, the possibility of computing the transform for a single transform variable value has been incorporated. The code has also been converted to Fortran 90 from Fortran 77.Restrictions: The procedure is most applicable to smooth functions defined on (0,∞) with a limited number of nodes.Running time: The example provided with the distribution takes a few seconds to execute.References:[1] J.D. Talman, J. Comp. Phys. 29 (1978) 35.     

numer, a code for Numerov integrations of Coulomb functions

The code numer is used for numerical integrations of Coulomb radial wave functions using the Numerov method. The input specifies a function and its derivative to start integrations, an integration range and an accuracy parameter ac such that the accumulated error is no larger than ac. Alternative input is initial function, integration step, and function after first step. For positive energies, options exist to use either the atomic-physics variables (?,r) or the nuclear-physics variables (η,ρ).     

Calculation of Fourier transform using symmetrical properties

A matrix method which computes discrete Fourier transforms using a digital computer program is presented in this paper. The proposed technique takes advantage of symmetry of the complex functions about the real and imaginary axes to reduce the number of calculations necessary in a given Fourier transform. Computationally the method described is not as efficient, especially for a large N, as the well-known Cooley-Tukey method. However, it differs from the Cooley-Tukey formulation in two notable respects: first, the present method is not as restrictive in the selection of values of N as the Cooley-Tukey method; and second, the calculations can proceed with the first value of the time function, thus eliminating the need for storing data before beginning with the transform calculations as is the case with the Cooley-Tukey method. In a number of applications, realizing these two conditions is more important than computational efficiency.The logic of a computer program which calculates the Fourier transform using the symmetry properties is described by a flow chart. This paper also includes numerical examples, using the computer program, of a Fourier transform from the time domain to the frequency domain. In addition, the program calculates the inverse Fourier transform, reconstructing the original time function from its frequency contents.     

An algorithm for the Fourier-Bessel transform

It is shown that the Fourier-Bessel (Hankel) transform may be calculated by means of a single one-dimensional Fourier transform followed by repeated summations of preselected Fourier components. For an N/2 Fourier/Bessel transformation of a complex function this algorithm involves of the order of N log₂N multiplications and N² additions and it requires about 4N storage allocation. Numerical test cases illustrate various features of the proposed method.     

Numerical computation of infinite Bessel transforms with high frequency

In this paper numerical evaluation of infinite Bessel transforms with high frequency is considered. We first derive an asymptotic formula by using the integration by parts. Next we use an interpolation formula to evaluate the infinite Bessel transforms by choosing the suitable basis and nodes. The corresponding error results are proved, and numerical examples are shown to illustrate the efficiency and accuracy of the presented formulae.     

A numerical transform with epsilon-quadrature

Early z-transform investigators originated tuned quadrature. Later, Smith applied tuned integration to dynamical system simulation. Herein, epsilon-quadrature extends trapezoidal quadrature to a numerical transform associated with tuned integration. Examples demonstrate instability for a digital Butterworth filter and the accuracy improvement scenario for the new numerical transform.     

The fast Fourier transform for general order

An elementary and transparent representation of the fast Fourier transform is given. Instead of using the usual and highly algebraic approach it is shown how a Fourier transform of the ordern=p·m can be reduced top Fourier transforms of orderm by performing essentiallym Fourier transforms of orderp on the data. The resulting process is discussed in more detail forn=3 ^q andn=5 ^q . The problem of retrieval of the wanted coefficients from the final data is solved by a simple argument. The generalization for an ordern equal to a product of powers of prime numbers is immediate.     

The circlet transform: A robust tool for detecting features with circular shapes

We present a novel method for detecting circles on digital images. This transform is called the circlet transform and can be seen as an extension of classical 1D wavelets to 2D; each basic element is a circle convolved by a 1D oscillating function. In comparison with other circle-detector methods, mainly the Hough transform, the circlet transform takes into account the finite frequency aspect of the data; a circular shape is not restricted to a circle but has a certain width. The transform operates directly on image gradient and does not need further binary segmentation. The implementation is efficient as it consists of a few fast Fourier transforms. The circlet transform is coupled with a soft-thresholding process and applied to a series of real images from different fields: ophthalmology, astronomy and oceanography. The results show the effectiveness of the method to deal with real images with blurry edges.     

The Partial Fast Fourier Transform

An efficient algorithm for computing the one-dimensional partial fast Fourier transform \(f_j=\sum _{k=0}^{c(j)}e^{2\pi ijk/N} F_k\) is presented. Naive computation of the partial fast Fourier transform requires \({\mathcal O}(N^2)\) arithmetic operations for input data of length N. Unlike the standard fast Fourier transform, the partial fast Fourier transform imposes on the frequency variable k a cutoff function c(j) that depends on the space variable j; this prevents one from directly applying standard FFT algorithms. It is shown that the space–frequency domain can be partitioned into rectangular and trapezoidal subdomains over which efficient algorithms can be developed. As in the previous work of Ying and Fomel (Multiscale Model Simul 8(1):110–124, 2009), the contribution from rectangular regions can be reduced to a series of fractional-phase Fourier transforms over squares, each of which can be reduced to a convolution. In this work, we demonstrate that the partial Fourier transform over trapezoidal domains can also be reduced to a convolution. Since the computational complexity of a dealiased convolution of N inputs is \({\mathcal O}(N\log N)\), a fast algorithm for the partial Fourier transform is achieved, with a lower overall coefficient than obtained by Ying and Fomel.     

Optimal Fourier filtering of a function that is strictly confined within a sphere

We present an alternative method to filter a distribution, that is strictly confined within a sphere of given radius rc, so that its Fourier transform is optimally confined within another sphere of radius kc. In electronic structure methods, it can be used to generate optimized pseudopotentials, pseudocore charge distributions, and pseudo atomic orbital basis sets.     

On the computation and the effectiveness of discrete sine transform

The effectiveness of the discrete sine transform (DST) in terms of residual correlation as developed by Hamidi and Pearl[1] is investigated. Based on this criterion, the DST is compared with other discrete transforms such as cosine (DCT)[2] and Fourier (DFT)[3]. The odd-even property of the DST is established. Based on this, the DST is computed through the Walsh-Hadamard transform (WHT). This development parallels that of Jones et al.[4] for other transforms such as slant (ST), DFT, DCT, etc.     

Numerical evaluation of the Hankel transform by using linear Legendre multi-wavelets

An efficient algorithm for evaluating the Hankel transform Fn(p) of order n of a function f(r) is given. As the continuous Legendre multi-wavelets forms an orthonormal basis for L²(R); we expand the part rf(r) of the integrand in its wavelet series reducing the Hankel transform integral as a series of Bessel functions multiplied by the wavelet coefficients of the input function. Numerical examples are given to illustrate the efficiency of the proposed method.     

The multiple-parameter fractional Fourier transform

The fractional Fourier transform （FRFT） has multiplicity, which is intrinsic in fractional operator. A new source for the multiplicity of the weight-type fractional Fourier transform （WFRFT） is proposed, which can generalize the weight coefficients of WFRFT to contain two vector parameters m,n ∈ Z^M . Therefore a generalized fractional Fourier transform can be defined, which is denoted by the multiple-parameter fractional Fourier transform （MPFRFT）. It enlarges the multiplicity of the FRFT, which not only includes the conventional FRFT and general multi-fractional Fourier transform as special cases, but also introduces new fractional Fourier transforms. It provides a unified framework for the FRFT, and the method is also available for fractionalizing other linear operators. In addition, numerical simulations of the MPFRFT on the Hermite-Gaussian and rectangular functions have been performed as a simple application of MPFRFT to signal processing.     

Error propagation in fourier transforms

The Fourier transform has long been of great use in simulating mathematical or physical phenomena, especially in signal theory. However the finite length representation of numbers introduces round-off errors in computing. Here, developing a new point of view on the topic, we give an evaluation of the total relative mean square error in the computation of direct and fast Fourier transforms using floating point artihmetic. Thus we show that in direct Fourier transforms the output noise-to-signal ratio is equivalent to N or N² according to whether the arithmetic is a rounding or a chopping one, whereas for fast Fourier transforms it is equivalent to log₂(N) or [log₂(N)]², with N being the number of points of the signal. Good agreement with numerical results is observed.     

Subspaces of FM~mlet transform

The subspaces of FMmlet transform are investigated. It is shown that some of the existing transforms like the Fourier transform, short-time Fourier transform, Gabor transform, wavelet transform, chirplet transform, the mean of signal, and the FM-1let transform, and the butterfly subspace are all special cases of FMmlet transform. Therefore the FMmlet transform is more flexible for delineating both the linear and nonlinear time-varying structures of a signal.     

Numerical inversion of z-transforms with application to polymerization kinetics

A numerical technique for determination of the polymer chain length distribution from mathematical models of polymerization systems is presented. The proposed method, which can be applied to the class of model equations whose z-transform exists, is based upon a z-transform inversion formula that has the form of an Inverse Discrete Fourier Transform. A computer algorithm that employs the FFT to efficiently obtain numerical results from the derived formulas is summarized and an error analysis is given. An example problem that describes the polycondensation polymerization kinetics of bisphenol-A to polycarbonate is introduced as a realistic application of the proposed method. It is shown that for this example, both the accuracy and computational efficiency of the proposed z-transform inversion technique are superior to another independent method for obtaining the polymer chain length distribution that is based upon direct numerical integration. A methodology for identification of polymerization kinetic rate parameters from experimental polymer molecular weight distributions is also suggested.     

Interactive noise-controlled boundary image matching using the time-series moving average transform

In this paper we propose a time-series matching-based approach that provides the interactive boundary image matching with noise control for a large-scale image database. To achieve the noise reduction effect in boundary image matching, we exploit the moving average transform of time-series matching. We are motivated by a simple intuition that the moving average transform might reduce the noise of boundary images as well as that of time-series data. To confirm this intuition, we first propose a new notion of k-order image matching, which applies the moving average transform to boundary image matching. A boundary image can be represented as a sequence in the time-series domain, and our k-order image matching identifies similar boundary images in this time-series domain by comparing the k-moving average transformed sequences. We then propose an index-based method that efficiently performs k-order image matching on a large image database, and formally prove its correctness. We also formally analyze the relationship of orders and their matching results and present an interactive approach of controlling the noise reduction effect. Experimental results show that our k-order image matching exploits the noise reduction effect well, and our index-based method outperforms the sequential scan by one or two orders of magnitude. These results indicate that our k-order image matching and its index-based solution provide a very practical way of realizing the noise control boundary image matching. To our best knowledge, the proposed interactive approach for large-scale image databases is the first attempt to solve the noise control problem in the time-series domain rather than the image domain by exploiting the efficient time-series matching techniques. Thus, our approach can be widely used in removing other types of distortions in image matching areas.     

Subspaces of FM<Superscript>m</Superscript>let transform

The subspaces of FM^mlet transform are investigated. It is shown that some of the existing transforms like the Fourier transform, short-time Fourier transform, Gabor transform, wavelet transform, chirplet transform, the mean of signal, and the FM⁻¹let transform, and the butterfly subspace are all special cases of FM^mlet transform. Therefore the FM^mlet transform is more flexible for delineating both the linear and nonlinear time-varying structures of a signal.