I am wondering how to prove the below Fourier transform is non-negative? I did much simulation and it seems to be non-negative.
$$\int_0^\inf (be^{-at^p}-ae^{-bt^p})\cos(tx)dt, 0<a<b, \frac{1}{2}<p<1$$
In the OP, it assumed that $1/2<p<1$. Let us first show that, actually, the desired conclusion holds for $p\in(0,1/2]$. Let $$h(a):=\int_0^\infty \frac{e^{-at^p}}a\,\cos(tx)\,dt.$$ As noted by Johannes Hahn, it suffices to show that $h$ is decreasing.
We have
$$-h'(a)=\int_0^\infty g(t)\cos(tx)dt,$$
where
$$g(t):=\frac{e^{-a t^p} \left(a t^p+1\right)}{a^2},$$
and it is easy to see that $g$ is convex and decreasing to $0$ on $(0,\infty)$. So, by the Pólya criterion (see e.g. page 2310), $-h'(a)\ge0$, so that $h$ is indeed decreasing.
Now for any $p\in(0,1]$, the above reasoning shows that it suffices to prove that the function $s\mapsto (1+|s|^p)e^{-|s|^p}$ is positive definite.
It appears that the function $0\le u\mapsto (1+u^{1/2})e^{-u^{1/2}}$ is completely monotone on and hence a mixture of exponential functions $0\le u\mapsto e^{-cu}$ with $c>0$. So, the function $s\mapsto (1+|s|^p)e^{-|s|^p}$ is a mixture of functions $s\mapsto e^{-c|s|^{2p}}$ with $c>0$, which are positive definite for any $p\in(0,1]$, as the characteristic functions of (symmetric) stable_distributions.
So, it remains to check that the function $0\le u\mapsto (1+u^{1/2})e^{-u^{1/2}}$ is completely monotone on $[0,\infty)$.
Actually, the function $0\le u\mapsto (1+u^{1/2})e^{-u^{1/2}}$ is the Laplace transform of the positive function $0<t\mapsto \dfrac{e^{-1/(4t)}}{4 \sqrt{\pi }\, t^{5/2}}$ and hence a mixture of exponential functions $0\le u\mapsto e^{-cu}$ with $c>0$; see details on this below. So, we are done.
Proof of the statement on the Laplace transform: Note that the second derivative of $(1+u^{1/2})e^{-u^{1/2}}$ in $u$ is $e^{-u^{1/2}}/(4u^{1/2})$. So, after a simple rescaling, it is enough to show that $$J(a):=\int_0^\infty\exp\Big\{-\frac1t-at\Big\}\frac{dt}{2\sqrt t}=\frac{\sqrt\pi}2 \frac{e^{-2\sqrt a}}{\sqrt a}, \tag{1}$$ where $a>0$.
Using substitutions $t=u^2$ and then $u=1/(x\sqrt a)$, we get $$J(a)=\int_0^\infty\exp\Big\{-\frac1{u^2}-au^2\Big\}\,du =K(a)/\sqrt a,$$ where $$K(a):=\int_0^\infty\exp\Big\{-ax^2-\frac1{x^2}\Big\}\,\frac{dx}{x^2}.$$ Note that $K'(a)=-J(a)$ and $K(a)=J(a)\sqrt a$. So, we get the differential equation $$J'(a)=-\Big(\frac1{\sqrt a}+\frac1{2a}\Big)J(a),$$ whose general solution is given by $$J(a)=\frac c{\sqrt a}\,e^{-2\sqrt a}$$ for a constant $c$. To determine $c$, note that $$K(a)=J(a)\sqrt a =\int_0^\infty\exp\Big\{-\frac1{u^2}-au^2\Big\}\,du\,\sqrt a =\int_0^\infty\exp\Big\{-\frac a{y^2}-y^2\Big\}\,dy$$ and $$c=K(0+)=\int_0^\infty\exp\{-y^2\}\,dy=\frac{\sqrt\pi}2.$$ So, (1) follows. $\quad\Box$