PDF of $ | \sum_{k=1}^{n}{|h_k||g_k|\exp\left( j \theta_k \right)} |^2$ for small values of $n$ and $Q$? Given the following function of random variables
$$f = \left|\sum_{k=1}^{n}{|h_k||g_k|\exp\left( j \theta_k \right)} \right|^2,$$
where $h_1, \cdots, h_n$ and $g_1, \cdots, h_n$ are i.i.d. random variables following the complex Gaussian distribution $\mathcal{CN}(0, \sigma_{h}^2)$ and $\mathcal{CN}(0, \sigma_{g}^2)$ and $\theta_1, \cdots, \theta_n$ $\in [\frac{-\pi}{Q}, \frac{\pi}{Q}]$ are i.i.d. uniformly distributed random variables with probability density function (PDF) given by $\frac{Q}{2\pi}$, where $Q$ is a integer number greater than 0. Additionally, we assume that $h_k$, $g_k$ and $\theta_k$ are independent for all values of $k$. 
UPDATED on 10/01/2020
Based on my simulations I know that the PDF of $r$, which is defined as
$$ r = \left|\sum_{k=1}^{n}{|h_k||g_k|\exp\left( j \theta_k \right)} \right|, $$ can be accurately approximated by a Gamma random variable even for small values of $n$, and $Q$, $e.g.$, $n=1$ and $Q=4$. Therefore, the PDF of $f = r^2$ can be approximated as $$ P_{n}(f) = \frac{1}{2 \Gamma(\kappa)\Theta^{\kappa}} f^{\left(\frac{\kappa-2}{2} \right)} \exp\left(-\frac{\sqrt{f}}{\Theta}\right), f > 0.$$
However, I'm not being able to find the $\kappa$ and $\Theta$, $i.e.$, shape and scale, parameters of the Gamma random variable.
This problem arises from the study on wireless communications channels and is of great importance to the research community.
 A: • n $\mathbf{\gg}$ 1: 
The variable $f=r^2$ is the square of the distance $r$ from the origin after $n$ steps of a random walk on the plane with random direction [$\theta$ uniformly in $(0,2\pi)$] and mean squared step size $s^2$ given by
$$s^2=\mathbb{E}(|h|^2||g|^2)=\mathbb{E}(|h|^2|)\mathbb{E}(|g|^2)=\sigma_h^2\sigma_g^2.$$
For $n\gg 1$ the distribution of $r$ is a Maxwell distribution and hence the distribution of $f$ is exponential,
$$P_n(r)=\frac{2r}{ns^2}e^{-\frac{r^2}{ns^2}}\Rightarrow P_n(f)=\frac{1}{ns^2}e^{-\frac{f}{ns^2}}.$$
So  for $n\gg 1$ the distribution $P_n(f)$ is a Gamma distribution with shape $k=1$ and scale $\theta=n\sigma_h^2\sigma_g^2$.

• n $=$ 1: 
For $n=1$ one can calculate directly
$$P_1(r)=\int_0^\infty d|g|\int_0^\infty d|h|\, \delta(r-|gh|)P(|h|)P(|g|)$$
$$=(\sigma_h\sigma_g)^{-2}\int_0^\infty dh\, \frac{4r}{h} e^{-h^2/\sigma_h^2} e^{-(r/h)^2/\sigma_g^2}=\frac{4r}{\sigma_g^2\sigma_h^2}K_0\left(\frac{2r}{\sigma_g\sigma_h}\right),$$
with $K_0$ a Bessel function.
(I have used that the absolute value $|h|$ when $h$ is $\mathcal{CN}(0, 1)$ has distribution $p(|h|)=2|h|e^{-|h|^2}$.)
So for $n=1$ the distribution of $f=r^2$ is
$$P_1(f)=\frac{2}{\sigma_g^2\sigma_h^2}K_0\left(\frac{2\sqrt{f}}{\sigma_g\sigma_h}\right).$$

• arbitrary n:
A closed form for $P_n(r)$ exists for arbitrary $n$, see The non-isotropic two-dimensional random walk, by B.C. Barber (1993):
$$P_n(f)=(\sigma_g\sigma_h)^{-n-1}\frac{2f^{(n-1)/2}}{\Gamma(n)}K_{n-1}\left(\frac{2\sqrt f}{\sigma_g\sigma_h}\right).$$
The large-$f$ decay is a stretched exponential $\propto e^{-2\sqrt{f}}$, so this cannot be precisely modeled by a Gamma distribution (which decays exponentially).

anisotropic random walk:
All of this was for an isotropic random walk, with $\theta_k$ at each step  uniformly distributed in $(-\pi,\pi)$. The parameter $Q>1$ restricts the angle to $(-\pi/Q,\pi/Q)$, so the scattering is peaked in the forward direction. For $n\gg 1$ the anisotropy can be accounted for by a rescaling of the mean squared step size,
$$s_Q^2=\frac{\sigma_h^2\sigma_g^2}{1-\langle \cos\theta\rangle},$$
where $\langle\cdots\rangle$ denotes the angular average. We have
$$\langle\cos\theta\rangle=\frac{Q}{2\pi}\int_{-\pi/Q}^{\pi/Q}\cos\theta\,d\theta=\frac{Q}{\pi}\sin(\pi/Q),$$
so we arrive at the distribution
$$P_n(f)=\frac{1}{ns_Q^2}e^{-\frac{f}{ns_Q^2}},\;\;s_Q^2=\frac{\sigma_h^2\sigma_g^2}{1-(Q/\pi)\sin(\pi/Q)}.$$
This is a Gamma distribution with shape $k=1$ and scale $\theta=ns_Q^2$.    
For $n=1$ the distribution $P_1(f)$ is $Q$-independent, so the result $\propto K_0$ above still applies. For $n>1$ I presume the decay will still be a stretched exponential, but I do not have an exact result as for $Q=1$.
A: I understand the answer I gave, which involves either exact or asymptotic results, is not quite what the OP was looking for, so let me try to suggest a different approach. The question is: please find parameters $k$ and $\theta$ such that 
$$p_{k,\theta}(r) = \frac{r^{k-1}}{ \Gamma(k)\theta^{k}}  e^{-r/\theta}$$
approximates the distribution of the distance from the origin after $n$ steps in a certain anisotropic random walk on the two-dimensional plane. The random walk has step size distribution 
$$F(r)=\frac{4r}{\sigma_g^2\sigma_h^2}K_0\left(\frac{2r}{\sigma_g\sigma_h}\right)$$
and scattering angle uniformly distributed in $(-\pi/Q,\pi/Q)$.
The first step would be to eliminate the parameters $\sigma_g$ and $\sigma_h$ by rescaling 
$$r=R\sigma_g\sigma_h$$ 
The distribution $P_{n,Q}(R)$ of $R$ no longer depends on $\sigma_g$ or $\sigma_h$. It only contains the two parameters $n$ and $Q$, both $\in\mathbb{N}$, and we wish to fit this with the two-parameter distribution $p_{k,\theta}(R)$, with $k,\theta\in\mathbb{R}$.
The next thing to note is that in the context in which the problem appears, only a few integer values of $Q$ are of relevance. We also have exact asymptotic results for $n\gg 1$, so if we would take 100 different combinations of $Q$ and $n$ we would have covered much of the relevant parameter space. The fitting procedure could be simplified by fitting first for large $R$, to extract $\theta$, and then for small $R$, to extract $k$. I guess this could be automated to a large extent, or perhaps even done by hand.

As an example, I tried this for $Q=1$, when I know the exact $P_{n,1}(r)$. The Mathematica command
Manipulate[
 LogPlot[{(4*r^n/Gamma[n])*BesselK[n-1, 2*r], 
(r^(k-1)/(Gamma[k]*t^k))*Exp[-r/t]},{r,0,20}],
 {n,2,10,1}, {k,0.1,5}, {t,0.1,5}]
provides for a convenient way to do the fit by sliding the parameters $n$, $k$, and $t\equiv\theta$. Here is a representative output for $n=3,5,8$:



I followed up on the suggestion by the OP to carry out the fit by adjusting the parameters $k$ and $\theta$ such that the first two moments of $R$ agree. For $Q=1$ we have the exact result
$$\mathbb{E}(R)=\frac{\sqrt{\pi } \,\Gamma \left(n+\frac{1}{2}\right)}{2 \Gamma (n)},\;\;\mathbb{E}(R^2)=n,$$
to be compared with the Gamma-distribution result $\mathbb{E}(R)=k\theta$,  $\mathbb{E}(R^2)=k(1+k)\theta^2$.
I show the comparison for $n=5$, when $k=2.955$ and $\theta=0.654$, see the plot below. It seems this fitting procedure is worse than the unconstrained fit (with $k=3.6$ and $\theta=0.55$) carried out above.

