Implementation of Chirp-z transform

[Besler]

I want to see if there is any interest to include an implementation of czt and iczt in this package. I'm happy to do the work so long as we can agree on the interface and you are willing to provide review.

Theory
The DFT computes the spectral components along the unit circle. The Chirp-z transform (czt) computes components along spiral arcs, which are sometimes aligned along the unit circle. The DFT is a sampling along the y = j \omega straight line in the s-plane. The czt is a sampling of arbitrary straight lines in the s-plane.

My specific interest is in radar, where data from vector network analyzers is measured along the unit circle but not starting at the original. We measure data only from, say, 1 GHz to 5 GHz, and assume the components outside this range are all zeros. We can do this because we provide a band limited signal as input. The iczt method is the primary means of recovering the time domain signal.

Implementation
I generally see the czt and iczt implemented from Bluestein's algorithm. The czt breaks down into convolution, which can be implemented efficiently with fft, multiply, ifft. And fft and ifft are already implemented in this crate.

Design Considerations
I can think of three design considerations:

1. Support arbitrary czt
2. Support czt aligned on the unit circle
3. Support signals not sampled at equidistances

I think we can rule out 3) immediately and implement 1) with a simple interface for 2).

Interface
You probably have a way you want this interface to look. In general, I would like to see something along the lines of the following:

/// Compute the Chirp-z transform
fn czt(buffer: &mut [Complex<T>], a: &Complex<T>, w: &Complex<T>);

/// Compute the inverse Chirp-z transform
fn iczt(buffer: &mut [Complex<T>], a: &Complex<T>, w: &Complex<T>);

/// Compute the Chirp-z transform along the unit circle
///
/// Sampling starts at `theta_start` on the unit circle
fn unit_czt(buffer: &mut [Complex<T>], theta_start: &T);

/// Compute the inverse Chirp-z transform along the unit circle
///
/// Sampling starts at `theta_start` on the unit circle
fn unit_iczt(buffer: &mut [Complex<T>], theta_start: &T);

The fft, multiply, and ifft are performed inside czt and izct using FftPlanner. The samples variable are the z-plane spiral samples as defined in the algorithm. Internally, unit_czt and unit_iczt just create samples and call czt and iczt. czt and iczt are "in place" as is done currently.

At a later date, and if needed, the unit_czt and unit_iczt implementations can be specialized for performance. But generally, we'll see decent performance just by using the excellent implementations in the current crate. However, this interface doesn't follow the style of FftPlanner.

Literature

• Wikipedia
• Press, Mit. "A Linear Filtering Approach to the Computation of the Discrete Fourier Transform." (1969): 171-172.
• Shilling, Steve Alan. A study of the chirp Z-tranform and its applications. Diss. Kansas State University, 1972.
• Rabiner, Lawrence R., Ronald W. Schafer, and Charles M. Rader. "The chirp z‐transform algorithm and its application." Bell System Technical Journal 48.5 (1969): 1249-1292.
• Sukhoy, Vladimir, and Alexander Stoytchev. "Numerical error analysis of the ICZT algorithm for chirp contours on the unit circle." Scientific reports 10.1 (2020): 1-17.
• Sukhoy, Vladimir, and Alexander Stoytchev. "Generalizing the inverse FFT off the unit circle." Scientific reports 9.1 (2019): 1-12.
• MATLAB czt function

[Besler]

Should note that the fast iczt is patented (prior art date 2017) through the Iowa State University Research Foundation, Inc. If we want the fast version, we can ask if they plan to enforce and what their licensing policy is. However, we can do a "not fast" iczt implementation for now.

[ejmahler]

Thanks for bringing this up. I'm definitely interested, but give me a couple days to think over what kind of API to expose.

fn unit_czt(buffer: &mut [Complex<T>], start: &Complex<T>, end: &Complex<T>);

If this is operating on a unit circle, are start and end just angles?

Could we accomplish the spiral behavior by lerping from start to end in polar coordinates?

[Besler]

Apologies, I was a bit sloppy in the prototype. I edited above to agree with the MATLAB implementation and reduced the unit_czt to the correct number of parameters.

theta_start is exactly that, an angle. In the unit circle case, we actually already know the end angle because we would know the number of samples.

[ejmahler]

I came to FFTs from a practical direction, rather than a theoretical one, so unfortunately a lot of this goes over my head.

The wiki article emphaizes that bluestein's algorithm is closely related to the chirp-z transform. Do you know, in very high-level terms, what differences bluestein's algorithm has from the ctz?

[ejmahler]

Also: Are you aware of any reference implementations I can use for testing? They don't have to be in rust, any language is fine.

[Besler]

Bluestein's algorithm (for czt) implements czt using the fft. The equation drops into convolution, so you can implement the convolution as fft, multiply, ifft. It's just one method, though.

I personally haven't implemented czt or iczt before, but I'll give it a shot today.

Examples:

• czt python I know this guy and he uses this all the time so I would trust the implementation.
• bluestein in C
• rust czt

[Besler]

Okay, this is a by-the-books implementation of czt. I wasn't able to implement iczt as it doesn't have a nice form besides the Sukhoy-Stoytchev algorithm. I see most people implement iczt through a complement and czt algorithm.

use rustfft::{FftNum, num_complex::Complex};
use std::ops::{AddAssign};

/// Chirp Z transform
/// 
/// Implementation is O(n^2) and is simply looping over the expression.
fn czt<T>(buffer: &[Complex<T>], a: &Complex<T>, w: &Complex<T>) -> Vec<Complex<T>> where
    T: FftNum,
    Complex<T>: AddAssign
{
    let m = buffer.len();
    let mut result = vec![Complex::new(T::zero(), T::zero()); m];

    // CZT
    for k in 0..m {
        let z = a * w.powi(-1 * (k as i32));
        for n in 0..m {
            result[k] += buffer[n] * z.powi(-1 * (n as i32));
        }
    }

    result
}

fn main() {
    let t = 0f64;
    let f = 0f64;
    let df = 1_000f64;
    let dt = 1./df;

    println!("CZT");
    let a = Complex::new(0., 2.*std::f64::consts::PI * f * dt).exp();
    let w = Complex::new(0., -2.*std::f64::consts::PI * df * dt).exp();
    let input = vec![Complex::new(1., 0.), Complex::new(2., 0.), Complex::new(3., 0.), Complex::new(4., 0.)];

    let czt_result = czt(&input, &a, &w);
    println!("\tInput: {:?}\n\tCZT: {:?}", input, czt_result);
}

[HEnquist]

I have a practical question. Would thos benefit from being part of RustFFT, or could it just as well be a separate crate? For real-to-complex transforms I have made a crate (RealFFT) that presents an api similar to RustFFT, and internally it uses RustFFT for all the heavy work. Would that approach work here too?

[Besler]

It could definitely be a separate crate for the basic algorithms! There will be some extension, but that works.