Chirp Z-transform
Compute the Chirp Z-transform along a spiral contour on the z-plane.
czt(x, m = NROW(x), w = exp(complex(real = 0, imaginary = -2 * pi/m)), a = 1)
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:
|
w |
ratio between spiral contour points in each step (i.e., radius
increases exponentially, and angle increases linearly), specified as a
complex scalar. Default: |
a |
initial spiral contour point, specified as a complex scalar. Default: 1. |
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].
Chirp Z-transform, returned as a vector or matrix.
Daniel Gunyan.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.
[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.
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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.