Is it possible to do this procedure with different functions...?
Yes.
...for what reason do we not (at least not in any standard introductory texts I have ever read) do so?
No fundamental reason. It's just convenient in some contexts, that's all. In flat spacetime, using solutions with translation symmetry is often convenient, but it's not necessary. The theory as a whole still has translation symmetry even if use mode functions that don't individually have translation symmetry.
I'll explain how to formulate things using generic mode functions in flat spacetime. This approach is also important for quantum field theory in curved spacetime, because the curved-spacetime version of the Dirac equation typically doesn't have plane-wave solutions.
Quantization with general mode functions
We want the field operators to satisfy these equal-time anticommutation relations:
$$
\newcommand{\bfx}{\mathbf{x}}
\newcommand{\bfy}{\mathbf{y}}
\newcommand{\bfp}{\mathbf{p}}
\newcommand{\dpsi}{\psi^\dagger}
\big\{\psi_a(\bfx,t),\,\dpsi_b(\bfy,t)\big\}
\propto\delta_{a,b}\delta^3(\bfx-\bfy).
\tag{1}
$$
We also want them to satisfy the Dirac equation
$$
(i\gamma\partial+m)\psi(\bfx,t)=0.
\tag{2}
$$
The Dirac equation implicitly tells us how to express the field operators at times $t\neq 0$ in terms of the field operators at time $t=0$, so we need to check that the left-hand side of equation (1) really is independent of time. This consistency-check is a by-product of the following analysis.
Let $f_a(\bfx,t)$ and $g_a(\bfx,t)$ be ordinary complex-valued functions, where $a$ is a spinor index, and define
$$
(f,g)
\equiv
\sum_a\int d^3x\ f_a^*(\bfx,t)g_a(\bfx,t).
\tag{3}
$$
If $f$ and $g$ are solutions of the Dirac equation, then $(f,g)$ is independent of time. To prove this, use the fact that equation (2) can be written as $i\dot f=hf$, where the differential operator
$$
h=-\gamma^0\left(i\sum_{k>0}\gamma^k\partial_k+m\right)
\tag{4}
$$
satisfies $(g,hf)=(hg,f)$. This shows that $(f,g)$ is independent of time. Now let $M$ be a set of mode functions, solutions to the Dirac equation that are linearly independent and that satisfy the completeness relation
$$
\sum_{f\in M} f_a(\bfx,t) (f,g) = g_a(\bfx,t)
\tag{5}
$$
and the non-overcompleteness relation
$$
(f,g) = 0
\hspace{2cm}
\text{if }f,g\in M\text{ and }f\neq g.
\tag{6}
$$
The set $M$ of mode functions is not necessarily countable, so the "sum" in (5) should be understood in a suitably general sense. We can use these mode functions to solve equations (1)-(2). Let $\psi(f)$ be a time-independent operator "indexed" by the mode functions $f\in M$, and impose the anticommutation relations
$$
\big\{\psi(f),\,\dpsi(g)\big\}
\propto (f,g).
\tag{7}
$$
The mode operators $\psi(f)$ are time-independent by definition. This does not contradict the fact that $f$ is time-dependent, because $f$ is being used here only as an "index" to distinguish one operator from another. Now equations (1)-(2) are automatically satisfied by
$$
\psi_a(\bfx,t)=\sum_{f\in M} f_a(\bfx,t)\psi(f).
\tag{8}
$$
This generalizes the equation shown in the question. In hindsight, we can use equations (6) and (8) to derive the satisfying relationship
$$
\psi(f)=\sum_a\int d^3x\ f_a^*(\bfx,t)\psi(\bfx,t).
\tag{9}
$$
As a consistency check, notice that the right-hand side is actually independent of time for the same reason that $(f,g)$ is independent of time.
Constructing a vacuum representation
To construct a Hilbert-space representation with a lowest-energy (vacuum) state, we can impose one more condition on the mode functions: each one of them should involve either only positive frequencies or only negative frequencies. This is desirable because the positive- and negative-frequency parts of a field operator (or any other Heisenberg-picture operator) act as energy-lowering and -raising operators, respectively, as explained in my answer to another question:
Write $M=M_+\cup M_-$, where $M_+$ and $M_-$ are the sets of positive- and negative-frequency mode functions, respectively. Define a state $|0\rangle$ by the conditions
\begin{gather}
\psi(f)|0\rangle = 0
\hspace{1cm}
\text{for }f\in M_+
\\
\dpsi(f)|0\rangle = 0
\hspace{1cm}
\text{for }f\in M_-.
\tag{10}
\end{gather}
All other states obtained by acting on $|0\rangle$ with the other mode operators have higher energy than $|0\rangle$, so $|0\rangle$ is the vacuum state. We say that the other mode operators $\dpsi(f)$ with $f\in M_+$ and $\psi(f)$ with $f\in M_-$ create particle and antiparticles.
Recovering the version shown in the equation
To recover the version shown in the question, take the set of mode functions to be a complete set of single-frequency plane-wave solutions of the Dirac equation. To conform to tradition, write $b$ instead $\psi(f)$ when $f$ has positive frequency, and write $d^\dagger$ instead of $\psi(f)$ when $f$ has negative frequency. Equations (10) say that the operators $b$ and $d$ annihilate the vacuum state. Their adjoints $b^\dagger$ and $d^\dagger$ create particles and antiparticles.
Spontaneous particle creation in curved spacetime
The same idea of using generic mode functions also works in curved spacetime, where plane-wave solutions might not be available. In a time-dependent curved spacetime (or in the presence of other time-dependent background fields), the equation of motion may not preserve separate positive- and negative-frequency parts. (These can still be defined in asymptotic regions where spacetime is essentially flat.) In other words, a solution that has positive frequency initially may end up having a mix of positive and negative frequencies, and this can have interesting effects. We can choose the state to be annihilated by the mode-operators that have positive frequency initially, but then it won't be annihilated by the mode-operators that have positive frequency later. The result is that the initial vacuum state evolves into a non-vacuum state. This is how Hawking radiation is derived, using the equation of motion for the field in the time-dependent spacetime of a collapsing star.
