wavelength from peaks

Question

I have a simple 1sec signal of 1Hz added with 10Hz. I can detect all peaks (I use Python Scipy find_peaks) and I'm trying to deduce the 10Hz signal from those peaks. I can do it but it's not clean since the peaks spacing are not exactly the wavelength. Is there a formula that would adjust the wavelength according to the peaks height difference ? I'm trying to design a filter that removes the highest frequencies at each pass through the signal until there is none left.

and here's the Python code (without signal import or plotting):

_signal_build = np.zeros(len(_signal))
_sig_max, _ = find_peaks(_signal)
_sig_min, _ = find_peaks(-_signal)
for _current in range(0, len(_sig_max) -1):
      _wav_len = int(_sig_max[_current +1] - _sig_max[_current])
      _sine_start = _sig_max[_current] - round(_wav_len / 4)
      _sine_amp = 0.25*_signal[_sig_max[_current]] + 0.25*_signal[_sig_max[_current+1]] - 0.5*_signal[_sig_min[_current]]
      _sine = np.sin(np.linspace(np.pi, -np.pi + 2*np.pi/_wav_len, _wav_len)) * _sine_amp
      _signal_build[_sine_start:_sine_start + _wav_len] = _sine
      print(f"{_sig_max[_current]}-{_sig_max[_current+1]}, _wav_len: {_wav_len}, _sine_start: {_sine_start}, _sine_amp: {_sine_amp}")
_signal_filt = _signal - _signal_build

Thanks, cheers :)

Guillaume

Answer 1

In your comments you stated the purpose of this is to design a higher resolution transform than the FT offers. A consideration for the case of multiple frequencies that are not closely spaced is to use the adjacent bins in the DFT to interpolate a more precise frequency value. The simplest way to do this is to zero pad the time domain signal prior to computing the DFT, which will provide more samples of the Discrete-Time-Fourier-Transform (DTFT), revealing the frequency location of a localized single tone with much higher precision (the max value of the interpolated DFT will be the frequency location of a single tone).

Some precautions with this approach: Frequency resolution is theoretically the inverse of the time duration of the signal; for a non-windowed (or optionally stated using a rectangular window) each DFT bin in frequency represents the integrated energy under a Sinc shaped frequency response centered on that bin, with the first null of the main lobe of the Sinc function spaced one bin away. Thus we immediately see how this will not work to discern two closely spaced bins that may both be within the same main lobe (zero padding interpolates the frequency response but does NOT increase frequency resolution). For the case of one single frequency tone with zero padding the DFT result will be this Sinc function showing the additive correlation that would result to frequencies at any other locations under this curve, the peak of this will indeed be sufficiently close to the exact frequency location of the tone (within the precision of the number of samples we have provided with zero padding). However for two closely spaced frequencies, the interpolated result will appear as one Sinc function but will be the weighted and summed contribution of the two individual complex Sinc functions that each tone would give. Since the Sinc function shows relatively strong correlations even for tones spaced further away (the peak of the sidelobes of the Sinc only go down at rate $1/f$), it also shows how multiple tones can influence each other in the result depending on their relative locations. To combat this latter point, windowing should be used, which serves to increase the main lobe (increases the minimum distance two tones can be to discern each of them) while significantly decreasing sidelobes (so removes sensitivity to multiple tones further away).

Another precaution is that real signals appear as two tones in the DFT as demonstrated by Euler's relationship for a real cosine:

$$cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$$

Where here each exponential frequency ($e^{j\omega t}$) Is what actually maps to an individual tone in the frequency domain.

So for the case of higher or lower frequency real tones or low number of samples these two tones that accompany a real cosine or sine function can also be too close to influence each other to not provide an accurate result.

This is demonstrated by the simple case of the DFT of 5 samples for a frequency that is exactly at 1.5 bins for both a cosine and single exponential frequency where the blue stem plot is the result for the DFT and the interpolated line plot is the samples of the DTFT from zero-padding:

And repeated for the case of a real tone using a cosine of the same frequency:

Here we see the influence of the two closely spaced tones (in terms of number of bins) since one would be located at 1.5 bins and the other at 3.5 bins so there respective correlations influence the actual location of the peak (the frequency locations are both 2 bins apart as well as 3 bins apart due to the cyclical nature of the DFT). This example is the case of both a low real frequency in terms of number of bins and low number of samples. Even increasing the number of samples won't decrease the influence of the two tones as the two frequency locations will still be 3 bins apart regardless of number of time domain samples, so even windowing would not help in this case.

With that precaution hopefully better understood, here is a practical example of the approach with the case of 50 time domain samples and two real frequency tones at f1 = 12.3 bins and f2 = 20.5 bins, for the case of no windowing and windowing with a Kaiser window using $\beta=8$ (kaiser(50,8)). Zero padding was out to 10,000 samples, so providing for a resolution of the max value of 0.005 bins. With that the actual max values bin locations in the DTFT for both cases was:

No window: 12.255, 20.485

Kaiser window:12.300, 20.500

Also @CedronDawg has worked out the formulas for precisely computing the frequency of single tones which perhaps can be expanded to your case. See his blog post here: https://www.dsprelated.com/showarticle/1284.php

Answer 2

Dan gave you the reference to my article so you should be able to calculate the parameters of each tone more accurately.

Your iterative approach is how I do it as well. Find a tone, remove it, find the next, etc. However, you can save yourself a lot of calculations by removing the tone's DFT contribution directly from the DFT using the formulas (eq 25) here: https://www.dsprelated.com/showarticle/771.php. (or https://www.dsprelated.com/showarticle/1120.php if the frequency is very close to a bin/whole integer)

After the first pass you have a collection of tones. You can do converge to a better solution by doing further passes. For each tone remove the effects of the other tones from the DFT and recalculate the parameters. These calculations only need to be done on the two or three bins you use for the parameter estimation, not the whole DFT. This makes it a very efficient algorithm for tonal decomposition.