This is a bit awkward to explain without a bunch of diagrams, but here’s an attempt (and this is just one way of thinking about it, there’s a few different ways of explaining it that end up at the same result):
Let’s say you sample real values at 20MHz. You can interpret that in the frequency domain as if it had frequency components at 0..10MHz (the signal you want), and a mirror image at 0..-10MHz. i.e. this is what you’d get if you did a DFT on your input time-domain data to get frequency-domain values. The mirrored frequencies are because the input is real-valued only (that’s just how the DFT behaves with real-valued input, the negative frequencies are required to cancel out any imaginary values that would otherwise appear)
It’s also useful to remember here that while conventionally a DFT is going to produce values for frequencies in the range (-10..+10), you could equivalently treat it as if it was (+10..+30) or (-30..-10) or any other multiple-of-20MHz offset; numerically, plugging any of those aliased frequencies into a complex-exponential expansion would produce the same time-series data.
Now, you multiply your time-domain data by a complex-valued signal with frequency -Fs/4, which (and this is why you pick that magic frequency) is just this simple repeating sequence:
1 + 0j
0 - 1j
-1 + 0j
0 + 1j
Multiplying by -Fs/4 in the time domain is equivalent to convolving by -Fs/4 in the frequency domain, i.e. it will shift all the input frequencies by -Fs/4.
This multiplication is, I suspect, where the idea of “you can just interpret it as pairs of IQ samples” comes from – that’s not really what’s happening here though. The result of this multiplication is, conveniently, a stream of complex values where every second value is real-valued-only and the other values are imaginary-only, which lets you further optimize calculations later because there are lots of zeros that you can skip over.
When you shift that spectrum by -Fs/4 i.e. -5MHz by multiplying by a complex-valued oscillator at -Fs/4, you end up with this:
the signal you want at (0..+10) shifts to (-5..+5) (this is what you want to end up in the I/Q data)
the mirrored signal at (-5..0) shifts to (-10..-5); you don’t want this, it’s duplication of the signal
the mirrored signal at (-10..-5) has an alias at (+10..+15), which shifts to (+5..+10); you don’t want this either.
The LPF discards those mirrored signals at >5 and <-5MHz, since they’re not really useful when you’re working directly with complex IQ data, you only want the one copy of the signal and not the mirrored versions.
Finally, since you’ve just applied a LPF that threw away half the bandwidth, and there is now little-to-no energy <-5 or >+5MHz, you can also safely decimate-by-2 by just discarding every second sample. Or, equivalently, you can just calculate every second output of the LPF.
If the LPF happens to be a halfband FIR, then it turns out that the structure of the halfband FIR constants together with the alternating-real-and-imaginary values from the original -Fs/4 multiplication means that you can just apply a modified filter to the real-only values, and a simple delay to the imaginary-only values. This looks really weird if you just look at the final code, but it really is doing a full LPF-and-decimate-by-2, it’s just that there are a lot of zeros that can be ignored along the way.
If you didn’t apply the LPF and decimate-by-2, then the resulting baseband signal would be at twice the sampling rate you needed, and the extra bandwidth would be taken up with those mirror-images of the signal you want. Usually you don’t want that.
If you didn’t apply the LPF but did decimate-by-2 anyway, then reducing the sample rate would make aliases of the mirror-image copies overlap and interfere with the signal you want. This is approximately what “just interpret the real-valued samples as I/Q pairs” would do.
(I don’t know if that helped!)