Fast algorithms for nonuniform Chirp-Fourier transform

1.   Introduction

The majority of digital signal processing theories currently are based on the uniform, ideal sampling signal model. With the development of science and technology, the range of processed signals has become very broad. The digital signals obtained in practical applications are generally nonuniform and nonideal signal sampling sequences. For example, in the process of radar signal acquisition, geographic data acquisition, and medical image signal acquisition, due to the presence of certain interference factors, it is not easy or impossible to obtain completely uniformly sampled signals. Meanwhile, in practical applications, due to the influence of signal receiving equipment and external interference factors, the sampling points of the obtained signal may not necessarily be the sampling points of the original signal. Therefore, nonuniform sampling has become one of the important hotspots in signal processing ^[1,2]. The Fourier transform (FT), which is very important in the Fourier analysis theory system, is the basis for the traditional signal processing. It has good performance in handling stationary signals. However, it has significant limitations to process nonstationary using FT. The Chirp signal is a common nonstationary signal. Therefore, numerous novel mathematical techniques, including the fractional Fourier transform (FrFT) ^[2,3,4], the Chirp-Fourier transform (CFT) ^[5,6,7,8], and the linear canonical transform (LCT) ^[9,10,11,12], have been proposed to process nonstationary signal. Among them, since the CFT has good mathematical qualities and fast algorithms, it is frequently used in nonstationary signal processing.

The application of the CFT relies on its discretization and fast algorithms. Discrete Chirp-Fourier transform (DCFT) was first proposed by Xia in 2000 ^[5]. However, this method requires that the number of sampling points is prime, and the parameters of frequency modulation must be integers or approximate integers, which has significant limitations in practical applications. Fan et al. derived modified DCFT (modified Chirp-Fourier transform, MDCFT) to improve the aforementioned limitations ^[6]. On the basis of these, Guo Yong et al. ^[7] gave a fast method based on CFT for estimating the parameters of Newton's ring. The segmented CFT fast algorithm was proposed by Zong Ying et al. in 2023 ^[8]. Simulation results demonstrate that the suggested algorithm can realize the fast computing of the local spectrum. Most of the existing CFT algorithms are based on the uniform sampling. For the nonuniform sampling, the above algorithms can't achieve good results. It is important to consider the nonuniform sampling associated with the CFT. Meanwhile there are fewer researches on fast algorithms for nonuniform CFT currently. In order to promote the CFT to be better applied in various fields, we will propose the nonuniform fast Chirp-Fourier transform (NUCFT) algorithm.

The rest of the paper is organized as follows: In Section 2, we give definitions of the CFT, Chirp-Fourier series, and associated lemmas. The main results are derived in Section 3. Finally, the conclusions are presented in Section 4.

2.   Preliminary

Definition 2.1. The CFT of a continuous signal $x(t)$ is defined as ^[13]

where $r = - \infty, \cdots, \infty.$

When $r = 0,$ the transform degenerates to the FT, i.e., ${X_r}\left(u \right) = \int_{ - \infty }^\infty {x\left(t \right){e^{ - j2\pi ut}}dt}.$

The inverse of the CFT for a continuous signal $x(t)$ is defined as

By discretizing the time and transform domains, different definitions of the discrete CFT can be obtained ^[7,14,15]. The most commonly used DCFT is the following.

Definition 2.2. The DCFT of a discrete sequence $x(n) = x(n\Delta t)$ is defined as ^[16]

where $m = 0, 1, \cdots, N - 1.$$\Delta t$ and $\Delta u$ are the time and frequency domain sampling intervals for time ${\text{t}}$ and transform domain $u$, respectively.

The inverse of the DCFT for a discrete sequence $x(n) = x(n\Delta t)$ is defined as

The definition of the Chirp-Fourier series for the signal $x(t)$ is obtained as follows.

Due to

The sequence $\left\{ { \cdots, {\varphi _{r, - 1}}(t), {\varphi _{r, 0}}(t), {\varphi _{r, 1}}(t), \cdots } \right\}$ is orthogonal in the interval $[-T/2, T/2]$. Since the CFS only employs finite-length functions and this basis function is an FM function, the CFS of a finite-length signal $x(t)$ is

where the coefficient of the CFS is ${{\text{c}}_k} = \frac{1}{{\sqrt T }}\int_{ - T/2}^{T/2} {x(t){e^{ - i(kt + r{t^2})}}} dt.$

The pertinent formulas that have been applied to derive the NUCFT are given as the following.

Lemma 2.1. For any real $c,$

Lemma 2.2. For any integer $k,$

Lemma 2.3. ^[17] For any real $b > 0$ and complex $z,$

Lemma 2.4. ^[17] For any real $b > 0$ and $a > 0,$

In the next section, the main results will be derived based on the above facts.

3.   Main results

3.1.   Exact statement of the problems

In the remainder of this paper, we will operate under the following assumptions.

(1) $\omega = \{ {\omega _0}, \cdots, {\omega _N}\}$ and $x = \{ {x_0}, \cdots, {x_{N - 1}}\}$ are finite sequences of real numbers.

(2) ${\omega _k} \in [-N/2, N/2-1]$, for $k = 0, \cdots, N - 1$.

(3) ${x_i} \in [- \pi, \pi]$, for $j = 0, \cdots, N - 1$.

(4) $\alpha = \{ {\alpha _0}, \cdots, {\alpha _{N - 1}}\},$$f = \{ {f_{ - N/2}}, \cdots, {f_{N/2 - 1}}\},$$\beta = \{ {\beta _{ - N/2}}, \cdots, {\beta _{N/2 - 1}}\},$$g = \{ {g_0}, \cdots, {g_{N - 1}}\},$$\gamma = \{ {\gamma _0}, \cdots, {\gamma _{N - 1}}\}$, and $h = \{ {h_0}, \cdots, {h_{N - 1}}\}$ are finite sequences of complex numbers.

We will derive three types of NUCFT regarding the sampling uniformity in the time domain and transformation domain. The description of the three problems is presented.

Problem 1: NUCFT(NUCFT-I): Uniform samples and non-integer frequencies:

We consider that given $\alpha$, to find $f = F(\alpha)$:

where $j = - N/2, \cdots, N/2 - 1$.

Problem 2: NUCFT(NUCFT-II): Nonuniform samples and integer frequencies:

We consider that given $\beta$, to find $g = G(\beta)$:

where $j = 0, \cdots, N - 1$.

Problem 3: NUCFT(NUCFT-III): Nonuniform samples and non-integer frequencies:

We consider that given $\gamma$, to find $h = H(\gamma)$:

The NUCFT algorithms for the above three problems are based on the following principal steps. Any function ${e^{i(cx + r{x^2})}}$ can be represented on any finite interval on the real line using a small number of the terms of ${e^{ - b{x^2}}}{e^{i(kx + r{x^2})}}$, and this number of the terms of $q$ is independent of the value $c$. For the efficient calculation (2.12)–(2.14), we approximate each ${e^{i({c_m}x + r{x^2})}}$ in terms of $q$-term CFS, and approximate the value of a CFS at each ${x_n}$ in terms of values at the nearest $q$ uniformly-spaced nodes.

3.2.   Derivation of algorithms

3.2.1.   Relevant facts from approximation theory

The main tool of this paper is to conduct a detailed analysis of CFS of $\phi :[- \pi, \pi] \to C$ given by the formula

where $b > \frac{1}{2}$ and $c$ are any real. We present this analysis in the lemmas and theorems of this subsection.

Lemma 3.1. For any real $b > \frac{1}{2}, c$ and any integer $k,$

Proof. Based on Lemma 2.4, we have

□.

Lemma 3.2. For any real $b > \frac{1}{2}, c$ and any integer $k,$

Proof.

Also,

Due to (3.8) and (3.9), we have

□.

Theorem 3.3. The function $\phi (x) = {e^{ - b{x^2}}}{e^{i(cx + r{x^2})}}$$(x \in (- \pi, \pi))$ can be approximated by the Chirp-Fourier series, and the approximation error can be obtained by the following inequality:

where $b > \frac{1}{2}$ and $c$ are constants

Proof. For $x \in (- \pi, \pi),$

The ${\sigma _k}$ is the $k$ th Chirp-Fourier coefficient for $\phi,$

where ${\rho _k} = \frac{1}{{2\sqrt {b\pi } }}{e^{ - {{(c - k)}^2}/4b}}, k = - \infty, \cdots, + \infty.$

Due to (3.13), we have

Due to (3.13), (3.14), and (3.15), we obtain

Owing to

and substituting (3.17) into (3.16), we have

Thus, combined with (3.18), we obtain

□.

According to Theorem 3.3, the functions ${e^{ - b{x^2}}}{e^{i(cx + r{x^2})}}$ can be approximated by the CFS whose coefficients are given analytically, and the error of the approximation decreases exponentially as $b$ increases. The coefficients ${\rho _k}$ have a peak at $k = [c]$($[\#]$ is the nearest integer to $\#$), and decay exponentially as $k \to \pm \infty$. We keep only the $q + 1$ largest coefficients, where the integer $q$ is chosen such as

and the following inequality is succeeded to

Theorem 3.4 is presented, which contains the truncation error of CFS of $\phi (x)$. Additionally, the functional approximation of ${e^{i(cx + r{x^2})}}$ is further derived.

Theorem 3.4. The function $\phi (x) = {e^{ - b{x^2}}}{e^{i(cx + r{x^2})}}(x \in (- \pi, \pi))$ can be approximated by the $q + 1$ term Chirp-Fourier series, and the truncation error can be obtained by the following inequality:

where $b > \frac{1}{2}$ and $c$ are constants and $q$ is an even integer such that $q \geqslant 4b\pi.$

Proof.

Owing to

Due to (3.24), (3.25), and Lemma 2.4, we have

and the combination of (3.20), (3.21), and (3.26), we obtain

Substituting (3.27) into (3.24) and (3.25), we have

Substituting (3.19) and (3.28) into (3.23), we obtain

Due to (3.22), we have

Corollary 3.5. For any integer $m \geqslant 2$, the conditions of Theorem 3.4 are satisfied, and we have

where $x \in [- \frac{\pi }{m}, \frac{\pi }{m}]$. □.

By extending further $x \in [- \frac{\pi }{m}, \frac{\pi }{m}]$ to $x \in [-d, d]$ in the Corollary 3.5, we obtain the following theorem:

Theorem 3.6. Let $b > \frac{1}{2}, c, d > 0$ be real numbers and $m \geqslant 2, q \geqslant 4b\pi$ be integers. Then, for any $x \in [-d, d],$ we have

where

□.

3.2.2.   Implementation of algorithms

We provide detailed algorithm descriptions in this subsection.

a) Problem 1: NUCFT(NUCFT-I): Uniform samples and non-integer frequencies

Setting $d = \pi, c = {x_j}, k = {\mu _k} + j$ in Theorem 3.6, we obtain that

where $m \in Z,$$k = 0, \cdots, N - 1, j = - q/2, \cdots, q/2 - 1, x \in [- \pi, \pi], \varepsilon = {e^{ - b{\pi ^2}(1 - 1/{m^2})}} \cdot (4b + 9),$ and ${P_{jk}}$ is defined as follows:

We will denote by ${\mu _k}$ the nearest integer to $m{\omega _k}.$

Setting $x = \frac{{2\pi l}}{N}, l \in Z,$ we have

so that

For $k = - \pi, \cdots, \pi$ and by $\{ {T_j}\}$, a set of complex numbers defined by the formula

for $j = - mN/2, \cdots, mN/2 - 1.$

Furthermore, ${\tilde f_j}$ will be denoted as

Combining (3.33)–(3.38), we obtain

where $j = - N/2, \cdots, N/2 - 1, \{ {f_j} = F{(\alpha)_j}\}.$

The implementation of NUCFT-I is presented as the following:

Step1: Input parameter $\{ {\alpha _0}, {\alpha _1}, \cdots, {\alpha _{N - 1}}\}, \{ {\omega _0}, {\omega _1}, \cdots, {\omega _{N - 1}}\},$ and then choose precision $\varepsilon, b$ and $q = [4b\pi];$

Step2: Compute ${\mu _k} = [m{\omega _k}], k = 0, 1, 2, \cdots, N - 1,$ and then calculate ${P_{jk}}$ and ${\tau _l}$ according to (3.34) and (3.36), respectively;

Step3: According to (3.37), the uniformly sampled Fourier series ${T_j}$ of ${\tau _k}$ in $[- \pi, \pi]$ is calculated using an inverse fast Fouier transform of length $mN, j = - mN/2, \cdots, mN/2 - 1.$

Step4: Select the $\{ {T_j}\}$ value at $j = - N/2, \cdots, N/2 - 1,$ and then get approximate values ${f_j}$ by calculating $\mathop {{f_j}}\limits^ \sim = {e^{(ir + b){{(\frac{{2\pi l}}{{mN}})}^2}}} \cdot {T_j}.$

b) Problem 2: NUCFT(NUCFT-II): Nonuniform samples and integer frequencies

Setting $d = N/2, c = {x_j}, k = {\nu _j} + k$ in Theorem 3.6, we have

where $m \in Z, j = 0, \cdots, N - 1, k \in [-N/2, N/2], \varepsilon = {e^{ - b{\pi ^2}(1 - 1/{m^2})}} \cdot (4b + 9),$ and ${Q_{jk}}$ is defined as follows:

We will denote by ${v_j}$ the nearest integer to ${x_j}mN/2\pi.$

Setting ${\omega _k} = n, n \in Z,$ we have

so that

where $j = - N/2, \cdots, N/2 - 1$.

For $n = - N/2, \cdots, N/2$, a set of complex numbers $\{ {U_l}\}$ is defined by the formula

where $l = - mN/2, \cdots, mN/2 - 1.$ Furthermore, ${\tilde g_j}$ will be denoted by the following formula:

Combining (3.40)–(3.45), we obtain

where $j = 0, \cdots, N - 1, \{ {g_j} = G{(\beta)_j}\}.$

The implementation of NUCFT-II is given as the following:

Step1: Input parameter $\{ {\beta _{ - N/2}}, \cdots, {\beta _{N/2 - 1}}\}, \{ {x_0}, {x_1}, \cdots, {x_{N - 1}}\},$ and then choose precision $\varepsilon, b$ and $q = [4b\pi];$

Step2: Compute ${v_j} = [{x_j}mN/2\pi], j = 0, \cdots, N - 1,$ and then calculate ${Q_{jk}}$ and ${\mu _k}$ according to (3.41) and (3.43), respectively;

Step3: According to (3.44), the uniformly sampled FT ${U_l}$ of ${u_k}$ in $[-N/2, N/2]$ is calculated using an inverse FFT of length $mN, l = - mN/2, \cdots, mN/2 - 1.$

Step4: Select the $\{ {U_l}\}$ value at $j = 0, \cdots, N - 1,$ and then get approximate values ${g_j}$ by calculating ${\tilde g_j} = e{}^{b{{(2\pi j/mN)}^2}} \cdot {U_{{v_j} + l}}.$

c) Problem 3: NUCFT(NUCFT-III): Nonuniform samples and non-integer frequencies

Setting $d = N/2, c = {t_j}/m, k = {v_j} + l, x = n$ in Theorem 3.6, we have

where $m \in Z, j = 0, \cdots, N - 1, k \in [-N/2, N/2], \varepsilon = {e^{ - b{\pi ^2}(1 - 1/{m^2})}} \cdot (4b + 9),$ and ${R_{jl}}$ is defined as follows:

We will denote by ${\eta _j}$ the nearest integer to ${x_j}N/2\pi.$

Due to (3.3) and (3.35), we obtain

so that

For $k = - N/2, \cdots, N/2$ and by $\{ {V_l}\}$, a set of complex numbers is defined by the formula

where $l = - {m^2}N/2, \cdots, {m^2}N/2 - 1.$

Furthermore, the ${\tilde h_j}$ will be denoted by the following formula:

Combining (3.47)–(3.52), we have

where $j = 0, \cdots, N - 1, \{ {h_j} = H{(\gamma)_j}\},$ and

The implementation of NUCFT-III is presented as the following:

Step1: Input parameter $\{ {\gamma _0}, \cdots, {\gamma _{N - 1}}\}, \{ {\omega _0}, {\omega _1}, \cdots {\omega _{N - 1}}\},$ and then choose precision $\varepsilon, b$ and $q = [4b\pi];$

Step2: Compute ${\eta _j} = [{x_j}N/2\pi], j = 0, \cdots, N - 1,$ and then calculate ${R_{jl}}$ and ${v_j}$ according to (3.48) and (3.50), respectively;

Step3: According to (3.51), the uniformly sampled Fourier series ${V_l}$ of ${v_k}$ in $[-N/2, N/2]$ is calculated using an inverse FFT of length ${m^2}N, l = - {m^2}N/2, \cdots, {m^2}N/2 - 1.$

Step4: Select the $\{ {V_l}\}$ value at $j = 0, \cdots, N - 1,$ and then get approximate values ${h_j}$ by calculating ${\tilde h_j} = e{}^{b{{({x_j}/m)}^2}} \cdot {e^{i\gamma {{({x_j}/m)}^2}}} \cdot \sum\limits_{l = - q/2}^{q/2 - 1} {{R_{jl}} \cdot {V_{{\eta _j} + l}}}.$

3.3.   Computation analysis of algorithms

Tables 1–3 display the comparison of the NUCFT computations with the direct computation for the three instances of uniform time but nonuniform transform domain, uniform transform domain but nonuniform time, and nonuniform time and transform domain, respectively.

The results show that the computational accuracy of the three NUCFT is similar, and the proposed algorithms have lower computational than direct calculation under the condition of guaranteeing the computational accuracy.

4.   Conclusions

In this paper, we have described three algorithms for computing DCFT for non-equispaced data, which is based on the interpolation formulae to transform function values from uniform to nonuniform points. The computational complexity of each algorithm is ${\text{O}}({\text{NlogN}} + {\text{Nlog(}}1/{\text{\varepsilon)}}),$ where $N$ is the data length and $\varepsilon$ is the computational accuracy. It shows that the derived approach is effective for computing nonuniform DCFT. The proposed algorithms can be viewed as generalizations of DCFT, and will have a broad range of applications in many branches of mathematics, science, and engineering. One of the specific applications is the inverse synthetic aperture radar (ISAR) imaging. In the signal extraction of ISAR imaging processing, the nonuniform rotation of the target can lead to nonuniformity of the processed data and also result in phase errors in the extracted signal. However, nonuniform CFT can process nonuniform data according to actual needs, so the algorithm proposed in this article can be used for signal extraction in ISAR to eliminate phase errors caused by nonuniform rotation. In future research, we will further investigate the application of the proposed algorithm.

Author contributions

Yannan Sun: Supervision, Resources, Conceptualization, Funding acquisition, Writing - review & editing. Wenchao Qian: Formal analysis, Methodology, Writing – original draft, Writing – review & editing. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors thank the editor and referees for their very helpful and detailed comments. This work is supported by the National Natural Science Foundation of China (No.62001193).

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.