I am working with a modular spectral processing STFT toolkit I have written in C. I often need to carry spectral peak detection and analisys (frequency estimation, peak displacement, etc) but I am into troubles when there are two or more peaks merging in a same apparent peak, i.e many overlapping peaks resulting in a single peak detectable as a single local maximum. I am aware that this problem has been addressed copiously in literature, in particular in chromatography or mass spectrometry. A commonly suggested approach is analysing the 2nd derivative looking for valleys, but this trick is not viable when the spectral resolution is insufficient, as when you have perhaps two or three peaks sharing few bins. Also I found many cases where a very low frequency sawtooth spectrum results in a smooth and large single peak, for which even the 2nd derivative of magnitude looks monotonous. So what to do ? I suspect that proper analisys of phases other than of magnitudes alone could give more clues of the actual underlieing peak structure but I am not a maths guru despite skilled enough. I would like to be able, once detected a whole peak "body", to decompose in its constituting overlapping peaks (if more than one) and possibly isolate them, or at least having a good estimate. Is there a way to do that, which is not too computationally intensive and which can be done in real time ? Thanks for any relevant answers!
-
$\begingroup$ Since you apply FFT in a part of your signal, are you able to get more samples of it? Increasing your "window size" will increase your spectral resolution and you might then be able to separate spectral peaks that are indistinguishable in a shorter window. $\endgroup$– GKHJan 12, 2020 at 20:52
-
$\begingroup$ Well that is an endless problem. Aside the fact that my toolkit is intended for real time music purposes so I prefer to use relatively small frame sizes to minimize latency, even by increasing the resolution you can always encounter overlapping peaks (as in cases of an audio signal with chorus and thus very close peaks) $\endgroup$– elenaJan 14, 2020 at 13:12
-
$\begingroup$ May I ask what is the target application then? Why do you need such a fine spectral analysis? It might help suggesting other solutions. $\endgroup$– GKHJan 14, 2020 at 22:13
-
$\begingroup$ Even just for accurate pitch shifting you often end up doing an approximate work because of carrying phase, frequency, magnitude estimation and peak displacement on apparent peaks which can actually contain more peaks overlapping, and as I said you cannot use frame sizes of 16384 or more for realtime applications because latency would become untolerable. $\endgroup$– elenaJan 15, 2020 at 14:03
-
$\begingroup$ But it is not even a matter of accurate analysis so much... as I said, I am encountering cases where even by using large framesizes, the spectrum of a harmonic signal of very low frequency like a sawtooth can easily appear as a * whole large peak*, if harmonic spacing is closer than fft resolution allows. $\endgroup$– elenaJan 15, 2020 at 14:03
5 Answers
A longer FFT on a longer window of actual data will better separate spectral peaks that share bins in a too short FFT.
so, adding to both Damien's a hotpaw's answer, i recommend using a Gaussian window and making your window as wide as practical. the Fourier Transform of a Gaussian function is a Gaussian function. so the wider your Gaussian window in the time domain, the narrower the Gaussian peaks will be in the frequency domain. you can use quadratic interpolation of the log of the peak magnitude to locate the true peak frequency (between FFT bins) more precisely.
the Gaussian function does not have sidelobe bumps. but, because a true Gaussian goes on forever, you actually have to truncate the Gaussian to a finite length, which does give it tiny sidelobes in the frequency domain.
if you want to estimate both frequency and frequency sweep rate, this paper of mine from 2001 discusses how to do that with a Gaussian window.
This answer and this answer and this answer have math regarding the Gaussian function and it's Fourier Transform.
Update: I found a trick which has the potential to work. It has still to be perfectioned though. I am explaining it here in the hope I can turn out useful or interesting for somebody else in future.
-Start from bin 0 of spectrum thru last bin;
-For every local maximum (peak tip) found, carry estimation of frequency (in fractionary bin units), phase and magnitude (such estimations can easily be made by parabola fitting or more advanced sinc-fitting methods or ok even by gaussian fitting but this is another topic);
-Basing on said estimates, resynthesize the peak by a suitable sinc function and subtract iteratively the whole synthesized peak (even just its right side) from the spectrum (for debugging purposes, the synthesized peaks can even be added to a separate spectral buffer);
-This way all hidden peaks will eventually come out.
This method can't by any means be defined mathematically sound and perfect of course, but is promising. A demonstration of that is that converting to audio all stft spectral frames obtained in the debug buffer said above, the result has minimal artifacts and is 99% close to the original (by ear). I also make the algorithm marking all detected peaks for debug purposes and I would say it does its job (it has still to be perfectioned though as I said but it basically works). It is even capable of separating peaks with very close frequency like a sinewave with chorus.
-
-
1$\begingroup$ I do everything in realtime. On modern multicore 64bits cpus there is no issue with that. I am always writing highly optimized C/Cpp code. Even the sinc generation can be optimized and made very fast. I am using typedef structs for complex numbers and inline complex operations exploiting vectors and achieving great speed. I am using tables for trigonometric stuff. There are still occasions where this method skips some hidden peaks though. But as I said it is a trick with the potential of being highly improved. $\endgroup$– elenaJan 17, 2020 at 13:01
Try a different windowing function - with windowing you can either have good spectral resolution or low spectral leakage.
https://en.wikipedia.org/wiki/Window_function
So if you look at frequency spectrum of each window, it has a mainlobe and many sidelobes. You can window with functions that have a large main lobe and small side lobes, thus reducing spectral leakage but at the cost of lower frequency resolution. Or if you desire greater frequency resolution you look for a window with a thinner main lobe at the cost of more spectral leakage.
The reason for this is when you truncate a signal for fourier analysis you are convoluting your signal with a window of that width - applying different windows will therefore affect your resultant output.
-
$\begingroup$ I am using hanning windowing by default (rised cosine) since I found out that is the most "natural" kind of window you can use (it has no sidelobes, which can make you detect false peaks where there are sidelobes instead; it gives the narrowest peaks of all possible window functions; it translates to a simple 3-points convolution in frequency domain). Thanks for the links, I can investigate the Gaussian window usage but that won't solve the issue of overlapping peaks imho. And as I said, I can't increase resolution too much because for real time applications latency would be unbearable.. $\endgroup$– elenaJan 14, 2020 at 13:20
-
$\begingroup$ Pedantic me just to point out: Hanning window does have sidelobes, every window function has sidelobes and some amount of spectral leakage. In your opinion the window wont make a difference but mathematically it does make a difference $\endgroup$ Jan 14, 2020 at 13:32
-
$\begingroup$ It is just a matter of terms Damien, as I said I am not a math guru, with sidelobes I mean the bumps typical of rectangular windows, the hanning window produces very narrow peaks without bumps which would challange any peak detection algorithm based on local maxima $\endgroup$– elenaJan 15, 2020 at 14:07
-
1$\begingroup$ Padding doesn't increase resolution. It only create a denser grid. Unless you need better accuracy in the value of the peak I don't see any sense in massive padding. $\endgroup$– RoyiJan 16, 2020 at 10:03
-
$\begingroup$ @Royi good point, hadn't said that correctly in my post $\endgroup$ Jan 16, 2020 at 13:32
Astronomers are talking about the Rayleigh distance. The (angular) distance between two point sources (stars) where they can be resolved with some robustness using diffraction limited optics: https://courses.lumenlearning.com/physics/chapter/27-6-limits-of-resolution-the-rayleigh-criterion/
This article by Julius Orion Smith discuses spectral interpolation of (well resolved) sinoidals: https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
If you want more frequency resolution and time resolution (as in: you cannot trust the input to be stationary for any longer), then you are in a tough situation. Perhaps fitting a smoothly sweeping sparse set of sinoids where you restrict amplitude or phase relationship?
