Accurately predicting the shape of the HI velocity function of galaxies is regarded widely as a fundamental test of any viable dark matter model. Straightforward analyses of cosmological $N$-body simulations imply that the $Lambda$CDM model predicts an overabundance of low circular velocity galaxies when compared to observed HI velocity functions. More nuanced analyses that account for the relationship between galaxies and their host haloes suggest that how we model the influence of baryonic processes has a significant impact on HI velocity function predictions. We explore this in detail by modelling HI emission lines of galaxies in the SHARK semi-analytic galaxy formation model, built on the SURFS suite of $Lambda$CDM $N$-body simulations. We create a simulated ALFALFA survey, in which we apply the survey selection function and account for effects such as beam confusion, and compare simulated and observed HI velocity width distributions, finding differences of $lesssim 50$%, orders of magnitude smaller than the discrepancies reported in the past. This is a direct consequence of our careful treatment of survey selection effects and, importantly, how we model the relationship between galaxy and halo circular velocity - the HI mass-maximum circular velocity relation of galaxies is characterised by a large scatter. These biases are complex enough that building a velocity function from the observed HI line widths cannot be done reliably.