# Period-based Fast Fourier transform

I have a list of events times $t_n$ in which I would like to find repeating patterns in the temporal density (how many events per second). And finally plot how likely it is to find a pattern of period $p$ to let the user decide what's worth considering.

Therefore I thought I would take a continuous approximation of the density function of my event dates with a Kernel Density Estimate (KDE). Then take its Fourier transform, and finally discretize it by period instead of frequency. (Longer backstory here. And this might also solve another issue I've had.)

So, here is the KDE function and it's Fourier transform if I'm not wrong. $$k_h(t) = \frac{1}{Nh\sqrt{2\pi}} \sum_{n=0}^{N-1} \exp\left(-\frac{1}{2}\left(\frac{t - t_n}{h}\right)^2\right)$$ $$\widehat{k_h}(\nu) = \frac{e^{-2 \pi^2 h^2 \nu^2}}{N} \sum_{n=0}^{N-1} e^{-2 i \pi t_n \nu}$$

And since I want to plot it in terms of period instead of frequency, I can define $l_h$ taking an integer argument (the period length) as follow.

$$l_h(\lambda) = \widehat{k_h}\left(\frac{1}{\lambda}\right) = \frac{e^{\frac{-2 \pi^2 h^2}{\lambda^2}}}{N} \sum_{n=0}^{N-1} \exp\left(\frac{-2 i \pi t_n}{\lambda}\right)$$

This actually seems to work. But it's slow for a large $N$. So I was wondering if this could be accelerated in a way similar to the FFT. I can quantify the $t_n$ by rounding them to the nearest multiple of $\Delta_t$ and define $x_n$ to be the number of $t_k$ between $n\Delta_t$ and $(n+1)\Delta_t$ and approximate $l_h$ by $l'_h$ like this.

$$l'_h(\lambda) = \frac{e^{\frac{-2 \pi^2 h^2}{\lambda^2}}}{N} \sum_{n=0}^{N'-1} x_n \exp\left(\frac{-2 i \pi n}{\lambda}\right)$$

Note that $\Delta_t$ doesn't appear in the equation because I can take it equal to 1 (e.g. 1 hour) and decide the argument of $l'_h$ is expressed in multiples of that unit ($\lambda$ hours).

If it weren't for the $\frac{1}{\lambda}$ in the complex exponential, it would definitely look like the definition of a discrete Fourier transform. I tried to apply the same technique as the Cooley-Tukey FFT algorithm but got stuck pretty fast when trying to find an periodicity in the complex exponential of the even indexed input. (Similar to $E_{k+\frac{N}{2}} = E_k$ on wikipedia.) The exponential is just not periodic in the period length.

I'm open to any suggestion, doesn't matter if $N'$ has to be prime or a power-of-two or anything else. I can always make up for that.

• Potentially relevant: Chirp $z$-transform. – ccorn Sep 30 '16 at 18:22
• I didn't know any of this before, it's a lot to take in, I may have missed something. But I'm not sure how this could help. The fact that I divide by $p$ means it's not even a Z-transform. And the Bluestein algorithm cannot be applied because it wouldn't transform to a convolution. (BTW, wikipedia talk more about the Bluestein algorithm than the Chrip Z-transform.) – Celelibi Oct 1 '16 at 23:34
• I see. No, that does not seem to help, sorry. – ccorn Oct 2 '16 at 20:00
• You might want to look at the autocorrelation function instead: FFT, squared magnitude thereof (spectral power density), then inverse FFT. Of interest would be local extrema at some nonzero times. – ccorn Oct 2 '16 at 20:17
• Well, I tried to use the autocorrelation (see the two linked questions in the begining of this question). But unfortunately the result is not satisfactory. It doesn't show sharp peak at the correct period (slow to rise and slow to fade) and repeated patterns show multiple local maxima. – Celelibi Oct 4 '16 at 1:44