Spreading sequence and equalization based on a paper
This is really a bit of a guess, since I haven't read the paper, but:

A spreading sequence can be designed/used to whiten a signal – that is the PSD of $y_n$ is the sum of a white process $w_n$ and a colored process $\sum\limits_{i=1}^3 h_i z_i$.

Based on observation of the receive signal you can estimate the PSD (e.g. by $\left|\mathrm{FFT}\right|^2$) and subtract your noise power density estimate (that you could get either from observing more bandwidth than your transmission, or by observing your band while you're not transmitting) for $w$; that way, you'd get the "spectral shape" of $h$, and thus, at least a start for equalization; for simple signals (like BPSK) equalizing the amplitude alone might suffice.

Now, for QAM etc., this might not really be the way to go. 

However, we know the following about the autocorrelation of both your $w$ and your $z$: 

Their autocorrelation is zero but for a single dirac impulse with the power of the signal at $\tau=0$. 
Note that autocorrelation is a property of a stochastic process that you, again, cannot observe directly, but estimate from a sample autocorrelation. That's why when you look at graphs that correlate a signal with itself, you don't see that perfect behaviour but non-zero: your observation is of finite length.

Now, with that knowledge, a discrete-time non-multipath reception of a signal should have exactly the same property: with increased duration of observation, the autocorrelation of the receive vector should converge to a dirac impulse.

With multiple paths, the "shortest" path's receive signal should, however, also correlate with the the other paths – and thus, you should get "peaks" for every $\tau$ that is a path length difference.

That's pretty nice!

The "boring" way of finding those peaks is actually computing the sample autocorrelation – but since correlation has time complexity $\mathcal O(N^2)$, and you might even want to interpolate your $N$ samples first to get higher resolution, you'd do it the fast correlation with the FFT. Notice that you can then comfortably increase the resolution by zero-padding.

The algorithm for that would be 
$$\begin{align}
 \mathbb R^K  \ni &amp;\quad \mathbf{\tau}_\text{peaks} \\
 &amp;=\text{largest}_K \big(R_{yy}(\tau)\big)\\
 &amp;=\text{largest}_K \big(\mathrm{IFFT}_N\big\{\mathrm{FFT}_N\left\{y\right\}(f)\cdot\mathrm{FFT}_N^*\left\{y\right\}(f)\big\}(\tau)\big)\\
&amp;=\text{largest}_K \big(\mathrm{IFFT}_{\text{hermitian},N}\left\{\left|\mathrm{FFT}_N\left\{y\right\}(f)\right|^2\right\}(\tau)\big)\\
\end{align}$$

The "$\text{hermitian}$" is but a hint – the FFT for real input data is faster to calculate in most cases than the general complex input case.

Notice how $\text{largest}_K \big(\mathrm{IFFT}\big)$ is a frequency estimation problem, where you're looking for $N$ frequencies – and looking for known-$N$ dominant frequencies especially well-solved by parametric estimators, where you can get very high resolution – without needing to pad the FFT in the first place. I'd like to point you to ESPRIT:
$$\begin{align}
 \mathbb R^N  \ni &amp;\quad \mathbf{\tau}_\text{peaks} \\
&amp;= \text{ESPRIT}_{N,K}\left(\left|\mathrm{FFT}_N\left\{y\right\}(f)\right|^2\right)
\end{align}$$