Knowing a little C, goes a long way in Python

I’ve been branching out and learning some C while working on the latest release for Spectre. Specifically, I was migrating from a Python implementation of the short-time fast Fourier transform from Scipy, to a custom implementation using the FFTW C library (via the excellent pyfftw).

What I thought was quite cool was that doing the implementation first in C went a long way when writing the same in Python. In each case,

You fill up a buffer in memory with the values you want to transform. You tell FFTW to execute the DFT in-place on the buffer. You copy the DFT out of the buffer, into the spectrogram.

Understanding what the memory model looked like in C, meant it could basically be lift-and-shifted into Python. For the curious (and critical, do have mercy – it’s new to me), the core loop in C looks like (see here on GitHub):

for (size_t n = 0; n < num_spectrums; n++) { // Fill up the buffer, centering the window for the current frame. for (size_t m = 0; m < window_size; m++) { signal_index = m – window_midpoint + hop * n; if (signal_index >= 0 && signal_index < (int)signal->num_samples) { buffer->samples[m][0] = signal->samples[signal_index][0] * window->samples[m][0]; buffer->samples[m][1] = signal->samples[signal_index][1] * window->samples[m][1]; } else { buffer->samples[m][0] = 0.0; buffer->samples[m][1] = 0.0; } } // Compute the DFT in-place, to produce the spectrum. fftw_execute(p); // Copy the spectrum out the buffer into the spectrogram. memcpy(s.samples + n * window_size, buffer->samples, sizeof(fftw_complex) * window_size); }

The same loop in Python looks strikingly similar (see here on GitHub):

for n in range(num_spectrums): # Center the window for the current frame center = window_hop * n start = center – window_size // 2 stop = start + window_size # The window is fully inside the signal. if start >= 0 and stop <= signal_size: buffer[:] = signal[start:stop] * window # The window partially overlaps with the signal. else: # Zero the buffer and apply the window only to valid signal samples signal_indices = np.arange(start, stop) valid_mask = (signal_indices >= 0) & (signal_indices < signal_size) buffer[:] = 0.0 buffer[valid_mask] = signal[signal_indices[valid_mask]] * window[valid_mask] # Compute the DFT in-place, to produce the spectrum. fftw_obj.execute() // Copy the spectrum out the buffer into the spectrogram. dynamic_spectra[:, n] = np.abs(buffer)