Efficient computation of dot product of convolutions

Question

I have two FIR filters $f_{1},f_{2}$ and two short time-slices of signals $s_{1},s_{2}$. I need to compute: $$\left<s_{1}*f_{1},\ s_{2}*f_{2}\right>$$

Which I compute using six FFTs: each $s_{i}*f_{i}$ is computed by pointwise multiplication of the (padded) FFT vector of $s_{i}$ by the corresponding FFT of $f_{i}$, then an IFFT is computed and the result is truncated (this truncation is the reason computing the dot product in the frequency domain does not work).

The $f_{i}$ are known ahead of time, so this takes four FFTs at each iteration of my code.

Is there a way to get away with computing less? If the conjugate transpose of $f_{1}$ could be composed on $f_{2}$ and a FIR filter was the result, it would be posible to manage with just 2 FFTs with something like: $$IFFT\left(FFT(s_{2})\cdot FFT(f_{2}^T\circ f_{1})\right)$$

(where $FFT(f_{2}^T\circ f_{1})$ is computed ahead of time,) but this doesn't quite seem to work out.

A team member of mine asked a similar question (on another stackexchange community) but did not get helpful replies.

edit (details about truncation in my code):

To compute $\left<s*f\right>$ (supposing that $s$ is of length 20000 and $f$ of length 4000, though these numbers are made up) both vectors are zero-padded to length $20000+4000-1$, FFTs are multiplied pointwise, IFFT is computed on the product, and the samples in indices $\left[3999, 19998\right]$ are taken from the result (this represents a closed range of indices on a vector indexed from 0.)

This gives us all of the "meaningful" result computable from the slice of the signal: all dot products between $f$ and a shift of the signal-slice $s$ such that $f$ is shifted onto a completely-known part of $s$ (this is done over successive time slices of a signal as it is recorded.) It's important to notice the result is of length 16000 (the length of $s$ minus the length of $f$).

Answer 1

After thinking for a very long time on this problem, this is the best I've come up with:

Proposition. Given two FIRs $f_1,f_2$ of length $k+1$ and two signals $s_1,s_2$ of length $n+1$, we can compute the dot product $$<f_1*s_1,f_2*s_2>$$ using 4 FFTs of size $n+k+1$ (*) and 10 FFTs of size $2k$ (*). If the $f_i$ are given ahead of time, the computation is reduced to 2 FFTs of size $n+k+1$ (*) and 8 FFTs of size $2k$ (*)

(* or larger, i.e. FFT can be done on padded signals for better complexity.)

If the $f_i$ are given ahead of time, this method trades 2 FFTs of size $n+k+1$ for 8 FFTs of size $2k$. For the made up numbers $n\sim4000$ and $k\sim20000$ this might not be worth it. I'd say that if $n$ is larger than $15k$ then this method is probably worth it (i.e. $n+k>8\times2k$).

Proof. Set $\tilde{f_1},\tilde{f_2},\tilde{f_1},\tilde{s_2}$ to be the vectors $f_1, f_2, s_1, s_2$ zero-padded to length $N\ge n+k+1$. Denote by $\widehat{v}$ the DFT of a vector $v$, of the same dimensions as $v$ (though belonging to a vector space over $\mathbb{C}$). Also, denote by $v\times u$ the coordinate-wise multiplication of two vectors $u,v$. We know that

$$c_i:=\widehat{\widehat{\tilde{f_i}}\times\widehat{\tilde{s_i}}}$$

is the cyclic convolution of $f_i$ with $s_i$. The subvector from position $k+1$ to position $n+1$, is what we're calling $f_i*s_i$ in this post. Write $c_i$ as a concatenation of 3 vectors: $$c_i=(a_i|f_i*s_i|b_i)$$

so that we have:

\begin{align} <f_1*s_1,f_2*s_2>=&<c_1,c_2>-<a_1,a_2>-<b_1,b_2>\\ =&<\widehat{c_1},\widehat{c_2}>-<a_1,a_2>-<b_1,b_2>\\ =&<\widehat{\tilde{f_1}}\times\widehat{\tilde{s_1}}, \widehat{\tilde{f_2}}\times\widehat{\tilde{s_2}}>-<a_1,a_2>-<b_1,b_2>. \end{align}

The first summand is computed with 4 FFTs, or 2 FFTs if $\widehat{\tilde{f_1}},\widehat{\tilde{f_2}}$ are already given, of length $N$. The $a_i$ are the non-truncated application of $f_i$ to the first $k$ places of $s_i$. I.e. set $\bar{s_i}$ to be the vector of length $2k$ consisting of the first $k$ coordinates of $s_i$, and then $k$ zeros, and $\bar{f_i}$ to be the zero-padding of $f_i$ to length $2k$. Then the first $k$ positions of $$\widehat{\widehat{\bar{f_i}}\times\widehat{\bar{s_i}}}$$ are equal to $a_i$. Similarly, $b_i$ is the non-truncated application of $f_i$ to the last $k$ positions of $s_i$. Thus in order to compute $a_1,a_2,b_1,b_2$, we can use $10$ FFTs of size $2k$, and if we precompute $\widehat{\bar{f_i}}$ we can use 8 FFTs.

Note that in practice you might have the $s_i$ being part of a stream. In that case, the beginning of one $s_i$ is the end of a previous slice, so $a_i$ is the previous $b_i$, cutting the number of FFTs to 4 (in the case of known $f_i$). Now the method is already effective when $n>7k$. I'm guessing this works well for your intended application ;)

Answer 2

Convolution is a specific type of linear operator. So in all generality, you can write your expression as $\langle F_1 s_1 | F_2 s_2 \rangle$ with two linear operators $F_{1,2}$. Using the rules of operators in bilinear forms, we can rewrite this as $\langle s_1|F_1^\dagger F_2 |s_2\rangle$ and that again as $\langle s_1|\left(F_1^\dagger F_2\right) s_2\rangle$.

In your case, the two operators have a pleasant property. They are toeplitz, and if we take your padding into account, even circulant. The product of two circulant operators is again circulant, as is the conjugate transpose. That means you can write the result of the operator product as a single convolution. Your convolution kernel is the correlation product $\mathrm{corr}({f_1}^*,f_2)$.

If you implement this using fast convolution, you need to pad this correlation product and the size of your kernel will increase to the sum of the sizes of $f_{1,2}$. Accordingly, the required padding for $s_2$ will also increase.