Inverse Fourier transform of $\text{sinc}(t)^{1/k}$

Question

This question is a follow-up to a statistics question on crossvalidated.SE where we are looking to see if it is possible for the difference of two independent random variables to be uniformly distributed from negative to positive unity. We have discovered that a necessary condition for this outcome is to have a distribution with a Fourier transform with squared-norm equal to a reciprocal power of the sinc function (the Fourier transform for the uniform distribution). The present question is seeking to find out if this is possible or not.

So, with that in mind, I would like to know if there is any probability density with a characteristic function that is a reciprocal power of the sinc function. For a density $f$ with this property, the Fourier transform $\hat{f}$ is the function:

$$\hat{f}(t) = \text{sinc}(t)^{1/k} \quad \quad \quad \text{for } k \in \mathbb{N}.$$

For $k=1$ the inverse Fourier transform $f$ is known to be a rectangular pulse function (i.e., a uniform density). I am mostly interested in the result when $k=2$, but the general result is also of interest to me.

Update: It has been pointed out in the comments that $z^{1/k}$ is a multi-valued function, and I have not specified which branch I am using. Given the above motivation for the problem, unless I am mistaken (and I might be, so please tell me if I am), it should not matter which branch is used.

Answer 1

This is not an answer to your question but rather a comment to your post on crossvalidated.SE.

In your post, you considered the case $k=2$. The question you tried to answer is if there are two i.i.d. distributions $A,B$ such that $A-B$ is uniformly distributed. It turns out that this is not possible. The condition for the characteristic function $\phi$ to be satisfied reads $$|\phi(t)|^2 = \operatorname{sinc}(\pi t)$$ which clearly has no solution.

A different question is if $A+B$ can be uniformly distributed. For this we need that $$\phi(t)^2 = \operatorname{sinc}(\pi t).$$ The solution is $$ \phi(t) = \sqrt{\operatorname{sinc}(\pi t)}$$ with a complex $\sqrt{\cdot}$ (whose branch cut is away from the real axis). However, in order that $\phi(t)$ is a characteristic function of a probability distribution, it has to fulfill (at least) $\phi(-t) = \overline{\phi(t)}$. Let us check $$ \phi(-t) = \sqrt{-\operatorname{sinc}(\pi t)} = \pm i \sqrt{\operatorname{sinc}(\pi t)} \neq \overline{\phi(t)}$$ for any choice of the square-root.

So it is not possible to have a uniform distribution as a sum of $k=2$ i.i.d. distributions.

Answer 2

As @adfriedman points out, your original problem actually has no solution in $\mathbb{C}$. Note that

$$z\overline{z} = |z|^2$$

is always real and non-negative. Your original problem

$$\varphi(t) \overline{\varphi(t)} = |\varphi(t)|^2 = \text{sinc}(t)$$

doesn't have any solutions, as $\mathrm{sinc}(t)$ has negative excursions. (Maybe if you add $0.217234..$ to the right hand side, you can find a solution.)

However, trying to find the Inverse Fourier Transform of $(\mathrm{sinc}(\omega))^{\frac{1}{2}}$ is interesting in its own right, so I took a look at what it might be.

As @Alex_Francisco notes, $(\mathrm{sinc}(\omega))^{\frac{1}{2}}$ is a complex, multi-valued function. Since $\mathrm{sinc}(\omega)$ is real and only the root of a negative real value introduces ambiguity, below I'll use the branch of $-1 = e^{i\pi + i2\pi n}$ for $n = 0$, so that $(-1)^{\frac{1}{2}} = e^{i\frac{\pi}{2}} =i$ unambiguously.

I briefly tried to figure out the Inverse Fourier Transform of $(\mathrm{sinc}(\omega))^{\frac{1}{2}}$ analytically. However, without even having a notion of what the solution should look like, that path was too difficult to get any useful results.

So, I resulted to using a (windowed) Inverse Discrete Fourier Transform to get a feel for what the solution should look like numerically, even if the analytic solution was not possible.

The following is my Octave (an open source MatLab clone) script I used:

N = 32768*8;
Ts = 0.0001; % sampling period
Fs = 1.0/Ts; % sample rate
f = [-(N/2):(N/2-1)]*Fs/N; % frequency axis
w = 2*pi*f; % radian frequency axis
t = [(-N/2):(N/2-1)]*Ts; % time axis

pkg load signal;  % An Octave-ism.  Not sure if needed by actual MatLab
W = blackmanharris(N, "periodic").';
%W = rectwin(N).';
% FIXME: compute window gain (aka window attenuation)

s = sin(w)./(w); % sinc(w) = sin(w)/w
s(N/2+1) = 1.0; % remove the singularity

% Window before inverse transforming to reduce aliasing effects
s = s .* W;
% FIXME: also scale for window gain (aka window attenuation)

figure(1);
plot (w, real(s), w, imag(s));
xlabel('Radian Frequency, \omega (rad/s)');
ylabel('Amplitude');
title('Windowed sinc(\omega) (real and imaginary parts)');

% Circularly shift the frequency domain samples, so when we take the DFT, we
% have a sinc() properly centered at 0 radians/sec.
% Take the (rectangularly windowed) IDFT and then circularly shift the time
% samples for nicer display on plots.
S = fftshift(ifft(ifftshift(s),N));
Smag = abs(S);
SmagdB = 10*log10(Smag);
Sphase = arg(S);

if 0
figure(2);
plot(t, Smag);
xlabel('Time, t (seconds)');
ylabel('Magnitude');
title('Linear Magnitude of IDFT of Windowed sinc(\omega)');

figure(3);
plot(t, SmagdB);
xlabel('Time, t (seconds)');
ylabel('Magnitude (dB)');
title('Logarithmic Magnitude of IDFT of Windowed sinc(\omega)');

figure(4);
plot(t, Sphase);
xlabel('Time, t (seconds)');
ylabel('Phase (radians)');
title('Phase of IDFT of Windowed sinc(\omega)');
endif

figure(5);
plot(t, real(S), t, imag(S));
xlabel('Time, t (seconds)');
ylabel('Amplitude');
title('IDFT of Windowed sinc(\omega) (real and imaginary parts)');

s2 = sqrt(s); % sqrt(sinc(w))

% Window before inverse transforming to reduce aliasing effects
s2 = s2 .* W;
% FIXME: also scale for window gain (aka window attenuation)

figure(6);
plot(w, real(s2), w, imag(s2));
xlabel('Radian Frequency, \omega (rad/s)');
ylabel('Amplitude');
title('Windowed (sinc(\omega))^{(1/2)} (real and imaginary parts)');


% Circularly shift the frequency domain samples, so when we take the DFT, we
% have a sinc() properly centered at 0 radians/sec.
% Take the (rectangularly windowed) IDFT and then circularly shift the time
% samples for nicer display on plots.
S2 = fftshift(ifft(ifftshift(s2),N));
S2mag = abs(S2);
S2magdB = 10*log10(S2mag);
S2phase = arg(S2);

if 0
figure(7);
plot(t, S2mag);
xlabel('Time, t (seconds)');
ylabel('Magnitude');
title('Linear Magnitude of IDFT of Windowed (sinc(\omega))^{(1/2)}');

figure(8);
plot(t, S2magdB);
xlabel('Time, t (seconds)');
ylabel('Magnitude (dB)');
title('Logarithmic Magnitude of IDFT of Windowed (sinc(\omega))^{(1/2)}');

figure(9);
plot(t, S2phase);
xlabel('Time, t (seconds)');
ylabel('Phase (radians)');
title('Phase of IDFT of Windowed (sinc(\omega))^{(1/2)}');
endif

figure(10);
plot(t, real(S2), t, imag(S2));
xlabel('Time, t (seconds)');
ylabel('Amplitude');
title('IDFT of Windowed (sinc(\omega))^{(1/2)} (real and imaginary parts)');

I used windowing to reduce the effect of aliasing on the IDFT, to get a better picture of what the analytic answer should look like. There are some artifacts due to the truncation of the infinitely long function that still remain though.

This figure shows a plot of the real and imaginary parts of the truncated, windowed $(\mathrm{sinc}(\omega))^{\frac{1}{2}}$ function:

This is that same figure zoomed in on the center:

which shows an alternation of the real and imaginary parts due to the positive and negative excursions of $\mathrm{sinc}(\omega)$ (before the square root was taken).

The next few figures plots of the real and imaginary parts of the IDFT of the truncated, windowed $(\mathrm{sinc}(\omega))^{\frac{1}{2}}$ function, with various regions zoomed in:

From these plots, it is fairly reasonable to conjecture that the Inverse Fourier Transform of $(\mathrm{sinc}(\omega))^{\frac{1}{2}}$ is an infinite sum of shifted, weighted $\delta(t)$ and $\delta'(t)$ impulses. Something like

$$\begin{align*}f(t) &= (a_0 + ib_0)\delta(t) \\ &+ \sum_{n=-\infty}^{-1} (-a_{2n+1} +i b_{2n+1})\delta'(t - [2n+1]) + (-a_{2n}-ib_{2n})\delta(t -2n) \\ &+ \sum_{n=1}^{\infty} (a_{2n-1} -i b_{2n-1})\delta'(t - [2n-1]) + (-a_{2n}-ib_{2n})\delta(t -2n) \\ \\ &= (a_0 + ib_0)\delta(t) \\ &+ \sum_{n=1}^{\infty} -\left(a_{2n} + ib_{2n}\right)\left[\delta(t+2n)+\delta(t-2n)\right] \\ &+ \sum_{n=1}^{\infty} -\left(a_{2n-1} - ib_{2n-1}\right)\left[\delta'\left(t+\left[2n-1\right]\right)-\delta'\left(t-\left[2n-1\right]\right)\right] \\ \end{align*}$$

with $|a_m| = |a_{-m}|$ and $|b_m| = |b_{-m}|$ and likely $|a_m| = |b_m|$.

I didn't bother to try and do a curve fit to figure out the envelope and the nominal values of the $a_m$'s and the $b_m$'s, especially since my windowing scales the values. I suspect they follow an envelope that is a power of a $\mathrm{sinc}()$ function, but that is pure speculation.