Signal Processing
Functions for frequency-domain analysis and filtering. The FFT implementation uses Maxima’s fftpack5 library (a Lisp translation of FFTPACK), which supports arbitrary input lengths (most efficient when the length has the form 2^r * 3^s * 5^t).
np_fft (a) — Function
Compute the discrete Fourier transform of a 1D ndarray.
Scaling convention (matches NumPy):
Y[k] = sum(x[n] * exp(-2 * %pi * %i * k * n / N), n, 0, N-1)
The forward transform does not include a 1/N factor. Use np_ifft to invert.
The input a must be a 1D ndarray (real or complex). The result is always a complex 1D ndarray of the same length.
Examples
(%i1) /* FFT of a real signal */
A : ndarray([1.0, 2.0, 3.0, 4.0]);
(%o1) ndarray([4], DOUBLE-FLOAT)
(%i2) F : np_fft(A);
(%o2) ndarray([4], COMPLEX-DOUBLE-FLOAT)
(%i3) np_to_list(F);
(%o3) [10.0, 2.0 %i - 2.0, - 2.0, - 2.0 %i - 2.0]
(%i4) /* DC component = sum of signal */
np_ref(F, 0);
(%o4) 10.0
(%i5) /* Round-trip: ifft(fft(x)) = x */
np_to_list(np_real(np_ifft(F)));
(%o5) [1.0, 2.0, 3.0, 4.0]
Frequency analysis of a sine wave:
(%i1) N : 128$
(%i2) t : np_linspace(0, 1, N)$
(%i3) /* 10 Hz sine wave, sampled at 128 Hz */
x : np_sin(np_scale(2 * %pi * 10, t))$
(%i4) F : np_fft(x)$
(%i5) /* Magnitude spectrum */
mag : np_abs(F)$
See also: np_ifft, np_fft2d, np_abs, np_real, np_imag, np_angle
np_ifft (a) — Function
Compute the inverse discrete Fourier transform of a 1D ndarray.
Scaling convention (matches NumPy):
x[n] = (1/N) * sum(Y[k] * exp(+2 * %pi * %i * k * n / N), k, 0, N-1)
The inverse transform includes the 1/N normalization factor.
The input a must be a 1D ndarray (real or complex). The result is always a complex 1D ndarray of the same length. For real-valued signals, use np_real to extract the real part of the result.
Examples
(%i1) /* Frequency-domain filtering */
x : np_arange(0, 64)$
(%i2) F : np_fft(x)$
(%i3) /* Zero out high-frequency components (low-pass) */
(%i4) y : np_real(np_ifft(F))$
(%i5) /* Verify round-trip */
closeto(np_ref(y, 0), 0.0);
(%o5) true
See also: np_fft, np_ifft2d, np_real
np_fft2d (a) — Function
Compute the 2D discrete Fourier transform of a 2D ndarray.
Scaling convention (matches NumPy’s fft2):
Y[k1,k2] = sum_n1 sum_n2 x[n1,n2] * exp(-2 * %pi * %i * (k1*n1/N1 + k2*n2/N2))
The forward transform does not include a 1/(N1*N2) factor. Use np_ifft2d to invert.
The input a must be a 2D ndarray (real or complex). The result is always a complex 2D ndarray of the same shape. Internally uses separable 1D FFTs along columns then rows.
Examples
(%i1) A : ndarray(matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]));
(%o1) ndarray([3, 3], DOUBLE-FLOAT)
(%i2) F : np_fft2d(A);
(%o2) ndarray([3, 3], COMPLEX-DOUBLE-FLOAT)
(%i3) /* DC component = sum of all elements */
realpart(np_ref(F, 0, 0));
(%o3) 45.0
(%i4) /* Round-trip: ifft2d(fft2d(x)) = x */
B : np_ifft2d(F)$
(%i5) realpart(np_ref(B, 0, 0));
(%o5) 1.0
See also: np_ifft2d, np_fft
np_ifft2d (a) — Function
Compute the inverse 2D discrete Fourier transform of a 2D ndarray.
Scaling convention (matches NumPy’s ifft2):
x[n1,n2] = (1/(N1*N2)) * sum_k1 sum_k2 Y[k1,k2] * exp(+2 * %pi * %i * (k1*n1/N1 + k2*n2/N2))
The inverse transform includes the 1/(N1*N2) normalization factor.
The input a must be a 2D ndarray (real or complex). The result is always a complex 2D ndarray of the same shape. For real-valued data, use np_real to extract the real part.
Examples
(%i1) A : ndarray(matrix([1, 2], [3, 4]));
(%o1) ndarray([2, 2], DOUBLE-FLOAT)
(%i2) F : np_fft2d(A)$
(%i3) B : np_ifft2d(F)$
(%i4) /* Recovered values match original */
realpart(np_ref(B, 0, 0));
(%o4) 1.0
(%i5) realpart(np_ref(B, 1, 1));
(%o5) 4.0
See also: np_fft2d, np_ifft, np_real
np_convolve (a, b) — Function
Compute the 1D discrete linear convolution of two real signals.
Returns a 1D ndarray of length len(a) + len(b) - 1 (“full” mode, matching NumPy’s default). Both inputs must be 1D real ndarrays. Uses direct computation with O(len(a) * len(b)) complexity.
Examples
(%i1) /* Convolve with impulse: identity operation */
A : ndarray([1.0, 2.0, 3.0]);
(%o1) ndarray([3], DOUBLE-FLOAT)
(%i2) np_to_list(np_convolve(A, ndarray([1.0])));
(%o2) [1.0, 2.0, 3.0]
(%i3) /* Moving average (box filter) */
signal : ndarray([1.0, 3.0, 2.0, 5.0, 4.0]);
(%o3) ndarray([5], DOUBLE-FLOAT)
(%i4) kernel : np_scale(1/3, np_ones(3));
(%o4) ndarray([3], DOUBLE-FLOAT)
(%i5) np_to_list(np_convolve(signal, kernel));
(%o5) [0.333, 1.333, 2.0, 3.333, 3.666, 3.0, 1.333]
(%i6) /* Output length = len(signal) + len(kernel) - 1 */
np_shape(np_convolve(signal, kernel));
(%o6) [7]
See also: np_fft, np_ifft, np_convolve2d
np_convolve2d (a, kernel) — Function
Compute the 2D discrete convolution of a matrix with a kernel (“valid” mode).
Both inputs must be 2D real ndarrays. Returns a 2D ndarray of shape (h - kh + 1, w - kw + 1), where (h, w) is the input shape and (kh, kw) is the kernel shape. Only positions where the kernel fully overlaps the input are computed (no padding).
Uses direct O(h * w * kh * kw) computation with flat column-major indexing for performance.
Examples
(%i1) /* Identity convolution with a 1x1 kernel */
A : ndarray(matrix([1, 2, 3], [4, 5, 6]));
(%o1) ndarray([2, 3], DOUBLE-FLOAT)
(%i2) np_to_matrix(np_convolve2d(A, ndarray(matrix([1]))));
(%o2) matrix([1.0, 2.0, 3.0], [4.0, 5.0, 6.0])
(%i3) /* 3x3 box blur kernel */
blur : np_scale(1/9, np_ones([3, 3]));
(%o3) ndarray([3, 3], DOUBLE-FLOAT)
(%i4) img : np_ones([5, 5]);
(%o4) ndarray([5, 5], DOUBLE-FLOAT)
(%i5) result : np_convolve2d(img, blur);
(%o5) ndarray([3, 3], DOUBLE-FLOAT)
(%i6) /* Sobel horizontal edge detector */
sobel_x : ndarray(matrix([-1, 0, 1], [-2, 0, 2], [-1, 0, 1]));
(%o6) ndarray([3, 3], DOUBLE-FLOAT)
See also: np_convolve, np_fft