Code Snippet to Convert ADC Samples to I/Q I'd like to pull a narrowband signal out of ADC samples, where the ADC rate is much higher than the bandwidth/central frequency of the narrowband signal.  I'd like to feed the narrowband signal into a software package that want I/Q samples, so I'd like to know if a code snippet is doing the right thing.  There are lots of descriptions online, which seem to be ambiguous about factors of 2, and I've so far failed at finding a code example.
Let's say the ADC operates at sampling frequency $f_{samp}$, and the narrowband signal has central frequency and bandwidth $f_0$ and $bw$.  I think what I'm supposed to do is take my ADC samples, multiply by cos/sin($2\pi f_0 t$), then lowpass/downsample that to get I/Q. What exactly does that lowpass/downsampling look like, though?  I would have naively expected to set the output sample frequency to be the Nyquist rate for the desired bandwidth for both I and Q.  But... that just seems crazy?  I can't imagine people would carry around two Nyquist-sampled timestreams for one signal, so I have to be missing something.  Since nothing is as concrete as code, here's a snippet (ignoring edge effects, windowing, etc.) that does what I think I've read but which I firmly believe must be wrong.  Any suggestions on what it should actually look like would be most appreciated!
import numpy as np

def adc2iq(samples,f_samp,f0,bw):
    #Convert real ADC samples taken at sampling frequency f_samp to I/Q                                                                                                            
    #data with full bandwidth bw centered at f0

    n=len(samples)
    #make vectors of sin/cos at frequency f0                                                                                                                                         
    t=np.arange(n)/f_samp
    mycos=np.cos(t*f0*2*np.pi)
    mysin=np.sin(t*f0*2*np.pi)

    #take the Fourier transform of the samples times the sin/cos vectors                                                                                                             
    ift=np.fft.rfft(samples*mycos)
    qft=np.fft.rfft(samples*mysin)

    #low-pass/downsample by inverse FFTing a suitable number of samples                                                                                                              
    keep_frac=2*bw/f_samp  #is this right?                                                                                                                                           
    kmax=np.int(keep_frac*len(ift))
    i=np.fft.irfft(ift[:kmax])
    q=np.fft.irfft(qft[:kmax])
    return i,q

(In case anyone is curious, my motivation is trying to process GPS signals inside 550 MHz worth of L-band radio astronomy data, and the GPS software seems happiest with I/Q samples)
Consider a real carrier at frequency $f_o$ with phase modulation given by
$$s(t) = 2\cos(2\pi f_o t + \phi(t))$$
Where for simplicity of explanation the amplitude was normalized to be $2$.
This real carrier is at positive and negative frequencies as given by Euler's formula, both of which contain the phase modulation (conjugate symmetric):
$$s(t) = 2\cos(2\pi f_o t + \phi(t)) = e^{j( 2\pi f_o t + \phi(t))}+  e^{-j (2\pi f_o t + \phi(t))}$$
If we multiply the signal by $e^{-j2\pi f_o t}$, which is by Euler's formula $\cos(2\pi f_o t) - j\sin(2\pi f_o t)$, it will move the positive frequency to baseband (DC) and the negative frequency to a negative frequency at twice the rate of the carrier (which we must filter out):
$$e^{-j2\pi f_o t} s(t) = e^{-j2\pi f_o t}e^{j( 2\pi f_o t + \phi(t))}+  e^{-j2\pi f_o t}e^{-j (2\pi f_o t + \phi(t))}$$
$$= e^{j \phi(t)} + e^{-j(4\pi f_ot + \phi(t))}$$
Thus the low pass filter is very wide, passing bandwidth $bw$ at baseband and filtering out the signal up at $-2f_o$.
This is the primary task in recovering the signal at baseband (here in this example the complex signal $e^{j\phi(t)}$ with $I = cos(\phi(t))$ and $Q= sin(\phi(t))$. Any further down-sampling to lower sampling rates after the filtering is just best practice since it is desirable in terms of cost, size, and power to run at the lowest sampling rate possible (minimizes resources needed).
In general, if a modulated signal does not have a symmetric (complex conjugate symmetric) spectrum about its carrier then a complex signal must be used to represent that signal at baseband (using two real signals as I and Q), since any real signal has a complex conjugate symmetric spectrum so one real signal alone would not be able to represent the asymmetric spectrum. The GPS C/A code waveform as the OP is using is however symmetric but if there is any residual small frequency offset (due to Doppler and local clock offsets) then there will be a small positive frequency spectrum offset mostly overlapping with the small negative frequency offset with only one real signal at baseband. By having a complex baseband we can represent that spectrum offset in frequency and therefore the receiver can properly remove that- so the complex baseband processing simplified the receivers ability to resolve frequency and phase offsets between the received signal and the receiver’s local reference.
Please see this additional post for further examples of translating to other frequencies with complex signals:
Frequency shifting of a quadrature mixed signal