Whittaker models for $GL_n$ and Fourier coefficients Let $G$ be a compact abelian group. Then we know, because of the Peter-Weyl theorem, that $L^2(G)$ decomposes as a Hilbert space direct sum of 1 dimensional representations of $G$. 
Let $\mathbb{A}$ denote the adeles for $\mathbb{Q}$. Suppose we are given an automorphic funtion $\phi:GL_2(\mathbb{A})\to \mathbb{C}$ (I am following the book of Gelbart, Automorphic forms on adele groups, definition 3.3). In particular, $\phi(\gamma g)=\phi(g)$, for $\gamma\in G(\mathbb{Q})$. 
Fix $g\in G(\mathbb{A})$, then we get a function on $\mathbb{Q}\backslash\mathbb{A}$ defined as 
$$x\mapsto \phi\left(\left(\begin{array}{cc}1 & x\\
                          0 & 1\end{array}\right)g\right)$$
This being a smooth function on the compact abelian topological group $\mathbb{Q}\backslash \mathbb{A}$ we can write its Fourier series, and the Fourier coefficient for the character $\psi$ will be 
$$W_{\phi,\psi}(g)=\int_{\mathbb{Q}\backslash \mathbb{A}}\phi\left(\left(\begin{array}{cc}1 & x\\
                          0 & 1\end{array}\right)g\right)\psi(x)^{-1}dx$$
Using the Fourier coefficients we may write 
$$\phi(g)=\sum_{\psi}W_{\phi,\psi}(g).$$
For $GL_n$ where $n>2$ one defines analogous objects 
$$W_{\phi,\psi}(g)=\int_{N(\mathbb{Q})\backslash N(\mathbb{A})}\phi(xg) \psi(x)^{-1}dx$$
However, the space $N(\mathbb{Q})\backslash N(\mathbb{A})$ is no more a topological group, though it is compact. Instead of doing this, why don't we simply repeat the process for the $n=2$ case, that is, for $g\in GL_n(\mathbb{A})$ consider, for example, the inclusion of $\mathbb{A}$ in $GL_n(\mathbb{A})$ given by $x\mapsto U(x):=I+xE_{12}$ where $E_{12}$ is the matrix whose $12$ entry is 1 and all other entries are 0. The function $x\mapsto \phi(U(x)g)$ is a smooth function on $\mathbb{Q}\backslash \mathbb{A}$ and we can define the Fourier coefficients as before, that is,
$$W_{\phi,\psi}(g)=\int_{\mathbb{Q}\backslash \mathbb{A}}\phi(U(x)g)\psi(x)^{-1}dx$$
EDIT 
Then one has 
$$\phi(U(x)g)=\sum_{\psi}W_{\phi,\psi}(g)\psi(x)$$
Once again, put $x=0$ to get 
$$\phi(g)=\sum_{\psi}W_{\phi,\psi}(g).$$
I have not checked this very carefully, but, I guess, as mentioned in the response below we will also have 
$$ \phi(g)=\sum_{\gamma\in \mathbb{Q}} W_\phi\left(\begin{pmatrix}\gamma&0\\0&1\end{pmatrix}g\right).$$
My question is why do we need to work with the more complicated object $N(\mathbb{Q})\backslash N(\mathbb{A})$. Moreover, if we have to work with this object, then is there an analogous theory of Fourier series behind this. 
 A: It is indeed reasonable to wonder what's going on with "Fourier expansions" along non-abelian (sub-) groups... since, among other things, any one-dimensional representation has to factor through the maximal abelian quotient, so must lose some information.
For contrast: the unipotent radical $N$ of the standard minimal parabolic (upper-triangular matrices) in $GL_3$ is just barely not abelian. It is an example of a sort of Heisenberg group, and the representation theory of Heisenberg groups is well understood, via forms of the Stone-vonNeumann theorem (or, more generally, "Mackey imprimitivity"...) That is, there are one-dimensional repns which are trivial on the center $\pmatrix{1 & 0 & * \cr 0 & 1 & 0 \cr 0 & 0 & 1}$, and for each non-trivial character on the center a unique isomorphism class of irreducible. But this is not what we are doing in writing "Fourier expansions" of automorphic forms.
The way to obtain the Fourier expension on $GL_n$ is iterative, and does use the fact that Fourier expansions along abelian subgroups lose no information. Thus, the literal Fourier expansion along $N^{n-1,1}(\mathbb Q)\backslash N^{n-1,1}(\mathbb A)$ loses no information (where $P^{i,n-i}$ is the standard maximal proper parabolic with diagonal blocks of sizes $i$ and $n-i$ and $N^{i,n-2}$ is its unipotent radical). For cuspforms, the $0$th Fourier coefficient is $0$. Now comes the trick: the upper $n-1$ block of the Levi component of $P^{n-1,1}$ acts transitively on non-trivial characters on $N^{n-1,1}$, and the isotropy subgroup of a suitable "base choice" is close to being the $(n-2,1)$ parabolic in that block of the Levi component. So the $n-1,1$ Fourier expansion of a cuspform is a sum of translates of a single Fourier component (note, not just "coefficient").
Then there's an induction, which I won't re-type here... Leading to
$$
f(g) \;=\; \sum_{\gamma \in U(\mathbb Q)\backslash H(\mathbb Q)} W_{f,\psi}(\gamma g)
$$
where $H$ is the upper left $n-1$ block of the Levi component of $P^{n-1,1}$ and $U$ is the unipotent radical of the standard minimal parabolic in $H$, and
$$
\psi\pmatrix{1 & x_{12} & x_{13} & \ldots \cr 0 & 1 & x_{23} & \ldots \cr
 & & \ddots & \cr
 & & & 1 & x_{n-1,n} \cr & & & & 1
} \;\;=\;\; \psi(x_{12}+x_{23}+\ldots+x_{n-1,n})
$$
is a character on the unipotent radical of the minimal parabolic in $GL_n$, obviously vanishing on the derived group of the unipotent radical.
The induction/iterative procedure to obtain this expansion uses the assumption that $f$ is a cuspform repeatedly to know that no information is lost in successively dropping the zero-th Fourier component in literal Fourier expansions along various abelian unipotent radicals.
As GH noted, Jim Cogdell's various notes on $GL_n$ go through this process in some detail.
Edit: to comment on "why not expand along this abelian subgroup..." as in the edited version of the question: first, it is certainly legitimate to do so. The merit of the slightly labyrinthine version (after Shalika and P-S) is that (because of "uniqueness of (local) Whittaker models") the resulting Whittaker function _factors_over_primes_ (for cuspform $f$ generating an irreducible repn of the adele group). That is, up to a single global normalizing constant, $W_{f,\psi}$ is entirely determined by the local data attached to $f$. 
This factoring-over-primes feature will not hold for arbitrary or capricious decompositions-along-subgroups, although it is very important to appreciate the situations where it does.
So, again, the virtue of the Shalika-PS version of a "Fourier expansion" is that the Fourier components are essentially products of local (Whittaker) functions. 
A: One needs to work with the more complicated object $N(\mathbb{Q})\backslash N(\mathbb{A})$ in order to be able to recover the automorphic form. This is called the Whittaker model for $GL_n$, and it goes as follows. Fix a nontrivial character $\psi$ on $\mathbb{Q}\backslash\mathbb{A}$, and consider the Whittaker function
$$ W_{\phi}(g):=\int_{N(\mathbb{Q})\backslash N(\mathbb{A})}\phi(xg) \psi(x_{1,2}+x_{2,3}+\dots+x_{n-1,n})^{-1}dx. $$
The we have the expansion
$$ \phi(g)=\sum_{\gamma\in N_{n-1}(\mathbb{Q})\backslash G_{n-1}(\mathbb{Q})} W_\phi\left(\begin{pmatrix}\gamma&0\\0&1\end{pmatrix}g\right).$$
Cogdell in his lecture notes (in Lectures on automorphic $L$-functions, AMS, 2004) says on page 31: "As I said, the proof is not hard. The difficult thing, if there is one, is in recognizing that this is what one needs. This was recognized independently by Piatetski-Shapiro and Shalika."
EDIT. The OP's edit and Paul Garrett's comment show that my answer is not really sastisfactory. Perhaps the real advantage of using $N(\mathbb{Q})\backslash N(\mathbb{A})$ instead of $U(\mathbb{Q})\backslash U(\mathbb{A})$ lies in the fact that the resulting functions $W_\phi$ transform nicely with respect to all unipotent elements, namely
$$ W_\phi(xg)=\psi(x_{1,2}+x_{2,3}+\dots+x_{n-1,n})W_\phi(g)\quad\text{for}\quad x\in N(\mathbb{A}). $$
