The Fourier Transform of a Bessel function: $J_a(be^{-x}), a,b>0$ Currently we require the Fourier transform  $g(p)$ of  $f(x)=J_a(be^{-x}), a,b>0$
$$g(p)=\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}e^{ipx} f(x) dx$$
We could not find it in Tables and Handbooks. Mathematica gives it in terms of Gamma functions and Gauss Hyper geometric function. A closer look suggested that Hypergeometric function was actually a hoax. As we were finally interested in $|g(p)|^2$, it was a pleasant surprise to get
$$|g(p)|^2=\frac{1}{2\pi(p^2+a^2)},~~~~~(1)$$
which is independent of $b$ !, even the verification of (1) by
Mathematica in analytic mode  for various sets of $a,b,p$ has been time consuming.
The question is how to prove (1) and get $g(p)$ by hand?
EDIT:
By Integrate of  Mathematica, I finally get
$$g(p)=\frac{1}{2\sqrt{2\pi}} \left(\frac{2}{b}\right)^{ip} \frac{\Gamma[(a+ip)/2]}{\Gamma[1+(a-ip)/2]}$$
EDIT:
Here is the Mathematica out put, see the pre=factor of hypergeometric function is 0.
g(p)=1/2 Gamma[
1/2 (a + I p)] ((2^(I p) (1/b^2)^((I p)/2))/
Gamma[1/2 (2 + a - I p)] +
2^-a (-(1/b^2)^(-a/2) + b^
a) HypergeometricPFQRegularized[{1/2 (a + I p)}, {1 + a,
1/2 (2 + a + I p)}, -(b^2/4)])
 A: A Bessel function (of any kind) $f(z)$ satisfies the ODE
$$
w^2 f''(w) + w f'(w) + (w^2 - a^2) f(z) = 0.
$$
Making the variable transformation $w = be^{-x}$ leads to the following ODE satisfied by $f(x) = J_a \left(be^{-x}\right)$:
$$
f''(x) + \left(b^2 e^{-2x} - a^2 \right) f(x) = 0
$$
Now we can take the Fourier transform of this ODE. Assuming that the Fourier transform $g(p) = \mathcal{F}\left[ f(x) \right]$ can be analytically continued into the complex plane as
$$
g(z) = \int_{-\infty}^{\infty} e^{izx} f(x) \, dx,
$$
the Fourier transform of $e^{-2x} f(x)$ is
$$
\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{ipx} e^{-2x} f(x) \, dx
= \int_{-\infty}^{\infty} e^{i(p+2i)x} f(x) \, dx = g(p+2i).
$$
With this ODE is transformed into
$$
-p^2 g(p) + b^2 g(p+2i) - a^2 g(p) = 0
$$
or more concisely
$$
g(p+2i) = \frac{p^2 + a^2}{b^2} g(p).
$$
This recurrence can be solved with a little guidance from Mathematica. Letting $g(p) = h\left(\frac{p}{2i}\right)$ and $a = 2c$, the recurrence for $g$ implies the following for $h$:
$$
h(z+1) = 4 \frac{c^2-z^2}{b^2} h(z)
$$
Assuming separability of $h(z) = h_1(z) h_2(z) h_3(z)$ where each satisfies a recurrence $h_1(z+1) = \frac{4}{b^2} h_1(z)$, $h_2(z+1) = (c+z) h_2(z)$, and $h_3(z+1) = (c-z) h_3(z)$,
each of these has the solution, for non-negative $z$ and $c$ and sufficiently large $z$ (i.e. $z\geq c+1$),
$$
h_1(z+1) \propto \left( \frac{4}{b^2} \right)^z \\
h_2(z+1) \propto (z+c-1)! \\
h_3(z+1) \propto (-1)^z (z-c-1)!
$$
or for non-positive integral $z$ and $c$ the alternate possibilities
$$
h_2(z+1) \propto (-1)^z \frac{1}{(-c-z)!}, \quad
h_3(z+1) \propto \frac{1}{(c-z)!}
$$
This implies that there is no unique solution to this recurrence relation. Though technically this gives two possible solutions for integral ranges of $z$ and $c$, one using the first pair of $h_{2,3}$ solutions and the second using the second pair, the possible forms for $h_2$ and $h_3$ can be mixed and matched to still yield an expression that formally satisfies the recurrence. Any resulting such expression can be analytically continued to the complex plane. Edit: the possibilities for $h_{2,3}$, including a term $(-1)^z$ that alternates in sign across the integers, are somewhat anomalous and have less natural or unique analytic continuations to the plane.
As an example, take the very non-arbitrary choice of the first $h_2$ solution and second $h_3$ solution:
$$
h(z) \propto \left( \frac{4}{b^2} \right)^z \frac{\Gamma(c+z)}{\Gamma(c-z+1)}
$$
and hence
$$
g(p) \propto \left(\frac{b}{2} \right)^{ip} \frac{\Gamma\left(\frac{a-ip}{2}\right)}{\Gamma\left(1+\frac{a+ip}{2}\right)}
$$
The proportionality factor can be fixed using the fact that $g(0) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} J_a(be^{-x})\, dx = \int_0^{\infty} \frac{J_a(u)}{u} \, du$. This is independent of $b$ and equal to $\frac{1}{a\sqrt{2\pi}}$, yielding the final result
$$
g(p) = \frac{1}{2\sqrt{2\pi}} \left(\frac{b}{2} \right)^{ip} \frac{\Gamma\left(\frac{a-ip}{2}\right)}{\Gamma\left(1+\frac{a+ip}{2}\right)}
$$
Note that this is identical to the OP's proposed solution up to the substitution $p \mapsto -p$; see below comment on whether there is a possible sign disagreement. This result also satisfies the cross check that the $L^2$ norm of $g(p)$ seems to equal to that of $f(x)$, based on some numerical integrals for a few different choices of $a$ and $b$.
However, a few major limitations are that

*

*The solutions to the recurrence are not unique, and mixing and matching the solutions for $h_2$ and $h_3$ could lead to well-defined alternatives. I happened (after trial and error and multiple edits to the answer to streamline it) to pick the one of the four possibilities that corresponded to the answer provided by the OP (up to a sign change for $p$). I also can't guarantee that there aren't other functional forms that formally solve the recurrence.

*The original ODE has two free parameters that are fixed to the specify the Bessel $J$ functions, which must translate to two degrees of freedom in $g(p)$ -- but I can only identify one of them (the proportionality factor), and I'm curious where the other one went.

*The OP has obtained an actual result for the Fourier transform which disagrees with my results, but does, in fact, satisfy the recurrence
$$
g(p-2i) = \frac{p^2+a^2}{b^2} g(p).
$$
This leads me wonder if I have a sign error in the above work, but I can't find any. However, the rescaling argument made by @skbmoore does imply that the entire dependence of $g(p)$ on $b$ should be through a term like $b^{ip}$ (and not $b^{-ip}$), by making the change of variables $x = x' + \log(b)$ that maps $J_a(be^{-x})$ to $J_a(e^{-x'})$, and making this substitution in the Fourier transform creates a constant term $e^{ip\log(b)}$.

A: The Fourier transform has a closed-form, but my derivation lacks rigor.
Eliminate the b-dependence by translation: $x \to x+\log{(1/b)}.$ Thus
$$ g(p):= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{ipx}J_a(be^{-x})dx=$$
$$=\frac{(\log{(1/b)})^{ip}}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{ipx}J_a(e^{-x})dx = \frac{(\log{(1/b)})^{ip}}{\sqrt{2\pi}}G(a,p) $$
where $G(a,p)$ has been defined so that we can deal with a 2-variable integral now.
Use the known expansion (Gradshteyn & Rhyzik 8.411.10)
$$J_a(z) = \frac{(z/2)^a}{\Gamma(a+1/2)\Gamma(1/2)} \int_{-1}^1e^{izt}(1-t)^{a-1/2}dt$$
Insert this into $G(a,p).$ Interchange integrals.  Do the inner integral explicitly,
$$ \int_{-\infty}^\infty e^{ipx}e^{-ax}e^{ie^{-x}} = (-it)^{-a+ip}\Gamma(a-ip)$$
Here, though, the second non-rigorous step has been taken (the first is that the intergrals can be intergchanged).  Mathematica says the previous eq. is true for t with a non-zero imaginary part.  I'm assuming that someone good with contour integration can justify assuming completely real t is O.K.  Mathematica can do the remaining integral as well.  Collecting,
$$G(a,p):=\int_{-\infty}^\infty e^{ipx}J_a(e^{-x})dx=$$
$$=\frac{ 2^{-a} (-i)^{-a+ip}\ \Gamma(1/2+(ip-a)/2)}{2\Gamma(1/2)\Gamma(1+(ip+a)/2)}  (-1)^{-a}\big((-1)^{a}+e^{-p\pi})\big)\Gamma(a-ip).$$
(Hopefully there are no typos.)
As a numerical check, I let a=3/10 and p=6/10.  For the integral I chopped the integral limits at $\pm$ 1200, used 30 digits precision, and MinRecursion of 5 in Mathematica's NIntegrate.  I got 14 digits of agreement from the integral, and the explicit form presented.
Out of time; hope this helps.
