What is the distribution of the Wilcoxon signed-ranked test statistic? Wikipedia notes, "$W$ follows a specific distribution with no simple expression". Unencouraging as that sounds, its mean and variance are respectively $0$ and $\sum_{k=1}^{N_r}k^2$. I'm thus guessing $W\sim\sum_{k=1}^{N_r}kY_k$ for iids $Y_k$ of mean 0, variance 1. However, if I'm right the $Y_k$ aren't Gaussian, as then $W$ would be too.
Despite the title, my question is really what distribution these $Y_k$ have. I realise $W$ won't have a "neat" pdf, though its characteristic function would be nicer; if each $Y_k$ has cf $\phi(t)=\exp(-\frac{t^2}{2}+o(t^2))$, $W$'s cf will be $\prod_{k=1}^{N_r}\phi(kt)$.
 A: I'll give an outline of what's involved.
Let's look at the original (no ties) case. Note that a specific rank either contributes $R_k$ or $-R_k$ to the sum (with sign chosen randomly under the null). Let $S_k$ be the sign of the $k$th observation (/difference in a paired test) -- i.e. $-1$ if the difference is $<0$ and $+1$ otherwise, and without loss of generality assume the data are in rank order ($R_k=k$). 
Then $W=\sum_k k S_k$. With a fixed set of ranks, each term in the sum is independent but is not identically distributed. 
It's easy enough to write down the distribution for small $n$ by simply considering the possibilities.
e.g. with $n=3$ observations (pair differences in a paired test) the signs are
+ + +
+ + -
+ - +
+ - -
- + +
- + -
- - +
- - -

(where - indicates $S_i=-1$, etc. You can then simply multiply each sign by the rank (1 for the first column, 2 for the second, etc) and record the values of the statistics:
1+2+3=6, 1+2-3=0, 1-2+3=2, 1-2-3=-4, -1+2+3=4, -1+2-3=-2, -1-2+3=0, -1-2-3=-6

which is $1/8$ probability on $\pm 2,\pm 4,\pm 6$ and $1/4$ probability on $0$.
Automated complete enumeration is fairly practical for small $n$ -- say up to 20-25 or so (which is already to the point where we can use the asymptotic approximation fairly happily). We can extend this a little further by focusing only on the tail.
However, efficient combinatorial algorithms exist for computing these distributions in samples that are not huge - n=1000 takes less than a second on a modest laptop, and n=3000 -- which has 4.5 million values in it -- takes less than ten seconds (i.e. they work well out to far, far beyond where the asymptotics "kick in" effectively). 

The moment generating function and characteristic function
You could indeed write the mgf or the cf of these $kS_k$ terms and hence of the sum (I think the mgf of $S_k$ is $\cosh(t)$, the cf would be $\cos(t)$ but I'll let you check that, and you can get the mgf or cf of $kS_k$ from $M_{bX}(t)=M_X(bt)$, so I think the mgf of $kS_k$ should be $\cosh(kt)$, and similarly $\cos(kt)$ for the cf). The mgf of the sum is the product of the mgfs (and similarly for the cf), so we can at least get to that point fairly easily.
Asymptotics
Under the null, the mean of your statistic $W$ is $0$ (it's a sum of terms whose means are all $0$) and the variances of the components are $k^2 \text{Var}(S_k)=k^2$ (because the variance of an equally probable $\pm 1$ is $1$), so $\text{Var}(W)=\sum_k k^2$.
However, you don't need to use the characteristic function to show asymptotic normality, since you can just rely on the Lindeberg version of CLT; showing the variance condition holds is easy.
A: For I.I.D $Y_{k}$, the Gaussian assumption will not matter much in the limit. Even when random variables are independent with different distributions, they converge to a normal distribution in the limit under some regularity conditions. Please see  any standard large sample theory book for a reference (My illustration comes from section 2.7 in 'Elements of Large Sample Theory' by E.L. Lehmann. A picture of the relevant theorem is attached below)

