Finite impulse response (FIR) filtering and convolution are core DSP tools: you shape a signal by summing delayed, weighted copies of it.

miniDSP provides five related APIs so users can start with direct time-domain sums, compare against the FFT overlap-add fast method, and design lowpass filters from scratch.

Time-domain convolution

For input \(x[n]\) (length \(N\)) and kernel \(h[k]\) (length \(M\)), full linear convolution is:

\[y[n] = \sum_{k=0}^{M-1} h[k]\,x[n-k], \quad n = 0,1,\ldots,N+M-2 \]

Out-of-range samples of \(x[\cdot]\) are treated as zero.

Reading the formula in C:

// n -> output index, k -> kernel index, x[n-k] -> signal[si], h[k] -> kernel[k]
for (unsigned n = 0; n < out_len; n++) {
    double acc = 0.0;
    for (unsigned k = 0; k < kernel_len; k++) {
        if (k > n) break;            // n-k would be negative
        unsigned si = n - k;
        if (si >= signal_len) continue;
        acc += signal[si] * kernel[k];
    }
    out[n] = acc;
}

API:

unsigned MD_convolution_num_samples(unsigned signal_len, unsigned kernel_len);
 
void MD_convolution_time(const double *signal, unsigned signal_len,
                         const double *kernel, unsigned kernel_len,
                         double *out);

Visuals — response and magnitude spectrum:

Quick example:

MD_convolution_time(signal, N, kernel, kernel_len, y_conv_time);

Moving-average FIR filter

The moving-average filter is the simplest FIR low-pass:

\[y[n] = \frac{1}{M}\sum_{k=0}^{M-1} x[n-k] \]

This implementation is causal and uses zero-padded startup (fixed divisor \(M\) at every sample).

Reading the formula in C:

// M -> window_len, x[n] -> signal[n], y[n] -> out[n]
double sum = 0.0;
for (unsigned n = 0; n < signal_len; n++) {
    sum += signal[n];
    if (n >= window_len)
        sum -= signal[n - window_len];
    out[n] = sum / (double)window_len;
}

API:

void MD_moving_average(const double *signal, unsigned signal_len,

unsigned window_len, double *out);

MD_moving_average

void MD_moving_average(const double *signal, unsigned signal_len, unsigned window_len, double *out)

Causal moving-average FIR filter with zero-padded startup.

Definition minidsp_fir.c:48

Visuals — response and magnitude spectrum:

Quick example:

MD_moving_average(signal, N, 8, y_ma);

General FIR filter

With arbitrary coefficients \(b[k]\) (length \(M\)), causal FIR filtering is:

\[y[n] = \sum_{k=0}^{M-1} b[k]\,x[n-k], \quad n = 0,1,\ldots,N-1 \]

This is the direct form used for most textbook FIR designs.

Reading the formula in C:

// b[k] -> coeffs[k], x[n-k] -> signal[n-k], y[n] -> out[n]
for (unsigned n = 0; n < signal_len; n++) {
    double acc = 0.0;
    for (unsigned k = 0; k < num_taps; k++) {
        if (k > n) break;    // zero-padded startup
        acc += coeffs[k] * signal[n - k];
    }
    out[n] = acc;
}

API:

void MD_fir_filter(const double *signal, unsigned signal_len,
                   const double *coeffs, unsigned num_taps,
                   double *out);

Visuals — response and magnitude spectrum:

Quick example:

MD_fir_filter(signal, N, coeffs, num_taps, y_fir);

FFT overlap-add fast convolution

For long filters, direct convolution is \(O(NM)\). Overlap-add computes the same full linear convolution in blocks using FFTs:

\[Y_b(k) = X_b(k)\,H(k), \quad y_b[n] = \mathrm{IFFT}(Y_b(k)) \]

Each block's time-domain result is added into the output with overlap between adjacent blocks, after an inverse FFT (IFFT) per block.

Reading the formula in C:

// x_b[n] -> block_time[n], H(k) -> (H_re[k], H_im[k]), y[n] -> out[n]
// Y_b(k) = X_b(k) * H(k), then overlap-add IFFT(Y_b) into out[]
for (unsigned start = 0; start < signal_len; start += block_len) {
    for (unsigned n = 0; n < nfft; n++) block_time[n] = 0.0;
    for (unsigned n = 0; n < block_len && start + n < signal_len; n++) {
        block_time[n] = signal[start + n];
    }
 
    // Educational DFT form for X_b(k) (FFT is used in production code).
    for (unsigned k = 0; k < nfft; k++) {
        double x_re = 0.0, x_im = 0.0;
        for (unsigned n = 0; n < nfft; n++) {
            double theta = 2.0 * M_PI * (double)k * (double)n / (double)nfft;
            x_re += block_time[n] * cos(theta);
            x_im -= block_time[n] * sin(theta);
        }
 
        // Complex multiply: Y_b(k) = X_b(k) * H(k)
        Y_re[k] = x_re * H_re[k] - x_im * H_im[k];
        Y_im[k] = x_re * H_im[k] + x_im * H_re[k];
    }
 
    // Educational inverse DFT + overlap-add.
    for (unsigned n = 0; n < nfft; n++) {
        double yb = 0.0;
        for (unsigned k = 0; k < nfft; k++) {
            double theta = 2.0 * M_PI * (double)k * (double)n / (double)nfft;
            yb += Y_re[k] * cos(theta) - Y_im[k] * sin(theta);
        }
        yb /= (double)nfft;
        if (start + n < out_len) out[start + n] += yb;
    }
}

API:

void MD_convolution_fft_ola(const double *signal, unsigned signal_len,
                            const double *kernel, unsigned kernel_len,
                            double *out);

Visuals — response and magnitude spectrum:

Quick example:

MD_convolution_fft_ola(signal, N, kernel, kernel_len, y_conv_fft);

Lowpass FIR filter design

The windowed-sinc method designs a lowpass FIR filter by sampling the ideal sinc impulse response and tapering it with a window function. Using a Kaiser window gives precise control over stopband attenuation via the \(\beta\) parameter.

For a normalized cutoff \(f_c = \mathrm{cutoff\_freq} / \mathrm{sample\_rate}\) and filter length \(N\):

\[h[i] = 2\,f_c\;\mathrm{sinc}\!\bigl(2\,f_c\,(i - (N-1)/2)\bigr) \;\cdot\; w_{\mathrm{Kaiser}}[i], \quad i = 0, 1, \ldots, N-1 \]

The coefficients are then normalized so they sum to 1.0 (unity DC gain).

Reading the formula in C:

// num_taps -> N, cutoff_freq -> Hz, sample_rate -> Hz
// fc -> normalized cutoff, center -> (N-1)/2, coeffs[i] -> h[i]
double fc = cutoff_freq / sample_rate;
double center = (double)(num_taps - 1) / 2.0;
 
double kaiser[num_taps];
MD_Gen_Kaiser_Win(kaiser, num_taps, kaiser_beta);
 
double sum = 0.0;
for (unsigned i = 0; i < num_taps; i++) {
    double x = (double)i - center;
    coeffs[i] = 2.0 * fc * sin(M_PI * 2.0 * fc * x)
              / (M_PI * 2.0 * fc * x) * kaiser[i];
    // (the sinc function handles x == 0 -> 1.0)
    sum += coeffs[i];
}
// Normalize for unity DC gain
for (unsigned i = 0; i < num_taps; i++) {
    coeffs[i] /= sum;
}

API:

void MD_design_lowpass_fir(double *coeffs, unsigned num_taps,
                           double cutoff_freq, double sample_rate,
                           double kaiser_beta);

Quick example:

double lpf[65];

MD_design_lowpass_fir(lpf, 65, 2000.0, sample_rate, 10.0);

Quick comparison

Method	Output length	Typical cost	Best use
Time-domain convolution	\(N+M-1\)	\(O(NM)\)	Teaching and short kernels
Moving average	\(N\)	\(O(N)\) (running sum)	First FIR low-pass example
General FIR	\(N\)	\(O(NM)\)	Arbitrary FIR taps
FFT overlap-add	\(N+M-1\)	\(O(N\log N)\) blocks	Long kernels / fast offline convolution
Lowpass FIR design	\(N\) taps	Design: \(O(N)\)	Anti-aliasing, band-limiting

Use direct methods to build intuition, then compare with overlap-add to see why FFT-based convolution matters for longer filters. Use the lowpass designer when you need a specific cutoff frequency with controlled stopband attenuation.

API reference

MD_convolution_num_samples() – compute full-convolution output length
MD_convolution_time() – direct time-domain full convolution
MD_moving_average() – causal moving-average FIR filter
MD_fir_filter() – causal arbitrary-coefficient FIR filter
MD_convolution_fft_ola() – full convolution via FFT overlap-add
MD_design_lowpass_fir() – design a Kaiser-windowed sinc lowpass FIR filter

Time-domain convolution

Moving-average FIR filter

General FIR filter

FFT overlap-add fast convolution

Lowpass FIR filter design

Quick comparison

Further reading

API reference