Open Besler opened 3 years ago
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.
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?
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.
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?
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.
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:
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);
}
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?
It could definitely be a separate crate for the basic algorithms! There will be some extension, but that works.
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:
The fft, multiply, and ifft are performed inside
czt
andizct
usingFftPlanner
. Thesamples
variable are the z-plane spiral samples as defined in the algorithm. Internally,unit_czt
andunit_iczt
just createsamples
and callczt
andiczt
.czt
andiczt
are "in place" as is done currently.At a later date, and if needed, the
unit_czt
andunit_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 ofFftPlanner
.Literature