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

Solution 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.

Solution 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.