text: Q:
  
  How do I find a 50Hz sin/cos signal in a measurement among other frequencies and noise?
  
  I measure magnetic fields with high precision, in a gradiometer setup with fluxgate magnetometers, sampling at 2000Hz. I can have big data sets if I want to - one is 600000+ samples long. I need to find the phase and amplitude of a 50Hz sine signal. However there are other (low frequency - 2.145Hz) stronger signals superimposed, along with noise. Because of those low frequency signals the DC offset is hard to eliminate, too. As a result, an FFT gives no results yet.
  Now it occured to me that it might be easier to find the 50Hz signal I look for (which is several orders of magnitude smaller then the DC offset and the other 2.1Hz signal) and determine its phase and amplitude, in stead of filtering out the rest. 
  What are good approaches and algorithms for that? I code this in python.
  edit: here is my python (with numpy) code for my carlos loop algorithm. 
     def decode_with_carlos_loop(self, frequency, periods_num):
          offsetfree_data = self.file_data
          shape = offsetfree_data.shape
          sampling_freq = 2000
  
          chunk_size = int(sampling_freq/frequency)
          alpha_step = .075
          indices = np.arange(chunk_size)
  
          for sensor_cnt in range(0, 4):
              for axis_cnt in range(1,3):
                  error_sin = error_cos = 0
                  cos_lo = sin_lo = []
                  theta_last = 0
                  theta = np.zeros(int(shape[1] / chunk_size) + 1)
  
                  for chunk_cnt, chunk_pointer in enumerate(range(0,shape[1],chunk_size)):
                      chunk = offsetfree_data[sensor_cnt, chunk_pointer:chunk_pointer+chunk_size, axis_cnt]
  
                      cos_lo = np.cos(2*np.pi*frequency*indices/sampling_freq + theta[chunk_cnt])
                      sin_lo = np.sin(2*np.pi*frequency*indices/sampling_freq + theta[chunk_cnt])
  
                      # mean is the lowpass filter
                      error_sin = np.mean(np.multiply(chunk, sin_lo))
                      error_cos = np.mean(np.multiply(chunk, cos_lo))
                      error = error_cos*error_sin
                      theta_last = theta[chunk_cnt+1] = theta[chunk_cnt] - alpha_step*error
                  print("sensor:axis:" + str(sensor_cnt) + ":" + str(axis_cnt) + "\tlast theta: " + str((theta_last + np.pi) % (2 * np.pi) - np.pi) + "\tcos_lo: " + str(cos_lo[0]) + "\tsin_lo: " + str(sin_lo[0]) )
  
  A:
  
  You can try tracking Phase-frequency of your $50 Hz$ signal using Costas Loop. Costas loop does not require the signal to be pre-processed in order to expose the desired frequency.
  I am not giving details of a Costas Loop because it can be found anywhere.It is pretty popular Carrier Recovery technique and a good starting point would be Wikipedia: CostasLoop 
  Instead, I would like to tell why I have chosen Costas Loop over other more common Squared-Diffrence Loop and Phase-Locked Loop :
  
  Squared-Difference Loop and PLL requires the signal to be pre-processed by a squaring non-linearity and Band-Pass filtering at $f_{center} = 2*f_c$. This is done to emphasize the desired frequency component. This preprocessing step is not required for Costas Loop.
  Phase Error sensitivity of Costas Loop is approximately double compared to PLL and Squared-Difference Loop. So, even smaller phase offsets in desired Frequencies are more accurately locked onto.
  There is only Low-Pass filtering required in Costas Loop which can be implemented pretty easily as a Moving Average implementation. As DC-Offset and other low frequency noise is very high which is leaving FFT technique not useful, you can use a low cutoff LPF of enough taps to get a sharp transition, in order to get accurate and almost noise free Error signal $cos(2(\phi - \theta))$.
  
  I have the MATLAB code implemented and customized for your situation at below path:
  CostasLoopPhaseRecovery
  The MATLAB Code runs n_runs number of Monty-carlo simulations to show that the algorithm will converge to the true phase of desired frequency eventually.
  There are few design parameters which will depend on scenario to scenario. Like, in your case since there is a Large DC offset and a very low frequency component of large magnitude, therefore, you will have to use a good Low-Pass Filter to filter out the Phase Difference signal $cos(2(\phi - \theta))$. I have used a Moving Average Filter for the purpose of Low-Pass Filtering, and so I had to increase the length of filter and increase step size $\alpha$ in order to get phase convergence accurate and faster.
  You will see a plot of Phase Convergence looking like below:
  
  The MATLAB code assumes a DC Offset muchlarger in magnitude than the desired signal component amplitude. I have added $\phi = 0.2$ in the desired signal and the Costas Phase recovery Loop converges to $\phi = 0.2$. There is an inherent ambiguity in phase recovery of $\hat{\phi} = \theta + n\pi$, which appears in the plot too, and it depends on the initial phase of the desired frequency which locally generated and multiplied to the incoming signal.
  Amplitude Estimation:
  Once you have a pretty good estimate of the phase $\phi$ of the sinusoid at frequency $50Hz$ you can generate a reference signal $x[n] = cos(2\pi. 50.nT_s + \hat{\phi}), \forall n \in \{ 0,1,2,3, \cdots , N-1\}$, where A is the parameter to be estimated. You can now use Least Squares technique to estimate the Amplitude as follows:
  $$\hat{A} = \frac{1}{N} <x,y>$$
  where, $<x,y>$ denotes inner product. The problem is that in doing so we have ignored the fact that our Noise is not uncorrelated (or White) but colored. So, this might result into very wrong estimates. The way to fix the estimation is by generating reference for whatever known frequencies are there in your signal and estimating Amplitude of them, and modeling other unknown frequencies as Colored Noise.
  So, in your case, you know that there is large DC offset in your signal and some small frequency component around $2.5Hz$. Assume DC, $2.5Hz$Sinusoid and $50Hz$Sinusoid amplitude as $A_o, A_{2.5} and A_{50}$. Let $y[n]$ be the measured signal, and then you can model $y[n]$ signal as : $$A_o.cos(2\pi 0.nT_s) + A_{2.5}cos(2\pi 2.5nT_s) + A_{50}cos(2\pi 50 nT_s + \hat{\phi}) + w(nT_s),$$where $w(nTs)$ is sampled colored noise (meaning Correlated).
  In Matrix Form it would be:
  $$\begin{pmatrix} y \end{pmatrix} = \begin{pmatrix} cos(2\pi 0.0T_s) & cos(2\pi 2.5.0T_s) & cos(2\pi 50.0T_s + \hat{\phi})\\ cos(2\pi 0.1T_s) & cos(2\pi 2.5.1T_s) & cos(2\pi 50.1T_s + \hat{\phi}) \\ \vdots&\vdots&\vdots\\cos(2\pi 0.N-1T_s) & cos(2\pi 2.5.N-1T_s) & cos(2\pi 50.N-1T_s + \hat{\phi})\end{pmatrix}. \begin{pmatrix}A_0\\A_{2.5}\\A_{50} \end{pmatrix} + \begin{pmatrix} w \end{pmatrix}$$ $$y = S.A+w$$
  Then the LS solution would be the following:
  $$\hat{A} = (S^HS)^{-1}S^H.y$$You can see that the noise is ignored again even though it is colored. The way to fix this is by estimating the Noise Covariance matrix and decolorising the noise or Whitening the noise, and then applying LS Technique. 
  You can also read about MAFI Algorithm to estimate Amplitudes of known Sinusoids in Colored Noise. That I hope will help you definitely. MAFI runs pretty close to Cramer-Rao Bound even at Low SNRs.
  
  A:
  
  Given that the FFT over that many samples does not show any results, your challenge may be in the overall spectral purity of the 50 Hz tone you seek. The FFT bin at 50 Hz is a correlation to that bin frequency which is the optimum detection in terms of SNR of a 50 Hz signal in the presence of white noise. The issue is the equivalent noise bandwidth of that bin is $1/T$ where $T$ is the total time duration of your signal. So if the energy of the signal you look for is wandering over several bins, then you will have reduced energy in each bin. Ultimately the spectral density of the signal you are looking for needs to be higher than the spectral density of the noise at whatever frequency you are searching in. If you don't see anything with the FFT it suggests that is not the case. 
  For non-stationary signals (such as the phase noise on a signal of interest at a specific tone) there is an optimum averaging time that will maximize the achievable SNR for the signal-- if your FFT is longer than this averaging time then your result will be degraded. With PLL approaches this same effect goes into choosing the tracking loop bandwidth. 
  
  