7

I am looking for a efficient way to find the value of k that minimizes

$\sum(s_t - b_{t+k})^2$

where $s$ and $b$ are N-dimensional vectors and the values are wrapped around like this:

$b_{t+k} := b_{(t+k)-N}$ if $(t+k) > N$

$b_{t+k} := b_{(t+k)+N}$ if $(t+k) < 1$

I need this for some algorithm and at the moment I do not have any better idea than to try all possible values of $k$.

463035818_is_not_a_number
  • 575

2 Answers2

8

First, note that if you multiply out the square, the two sums of squared terms don't depend on $k$, so you're effectively maximizing $\sum_ts_tb_{t+k}$. If you read up a bit on convolutions and discrete Fourier transforms, you'll find that they allow you to find that maximum in $N\log N$ time instead of the $N^2$ time it takes you to try out all values of $k$.

Edit with a rough sketch of what to do: Perform a Fourier transform on both sequences, multiply the one transform with the complex conjugate of the other (conjugating the transform corresponds to inverting the original sequence in order to write your sum as a convolution), then perform an inverse Fourier transform on the product. That gives you your cross-correlation sum as a function of $k$, and you can read off the optimal $k$ from the value with the highest magnitude.

joriki
  • 238,052
  • I didnt want to make the question too long... Once I noticed that the only part that depends on k is the scalar product between the two vectors I thought rotations will help. I will have to do some reading on convolutions and fourier, because atm I dont see how this would work. – 463035818_is_not_a_number Jul 07 '15 at 12:20
  • The point is that a cyclic convolution will give you the value of this scalar product for all possible shifts, so the maximization is $O(n)$ after that. Using Fourier methods you can compute this full convolution faster. – Michael Grant Jul 07 '15 at 12:24
  • It is also possible to try and align the phase while still in the Fourier domain. A translation in the spatial domain is equivalent to adding a phase $k\omega$ increasing linearly with frequency in the Fourier domain, which can then be used to solve a linear equation system for the phase. – mathreadler Jul 07 '15 at 14:11
  • @mathreadler: It's a good point, but I don't see how it results in a linear system. – joriki Jul 07 '15 at 14:12
  • 1
    Ok, I will just make a minimal example to convince myself it actually works, then I can show it in a separate answer, ok? ;) – mathreadler Jul 07 '15 at 14:14
1

Another approach is to use the translation or time-shift theorem in Fourier analysis:

$$h(x) = f(x-x_0) \Leftrightarrow \mathcal{F}\{h\}(\omega) = e^{i2\pi x_0 \omega} \mathcal{F}\{f\}(\omega)$$ We see that the phase increase is a linear function of $\omega$.

  1. Do a one way FFT for each signal,
  2. Calculate the phases.
  3. Calculate the derivative w.r.t. $\omega$ of the difference of the phases.
  4. Calculate a certainty measure from the absolute values.
  5. Calculate weighted mean of step 3 using certainty measure above.

The certainty measure is to avoid noise or numerical error to influence the solution. Parts of the Fourier domain which have low absolute value are much more noise sensitive, therefore the need to have a certainty measure.

Here is a sample 3rd order polynomial being displaced, noised and then corrected: enter image description here Here is phase estimate in blue and certainty measure in red. Y axis is estimated displacement for blue and x is frequency with DC being in the middle. enter image description here

mathreadler
  • 25,824
  • I think this answer is missing a lot of details. What do you take the derivative with respect to? How do you calculate a "certainty measure"? What does the weighted mean give you? –  Jul 07 '15 at 15:19
  • The derivative is with respect to the frequency variable in the Fourier domain. The certainty I chose to calculate as $\sqrt{|\mathcal{F}{h}(w)|}$, but any monotonically increasing function of the absolute value is a candidate - can depend lots on the application. The weighted mean is to trust parts of the fourier domain where we know that the signals reside. Parts of the fourier domain with low absolute value are more noise sensitive, therefore we should not trust them as much. – mathreadler Jul 07 '15 at 15:39
  • I find this very interesting but I will need some time to reproduce and fully understand it. – 463035818_is_not_a_number Jul 08 '15 at 09:21
  • To get it to be more robust to small shifts we need to perform the weighted averaging on a circle, but I guess that is over-course. – mathreadler Jul 08 '15 at 19:13