I have a function $x(t)$, that is the periodic concatenation of truncated Lorentzians.
$$x(t) = \left(\sum_{n=-\infty}^\infty \delta(t-nT)\right) * \left(\frac{1}{t^2+a^2} \cdot \mathrm{rect}_T(t)\right)$$
The $*$ is the convolution operator, $\delta$ is the Dirac delta, and $\mathrm{rect}_T(t) = 1$ if $|t|\le T/2$ and 0 otherwise.
The factor on the right is a Lorentzian, truncated at a width $T$. Convolving it with the Dirac comb makes it periodic.
For $a=1$ and $T=5$, $x(t)$ looks like this:
I want to compute the Fourier transform. Easy enough:
- The convolution becomes a multiplication and vice versa;
- the comb becomes another comb $\mathcal F[\sum_{n=-\infty}^\infty\delta(t-nT)](f) = \frac{1}{T} \sum_{n=-\infty}^\infty\delta(f-n/T)$;
- the FT of the Lorentzian is a double sided exponential: $\mathcal F[\frac{1}{t^2+a^2}](f) \propto \exp(-2\pi a|f|)$;
- the FT of the rectangular window is of course a sinc function: $\mathcal F[\mathrm{rect}_T(t)](f) = T \mathrm{sinc}(fT)$.
Putting it all together, I find:
$$\mathcal F[x(t)] \propto \frac{1}{T} \sum_{n=-\infty}^\infty\delta(f-n/T)) \cdot \left(\exp(-2\pi a|f|) * T \mathrm{sinc}(fT)\right).$$
(I have left out some of the normalization factors that I would anyway have gotten wrong.)
The Dirac comb just turns the continuous Fourier transform into a discrete Fourier series. The Fourier series is then just a convolution of an exponential and a sinc. All makes sense to me.
I wanted to quickly double check this by running a numeric test. For some $a$ and $T$, I sample $x(t)$, and compute the FFT. AFAICT, the FFT should match my analytical derivation.
But it doesn't.
In particular, for high frequencies I expect that the FT decays asymptotically like the sinc (which decays much slower than the exponential). The (envelope of the) sinc decays like $1/f$, so that was what I was expecting. Neither solution reproduces this: the analytical evaluation is basically constant, while the FFT behaves like $1/f^2$.
Where am I wrong?
My test code:
import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import convolve
T = 10
a = 1
t = np.linspace(-T/2, T/2, 2**12)
yy = 1/(t**2 + a**2)
plt.figure()
plt.plot(t, yy)
plt.title('x(t)')
plt.show()
plt.close()
fft = np.fft.fftshift(np.fft.fft(yy))
ff = np.fft.fftshift(np.fft.fftfreq(yy.size, t[1] - t[0]))
ft_ana = convolve(np.pi/a * np.exp(-2*np.pi * a * np.abs(ff)), T*np.sinc(ff*T), 'same')
plt.figure()
plt.loglog(ff, np.abs(fft), label='FFT')
plt.loglog(ff, np.abs(ft_ana), label='analytical')
plt.legend(loc='best')
plt.show()
plt.close()
The output:
