# Efficient computation of dot product of convolutions

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$).

• How exactly are you "truncating"? Can you show matlab, python or pseudo code? – Jazzmaniac Sep 6 '16 at 14:08
• @Jazzmaniac - edited the question to add an explanation. – Geva Yashfe Sep 6 '16 at 14:28
• Why do you pad with 2*4000? If you just wanted to avoid aliasing 1*4000 would have been sufficient. – Jazzmaniac Sep 6 '16 at 14:43
• You're right, this is a mistake in my description (the code handles the length of the result - in the described scenario 16000, in which case padding of 2*4000-1 is needed - and I had 20000 in mind instead). Fixed the question (including the truncated range)... – Geva Yashfe Sep 6 '16 at 14:57

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 $$$$ 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} =&--\\ =&<\widehat{c_1},\widehat{c_2}>--\\ =&<\widehat{\tilde{f_1}}\times\widehat{\tilde{s_1}}, \widehat{\tilde{f_2}}\times\widehat{\tilde{s_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 ;)

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.

• I'm not sure why the the transpose is circulant in this case. The filters are circulant in the sense the each row is a shift of the previous, but not all shifts are there because of the padding/truncation. – Geva Yashfe Sep 7 '16 at 6:15
• If you have padding, you can make the matrix that acts on the padded signal circulant because it will hit zeros anyway. – Jazzmaniac Sep 7 '16 at 12:42
• Can you give pseudocode? As described in the original question, I've already tried this idea and I doubt it works in this case. The padding and truncation are problematic and the relevant product operator is not circulant. You can't make it circulant by adding rows or columns at the beginning or end. – Geva Yashfe Sep 8 '16 at 6:24
• You don't pad the operator. You pad the signal. If the signal is zero, then it doesn't matter how the operator operates on that part of the signal. With that, you can make the operator circulant and it does to your specific padded signal what it is supposed to do. And I'm sorry to say that, but I'm too busy to write code now. If you require more of my time, you can always hire me. My answer takes you nearly all the way, so I doubt that is necessary if you try. – Jazzmaniac Sep 8 '16 at 16:40
• Of course you need to change the size of the matrix representations if you pad the signal. And it is you who is condescending. You demand help, I give it to you. You say I'm wrong without providing evidence and then attempt to force me into helping you more by accusing me of not helping you. It's your position to provide the code that you say is not working and to ask specifically if you don't understand my answer. And if you want professional help then better be prepared to pay for it. I'm done with you. My answer stands and is correct, if you have trouble implementing it, don't blame me. – Jazzmaniac Sep 8 '16 at 18:18