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]) )
  • 1
    $\begingroup$ Is the signal completely buried in the noise? Can you determine the signal to noise ratio or include a plot of what you're seeing? $\endgroup$
    – Engineer
    Apr 14 '20 at 13:49
  • $\begingroup$ Are any of the signal components synchronized to the sample clock? $\endgroup$ Apr 14 '20 at 16:29
  • 1
    $\begingroup$ “As a result, an FFT gives no results yet.” That sounds dubious, this problem is exactly what FT should excel at right of the box, nothing else needed. (Assuming it really is a perfect 50 Hz you're talking about, and it's a perfect linear superposition.) What window function have you been using? $\endgroup$ Apr 14 '20 at 22:19
  • $\begingroup$ FFT should work. I've used it to extract signal from very noisy sources. $\endgroup$ Apr 15 '20 at 0:28
  • $\begingroup$ Are you able to change/augment the measurement hardware? This sounds like exactly the problem a lock-in amplifier solves. $\endgroup$
    – Matt
    Apr 15 '20 at 2:25

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 :

  1. 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.

  2. 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.

  3. 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:


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.

  • $\begingroup$ He needs to acquire the phase of that known sinusoid too. That is why I have suggested a phase recovery Loop. It will be easier to converge to the true phase of desired frequency using a Phase recovery Loop. $\endgroup$
    – DSP Rookie
    Apr 14 '20 at 15:07
  • 1
    $\begingroup$ @DanBoschen Even I cannot see his comment now. But, yes we would need a second order PLL to track a ramping phase. $\endgroup$
    – DSP Rookie
    Apr 14 '20 at 17:20
  • 1
    $\begingroup$ I deleted it because I thought it was irrelevant to the question! (the OP just wanted to get the phase and amplitude which you have explained in your answer). But good to see the further discussion regarding 2nd order PLL. $\endgroup$
    – jithin
    Apr 15 '20 at 1:55
  • 1
    $\begingroup$ @DanBoschen The 2 PLLs are not in cascade actually, but the output $\theta_1[k]$ is added to the feedback $\theta_2[k]$ of the second PLL and the sum is fed to the frequency generator with phase = $\theta_1[k]+\theta_2[k]$. The same frequency generator is used for phase difference of both the PLLs. This is in context with digital PLL to track constant frequency offset and a phase offset in carrier frequency, which is equivalent to having a ramped phase offset. $\endgroup$
    – DSP Rookie
    Apr 15 '20 at 12:57
  • 1
    $\begingroup$ great, thanks. (it is 2.145Hz, not 2.5Hz, but whatever) $\endgroup$ Apr 17 '20 at 10:50

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.


One idea that comes to mind is to use a Hilbert filter. The lower frequency region needs to be tweaked in order to remove the offending 2.145Hz. Our power metering reference design example should help you get started, but you need to download the tool first.

I ran a quick test in the tool, and achieved the following:

enter image description here

  • $\begingroup$ are you trying to sell your product here? $\endgroup$ Apr 29 '20 at 8:01
  • $\begingroup$ Is that all you have to say? I would be a little more grateful that someone has taken timeout to answer your question. $\endgroup$
    – ASN
    Apr 29 '20 at 10:27
  • 1
    $\begingroup$ it was just a question to start with. i would need to download a tool, you dont specify how i would go about using it, and searching for "hilbert filter" left me guessing what it does. The website you link to does not explain it either. To me, it truely feels as if you wanted to place that link, most of all. If you actually tried to help, you didnt put a lot of effort into conveying your approach and much more into presentation and linking. Perhaps you could expand a little on your idea and the hilbert filter to balance that out. $\endgroup$ Apr 30 '20 at 11:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.