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.

10Hz not cleanly removed 10Hz fully removed

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 :)


  • 1
    $\begingroup$ Hi Hastoy and welcome to DSP.SE! It isn't entirely clear what you are asking-- could you add some plots and possibly formulas to show specifically the waveform and what you are trying to do? If you "mix" 1 Hz with 10 Hz that typically means multiply in time in which case the result would only be the sum and the difference (9 Hz and 11 Hz), except in practical implementation you can get some carrier feedthrough of 10 Hz but this would typically be 20 to 30 dB lower. And as far as your filtering, what exactly is the result you would like to achieve? None left of what? Plots will help. $\endgroup$ – Dan Boschen Mar 24 '20 at 2:59
  • $\begingroup$ thanks Dan for the comment, as you asked I added plots and a bit of Python code to make my question clearer. The idea of the filter is that I detect all peaks which gives me the highest frequencies present in the signal (one peak to the next). Then I generate sine waves and substract them from the signal. Then I do a second pass and null the "new" highest freqs, 3rd pass etc. until there is no more peak present in the signal so all sines will have been extracted. Trying to design a high resolution transform because Fourier is not accurate with short signals even in STFT. $\endgroup$ – hastoy Mar 24 '20 at 17:22
  • $\begingroup$ You can accomplish that from the DFT directly without too much trouble. @CedronDawg can help you- see his blog posts where he worked out the formulas for the exact frequency using the adjacent DFT bins. I can link his post later if her doesn’t see this first. $\endgroup$ – Dan Boschen Mar 25 '20 at 1:29

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:

Exponential Tone

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

Cosine Tone

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

Two tones, no window

Two tones, kaiser window

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

  • $\begingroup$ wow thanks for the "monster" answer Dan, much appreciated :) It'll take some time for me to digest it, I'll come back when I'm ready to post some results. Cheers, I wish you a bearable confinement ;) (C19 pandemic for guys reading this in the future) $\endgroup$ – hastoy Mar 25 '20 at 17:40
  • $\begingroup$ Bottom line: Window the time domain data, then zero pad the data then take the FFT, your max values will indicate the precise frequency locations. Follow the precautions listed regarding closely spaced tones. $\endgroup$ – Dan Boschen Mar 25 '20 at 17:43

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.

  • $\begingroup$ Hi Dan, thanks for the answer. I'm currently exploring wavelets but soon I will go back to Fourier strategies as they are so fast ! I will then read your blog and try your code, thanks again :) $\endgroup$ – hastoy Apr 25 '20 at 16:49

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.