FIR Filters and Convolution

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

miniDSP provides four related APIs so students can start with direct time-domain sums and then compare against the FFT overlap-add fast method.

---

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);

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);

---

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

Use direct methods to build intuition, then compare with overlap-add to see why FFT-based convolution matters for longer filters.

Further reading

• Convolution – Wikipedia
• Finite impulse response – Wikipedia
• Overlap-add method – Wikipedia

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