In signal processing it is often that we need to compare two finite impulse responses (FIR). It is also often that we don’t care to know the difference between two FIRs outside a specific frequency band, and only care about the difference within that band. Given the difference FIR $latex p$, given as $latex n+1$ coefficients for a sample rate $latex N$, and the frequency band of interest being $latex [f_1, f_2]$, it is common to simply assume the sample rate is actually $latex n+1$ and define:
$latex power_{[\frac{f_1}{N}, \frac{f_2}{N}]}^{approx}(p):=\frac{1}{n+1}\sum_{\substack{f\in\mathbb{N} \\ \frac{f_1}{N}(n+1)\le f \le \frac{f_2}{N}(n+1)}} w(\frac{f}{n+1})|\hat{p}(\frac{f}{n+1})|^2$
where for $latex p=(p_0, …, p_n)$ we put $latex \hat{p}(z)=p_0+p_1e^{-2\pi\iota z}+…+p_n e^{-2\pi\iota nz}$, and we have some coefficients $latex w(\frac{f}{n+1})$.
We want to preserve Parseval’s Theorem, namely we would like the following to hold:
$latex power_{[\frac{0}{N}, \frac{N/2}{N}]}^{approx}(p)=p_0^2+…+p_n^2.$
A bit of algebra, that I’m not interested in showing here, gets us:
$latex
w(x) = \begin{cases}
1 & x=0\ or\ x= \frac{1}{2}\\
2 & otherwise
\end{cases}.
$
In order to compute $latex power_{[\frac{f_1}{N}, \frac{f_2}{N}]}^{approx}(p)$ we recognise that the numbers $latex \hat{p}(\frac{f}{n+1})$ are a subsequence of the discrete Fourier transform (DFT) of $latex p=(p_0, …, p_n)$. We compute the DFT of $latex p$, probably using an implementation of FFT, and sum as needed.
Let’s instead try compute the true power of $latex p$ in the band $latex [\frac{f_1}{N}, \frac{f_2}{N}]$. Instead of declaring the definition, I will first motivate it a bit. In the approximation above, we can double the sample rate we’ve assumed:
$latex power_{[\frac{f_1}{N}, \frac{f_2}{N}]}^{approx}(p;2):=\frac{1}{2n+2}\sum_{\substack{f\in\mathbb{N} \\ \frac{f_1}{N}(2n+2)\le f \le \frac{f_2}{N}(2n+2)}} w(\frac{f}{2n+2})|\hat{p}(\frac{f}{2n+2})|^2.$
The sum is twice as large so we get “more” frequencies involved. On the other hand, if we use the same computation method as above then we pay twice as much for the larger FFT. Now, instead of doubling, we can multiply the sample rate by a number $latex M$ to get the sum:
$latex power_{[\frac{f_1}{N}, \frac{f_2}{N}]}^{approx}(p;M):=\frac{1}{M(n+1)}\sum_{\substack{f\in\mathbb{N} \\ \frac{f_1M}{N}(n+1)\le f \le \frac{f_2M}{N}(n+1)}} w(\frac{f}{M(n+1)})|\hat{p}(\frac{f}{M(n+1)})|^2.$
The larger $latex M$ is, the more our sum should become independent of $latex M$ – we’re taking more and more frequencies. But now to compute this the FFT needed is pretty large. So instead of simply taking a large $latex M$, let’s take it to be even larger, and define:
$latex power_{[\frac{f_1}{N}, \frac{f_2}{N}]}(p):=\lim_{M\rightarrow\infty} \frac{1}{M(n+1)}\sum_{\substack{f\in\mathbb{N} \\ \frac{f_1M}{N}(n+1)\le f \le \frac{f_2M}{N}(n+1)}} w(\frac{f}{M(n+1)})|\hat{p}(\frac{f}{M(n+1)})|^2.$
Wait what? The FFT is big, so make it bigger??? Hold on: we can now recognise the limit on the right to be a Riemann integral. In fact, we have:
$latex power_{[\frac{f_1}{N}, \frac{f_2}{N}]}(p)=2\int_{\frac{f_1}{N}}^{\frac{f_2}{N}} |\hat{p}(f)|^2df.$
This is the true power of $latex p$ in the band $latex [\frac{f_1}{N}, \frac{f_2}{N}]$. It is independent of the sample rate, so from now on we’ll simply write, for two real numbers $latex f_1, f_2\in[0, \frac{1}{2}]$:
$latex power_{f_1, f_2}(p):=2\int_{f_1}^{f_2} |\hat{p}(f)|^2df.$
How do we compute this? Let’s apply some algebra.
We’ll perform some abuse of notation and use $latex p$ to also denote the polynomial $latex p(x):=p_0+p_1x+…+p_nx^n$. Putting the definition of $latex \hat{p}(f)$ into the definition of power we have:
$latex
\begin{aligned}
power_{f_1, f_2}(p)&=2\int_{f_1}^{f_2}|p(e^{-2\pi\iota f})|^2df \\
&=2\int_{f_1}^{f_2}p(e^{-2\pi\iota f})\overline{p(e^{-2\pi\iota f})}df \\
&=2\int_{f_1}^{f_2}p(e^{-2\pi\iota f})p(e^{2\pi\iota f})df
\end{aligned}
$
where we have used that, since $latex p_i$ are real, $latex \overline{p(x)}=p(\overline{x})$. We continue by setting:
$latex
\begin{aligned}
q(x):&=p(x^{-1})p(x)=\sum_{i=0}^{n} p_ix^{-i} \sum_{j=0}^{n} p_jx^j \\ &=\sum_{k=-n}^{n}\big(\sum_{\substack{0\le i,j\le n\\ i-j=k}}a_ia_j \big)x^k=\sum_{k=-n}^{n}q_kx^k
\end{aligned}
$
so that
$latex
\begin{aligned}
power_{f1,f2}(p) &=2\int_{f_1}^{f_2}q(e^{2\pi\iota f})df \\
&= 2\int_{f_1}^{f_2}\sum_{k=-n}^{n}q_ke^{2\pi\iota k f}df \\
&= 2\sum_{k=-n}^{n}q_k\int_{f_1}^{f_2}e^{2\pi\iota k f}df \\
&= 2\sum_{k=-n}^{n}q_k \frac{e^{2\pi\iota k f}}{2\pi\iota k}\biggr|_{f_1}^{f_2} \\
&= 2\sum_{k=-n}^{n}q_k \frac{e^{2\pi\iota k f_2} – e^{2\pi\iota k f_1}}{2\pi\iota k}.
\end{aligned}
$
With more abuse of notation, set $latex q=(q_{-n}, …, q_n)$. Also, set
$latex v_{f_1, f_2}^n := (2\frac{e^{2\pi\iota k f_2} – e^{2\pi\iota k f_1}}{2\pi\iota k})_{k=-n}^n.$
So we have:
$latex power_{f_1, f_2}(p)=(q, v_{f_1, f_2}^n)$
where $latex (\cdot, \cdot)$ is the usual dot product. Note that the coefficients of $latex q(x)$ are also the coefficients of $latex p(x)\cdot x^np(x^{-1})$, which is a multiplication of two degree $latex n$ polynomials and hence can be computed using two FFTs. This shows that the power of a polynomial in a given band can be computed using two FFTs of size $latex 2n+1$ (number of coefficients of $latex q$), some exponentiating, and a dot product. Even though we took our sample rate to infinity, the computation turned out pretty simple.
One more improvement before you go: instead of computing the coefficients of $latex q(x)$, we can use the identity
$latex (q, v_{f_1, f_2}^n)=(Uq, Uv_{f_1, f_2}^n)$
where $latex U$ is the unitary discrete Fourier transform operator. Since $latex q$ is a multiplication of polynomials, so a convolution, $latex Uq$ can be expressed in terms of $latex Up$:
$latex Uq=U(p(x)*x^np(x^{-1}))=Up\times\overline{Up}$
where the $latex \times$ means component wise. If we already have $latex Uv_{f_1, f_2}^n$ computed, this means that computing $latex power_{f_1, f_2}(p)$ takes only a single FFT of size $latex 2n+1$. Amazing!
Ok, that’s it.