2
$\begingroup$

The standard algorithm for discrete convolution of two vectors $x\in \mathbb{R}^{n}$ and $y \in \mathbb{R}^{m}$ is (in essence) a FFT of the two input vectors, multiplication of the two elementwise, and then an inverse FFT of the result. In practice there are complications, for example if $m \ne n$ then we need might need some padding of one of the inputs, and padding of the outputs.

From the following graph, you can see that the length of the zero padding affects the results, especially at the boundary. I have searched through the literature, but I could not find anything on the error analysis of the fast convolution algorithm (like Higham's analysis of the FFT in Accuracy and Stability). What error level should I expect for the fast convolution algorithm, and should the padding have this large of an affect?

enter image description here

For reference, the code to generate the figure is here:

#!/usr/bin/env python3
import numpy as np
from numpy import fft as np_fft
from matplotlib import pyplot as plt


def convolve_by_fft(a, b, len_after_padding=None):
    # scipy.fftconvolve doesn't allow control of the padding length,
    # so make this function mimics what is done in `scipy.fftconvolve`
    # i.e., ifft(fft(a) * fft(b))
    fft, ifft = np_fft.rfftn, np_fft.irfftn
    result_len = a.size + b.size - 1
    if len_after_padding is None:
        len_after_padding = result_len
    assert len_after_padding >= result_len
    sp1 = fft(a, (len_after_padding,))
    sp2 = fft(b, (len_after_padding,))
    ret = ifft(sp1 * sp2, (len_after_padding,))
    return ret[:result_len]


def next_power_of_2(x):
    return 1 << (x - 1).bit_length()


rng = np.random.default_rng(seed=10000)
a, b = rng.random((800,)), rng.random((800,))

result_direct = np.convolve(a, b)  # direct
result_fft_no_padding = convolve_by_fft(a, b)  # fft no padding
len1 = next_power_of_2(a.size + b.size - 1)
result_fft_padding_len1 = convolve_by_fft(a, b, len1)  # fft padding to next power of 2
len2 = 2 * len1
result_fft_padding_len2 = convolve_by_fft(
    a, b, len2
)  # fft padding to next next power of 2

# verify the error of last element
print(f"last element golden: {a[-1] * b[-1]}")
print(
    f"last element direct: {result_direct[-1]}, error = {(result_direct[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)
print(
    f"last element fft_no_padding: {result_fft_no_padding[-1]}, error = {(result_fft_no_padding[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)
print(
    f"last element fft_padding_len1: {result_fft_padding_len1[-1]}, error = {(result_fft_padding_len1[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)
print(
    f"last element fft_padding_len2: {result_fft_padding_len2[-1]}, error = {(result_fft_padding_len2[-1] - a[-1] * b[-1]) / (a[-1] * b[-1])}"
)

# plot the error compare to the direct method
plt.plot(
    (result_fft_no_padding - result_direct) / result_direct,
    label="error(fft_no_padding)",
)
plt.plot(
    (result_fft_padding_len1 - result_direct) / result_direct,
    label="error(fft_padding_len1)",
)
plt.plot(
    (result_fft_padding_len2 - result_direct) / result_direct,
    label="error(fft_padding_len2)",
)
plt.legend()
plt.show()
```
$\endgroup$
2
  • $\begingroup$ I think the error level is expectable, the zero frequency corresponds to the integral of the function, and with 1600 evaluations at machine precision $10^{-16}$ you'll end up in the order $10^{-13}$. $\endgroup$
    – davidhigh
    Aug 11 at 20:02
  • $\begingroup$ @davidhigh: Forgive me for lack of clarity; I agree that the error level is acceptable, but I still was interested in the error analysis . . . . $\endgroup$
    – user14717
    Aug 12 at 16:47

1 Answer 1

1
$\begingroup$

The fourier transform is only (strictly) applicable on periodic and continuous functions. Depending on your data in the two arrays, they might be a good approximation to a periodic signal when they start and end on similar values and each adjacent entry is not too divergent. Padding data arrays for FFT is essentially an ugly patch to "make it be transformable" (periodic). When you have strong jumps in your data array, then using the FFT algorithm as you described is a bad match. When you have discontinuous functions you can not properly express them in a continous base (sine/cosine). Any error estimate is then more ...academic. On the other hand, when your input arrays meet these conditions, then you do not need any padding at all.

fourier transform of rect

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.