czt: Chirp Z-transform

In gjmvanboxtel/gsignal: Signal Processing

Chirp Z-transform

Description

Compute the Chirp Z-transform along a spiral contour on the z-plane.

Usage

czt(x, m = NROW(x), w = exp(complex(real = 0, imaginary = -2 * pi/m)), a = 1)

Arguments

x	input data, specified as a numeric vector or matrix. In case of a vector it represents a single signal; in case of a matrix each column is a signal.
m	transform length, specified as a positive integer scalar. Default: NROW(x).
w	ratio between spiral contour points in each step (i.e., radius increases exponentially, and angle increases linearly), specified as a complex scalar. Default: exp(0-1i * 2 * pi / m).
a	initial spiral contour point, specified as a complex scalar. Default: 1.

Details

The chirp Z-transform (CZT) is a generalization of the discrete Fourier transform (DFT). While the DFT samples the Z plane at uniformly-spaced points along the unit circle, the chirp Z-transform samples along spiral arcs in the Z-plane, corresponding to straight lines in the S plane. The DFT, real DFT, and zoom DFT can be calculated as special cases of the CZT[1]. For the specific case of the DFT, a = 0, m = NCOL(x), and w = 2 * pi / m[2, p. 656].

Value

Chirp Z-transform, returned as a vector or matrix.

Author(s)

Daniel Gunyan.
Conversion to R by Geert van Boxtel, [e-mail].

References

[1] https://en.wikipedia.org/wiki/Chirp_Z-transform
[2]Oppenheim, A.V., Schafer, R.W., and Buck, J.R. (1999). Discrete-Time Signal Processing, 2nd edition. Prentice-Hall.

Examples

fs <- 1000                                           # sampling frequency
secs <- 10                                           # number of seconds
t <- seq(0, secs, 1/fs)                              # time series
x <- sin(100 * 2 * pi * t) + runif(length(t))        # 100 Hz signal + noise
m <- 32                                              # n of points desired
f0 <- 75; f1 <- 175;                                 # desired freq range
w <- exp(-1i * 2 * pi * (f1 - f0) / ((m - 1) * fs))  # freq step of f1-f0/m
a <- exp(1i * 2 * pi * f0 / fs);                     # starting at freq f0
y <- czt(x, m, w, a)

# compare DFT and FFT
fs <- 1000
h <- as.numeric(fir1(100, 125/(fs / 2), type = "low"))
m <- 1024
y <- stats::fft(postpad(h, m))

f1 <- 75; f2 <- 175;
w <- exp(-1i * 2 * pi * (f2 - f1) / (m * fs))
a <- exp(1i * 2 * pi * f1 / fs)
z <- czt(h, m, w, a)

fn <- seq(0, m - 1, 1) / m
fy <- fs * fn
fz = (f2 - f1) * fn + f1
plot(fy, 10 * log10(abs(y)), type = "l", xlim = c(50, 200),
  xlab = "Frequency", ylab = "Magnitude (dB")
lines(fz, 10 * log10(abs(z)), col = "red")
legend("topright", legend = c("FFT", "CZT"), col=1:2, lty = 1)

gjmvanboxtel/gsignal documentation built on Nov. 22, 2023, 8:19 p.m.