Improving template with instantaneous amplitude

Given some recording (containing both the signal of interest and some Gaussian noise), can we improve our template (or model) of the signal using the analytic signal of the recording?

For example, a signal (red) with Gaussian noise results in some noisy recording (black). The instantaneous amplitude of the recording (black) deviates away from the true signal's instantaneous amplitude (red) because of this noise:

This got me thinking: Since we know the instantaneous phase of the true signal will be linear (when unwrapped), we can reconstruct this from the (noisier) instantaneous phase of the recording, say using OLS. For argument's sake, let's assume a fixed unit amplitude (even though in this example the amplitude is actually modulated at 2 Hz) and reconstruct the recordings either from the raw instantaneous phase (i.e., with noise) or from the reconstructed instantaneous phase (i.e., hopefully without noise) -- and compare these to the reconstructed signal from the instantaneous phase of the true signal (i.e., without noise).

Clearly, the reconstructed (OLS) instantaneous phase better approximates the true signal's instantaneous phase than the instantaneous phase we yield from the raw recording (error = true signal – reconstruction; figure below).

Finally, to see how close the reconstructions (either using raw or OLS instantaneous phase) approximate the true signal, I've computed these using the raw instantaneous amplitude for both. While (again) the OLS does better (below), there is still so much noise (because the instantaneous amplitude is computed from the raw recording).

Is there a way to better approximate the true signal (e.g., by low-pass filtering the instantaneous amplitude)?

• um, what's OLS? – Marcus Müller Dec 16 '17 at 11:16
• Thanks for your edit. OLS = ordinary least squares. Basically just estimating the instantaneous phase by fitting a line of best fit to the (noisier) raw instantaneous phase. – NBland Dec 16 '17 at 11:18
• so, do you have a latency or other constraint? In other words, do you need a "clean signal" estimate before you got all the recording, or is after ("offline") OK? – Marcus Müller Dec 16 '17 at 11:21
• because: making an assumption on the phase being linear and feeding that back is pretty much what every PLL does – and that works live, but the lower the filter bandwidth, the less noise you get at the output. Also, if you're really just after a single oscillation, you're just building a parametric spectrum estimator – here, superresolution technologies like ESPRIT might really work out great for you! – Marcus Müller Dec 16 '17 at 11:24
• It would all be offline. The trouble is all I know is the fundamental frequency of the signal -- I do not know the frequencies of the amplitude modulation. – NBland Dec 16 '17 at 11:26

Since I don't really know what you're using software-wise, I just whipped up a really short demo in GNU Radio Companion. This uses a PLL in a live processing demo:

If you're inclined to play around with it, here is the flow graph ready to be run if you've got a running installation of GNU Radio. The flow graph looks like this:

The left part up to the conversion to complex (which actually is a bit superfluous) is just the generation of a noise cosine.

The "PLL Freq Det" Block does pretty much what you describe: it compares the phase of the signal it observes with an internally generated oscillator, and generates a "phase signal" that is used to tell the internal frequency reference "how fast" it should oscillate. That signal is just a phase difference to be compensated, and needs to be internally filtered (the loop filter does that) before being applied to the internal oscillator.

At the same time, that signal is actually also a phase per sample - i.e., a frequency. So, you get a frequency estimate out of this.

Now, this is really practical for things like radio receivers that need to update their estimate of the received signal's frequency while operating, but in your case, you can do all this much more accurate, since you can actually work on the full signal.

A very popular (albeit, seldom optimal) estimator for the spectrum of a signal (and you're looking for which frequencies are making up your signal, so that's a question about the spectrum of your signal) is in fact the (magnitude of the) discrete fourier transform (DFT), often implemented by an FFT.

If you, for example, have your data in a Python program as vector signal, the following (untested) code should do interesting things to your signal:

import numpy
from matplotlib import pyplot

## Replace with your own samples and rate
samples = …
sampling_rate = 32e3

## Do a fourier transform of the whole data set
fourier_transform = numpy.fft.fft(samples)
positive_freq_components = fourier_transform[0:len(fourier_transform) / 2]
positive_freq_components_mag_squared = numpy.abs(positive_freq_components)**2
frequencies = numpy.linspace(0, sampling_rate / 2, len(fourier_transform) / 2)

## Plot that stuff!
pyplot.plot(frequencies, positive_freq_components_mag_squared)
pyplot.xlabel("Frequency (Hz)")
pyplot.ylabel("Digital Power")
pyplot.tight_layout()
pyplot.show()


This should give you a plot of power over frequency. Your noise will probably be relatively white, i.e. it will not have a lot of concentration in frequency domain. The sine, on the other hand, is very periodic and thus has all its energy at one frequency.

Note that I said "seldom optimal" about the DFT as spectrum estimator: It gives you the same frequency resolution all over your bandwidth, and you don't care about most frequencies. The (usable) frequency resolution depends linearly on the observation length, and you can't just get arbitrary frequencies by looking for peaks – all frequencies represented are sampling rate / length of DFT.

Other, "cooler", methods can give you a spectrum estimate that you can evaluate at arbitrary frequencies (such as MUSIC¹), or directly an estimate for the "strongest" frequency (or the strongest N frequencies) (such as ESPRIT¹). However, using these often requires a few strong assumptions, so I won't go deeply into that: If you want to know more about parametric estimators for fundamental frequencies, feel encouraged to research and ask a new question – there's folks on here way more experienced with such estimation tasks than I am.

• Thanks for this information. I am comfortable estimating the fundamental frequency of the signal -- and I have avoided DFT methods to do this for the same reasons as you say. I suppose my question is more on how we might separate noise in the recording from true amplitude modulations in the signal. Stabilising the instantaneous phase gets us a little bit closer, but all the progress is essentially lost because the instantaneous amplitude is computed from the (noisy) raw recording. – NBland Dec 17 '17 at 1:18
• The nice thing about the dft really is that it correlates the signal against all the oscillations with the possible frequencies. Noise doesn't correlate well - and thus, the longer your dft gets, the stronger the present frequencies are going to be in relation to the noise floor. – Marcus Müller Dec 17 '17 at 9:59
• And I suppose this would also be true of an LSSA approach. I have just avoided using the DFT and LSSA methods because I need frequency accuracy around the millihertz. – NBland Dec 17 '17 at 12:00
• I haven't tried this specific code, but what about: github.com/scivision/hires_spectrum/blob/master/matlab/esprit.m set n=1 for a single (the strongest) frequency. Try m in steps like 10, 20, 50, 64, 128, 1000, and watch out for the best results. (the return value is in radians per sample) – Marcus Müller Dec 17 '17 at 12:17