# How to get coordinates(x(time),y(frequency)) of spectral peaks from spectrogram?

I have created this spectrogram from a wav file. Please have a look:

As it can be observed clearly that my somewhat spectral peaks are visible. I want to get the x = time, y= frequency of these spectral peaks in order to find temporal displacement between the peaks. How can I get these coordinates?

I am looking for a programmatic approach (using python,java etc).

Regards,

Khubaib

• did you create the spectrogram from scratch? I'm asking this because if you did, it should be fairly obvious how to get the peak coordinates then. Commented Apr 14, 2020 at 10:25
• I created it by myself. But I want to use it in real time scenerio where two different smartphones play their specific frequency beeps just as above and I find the temporal displacement. Commented Apr 14, 2020 at 10:45
• Do you need the temporal displacement for the beeps on each row, or do you need the displacement between rows as well? Commented Apr 14, 2020 at 11:17
• @DanBoschen the temporal displacement between one beep from upper row and one beep for bottom row. Commented Apr 14, 2020 at 11:25
• ok if the beeps are of constant duration then what I outlined would be a good opproach since the FFT will give you the best estimate of the average pulse frequency and phases slope (delay) for each row. Commented Apr 14, 2020 at 11:29

Under the hood, your spectrogram consists of multiple, consecutive FFT frames (of length N) and each FFT frame contains N/2 frequency bins (even if you have N frequency bins, the upper part is discarded).

An FFT frame denotes a temporal position from the beginning of a signal and this temporal position can be computed as

time_pos = frame_number *  N / sampling_rate; //frame_number runs from 0......N-1


As an example, if your signal is sampled at 8000 Hz and N is 1024, the 10th frame will be at 10 * 1024 / 8000 = 1.28 sec from the beginning of the signal.

Just like every FFT frame has a fixed length (N samples long), each frequency bin has a fixed width in Hz (frequency resolution). This frequency resolution is given by the sampling rate divided by the frame/fft length so

freq_res = sampling_rate / N;


In order to find the frequency for a particular bin, simply multiply the bin index by the frequency resolution.

freq_10bin = 10 * freq_res;


In order to get the coordinates of the spectral peaks, you should iterate through all the frames and find those bins that have a value above some threshold thus indicating a peak. The frame indexes and bin indexes respectively are then used to compute the temporal and frequency coordinates for every peak (and the corresponding displacements).

• Thank you so much for a simpler and clear explanation. I have created this spectrogram using Python Scipy as : frequencies, times, spectrogram = signal.spectrogram(samples, sample_rate,nfft=256) Could you please explain a little bit more that how should I get the values from this available data? Commented Apr 14, 2020 at 12:17
• @Khubaib Ahmad, from what I can see in the scipy documentation, the spectrogram method returns an array of frequencies (bin frequencies) as well as an array of segment times (fft frames) so you should be able to get the coordinates from those 2 arrays. What I've written should be easy to understand. Commented Apr 14, 2020 at 12:45

Update for low latency real-time approach:

Given the two frequencies of interest $$f_1$$ and $$f_2$$, use a simple frequency rotated moving average filter such that the null of each filter is located at the other frequency location--- this is basically computing the two DFT bins directly using $$N =f_s/(f_2-f_1)$$. Implemented as a moving average filter, each output can trigger a threshold detector such that it would be easy to count the time between detections.

For example if $$f_s$$ is $$26$$ KHz and $$f_1$$ was at $$8$$ KHz, $$f_2$$ at $$10$$ KHz, then you would only need a moving average over 13 samples to compute this. For a non-integer case such as $$f_s$$ = $$25$$ KHz, then compute the moving average filter coefficients directly using $$c_n = e^{-j\omega_n}$$ where $$\omega_n$$ is the desired center frequency.

If the pulses are of fixed duration over the time block as shown, one approach would be to average all the y-axis values to get a single vector of magnitude vs frequency (averaged for all the time samples).

From this select the windowed maximum values over a frequency range using a threshold.

Then from the original data select the y row for each maximum value and take a zero-padded FFT of that row data. The zero-padded FFT will give you the best estimate of the average frequency over that row based on the lowest and strongest FFT bin. (Zero-padding will interpolate more samples of the DTFT which in an example like this will provide higher precision of the actual frequency.

The temporal displacement for the beeps on one row is the reciprocal of the frequency. For the temporal displacement for beeps on different rows you can derive this from the phases of the FFT results. For this you may need to decimate the row data first if the phase cannot be uniquely described with the data prior to decimation: Given a time delay between two rows in samples, the phase difference between the two will be a negative linear phase slope versus frequency extending from $$0$$ to $$-2\pi$$ for every one unit delay:

$$\phi[k] = -2\pi km/N$$

Where

$$k$$: Frequency index

$$m$$: delay in samples

$$N$$: total number of samples in row.

So you can use the above relationship to compare the phases of the frequency bins, but you must confirm using the above relationship that the frequency is low enough such that the phase does not cycle more than $$2\pi$$ by the time it has reached that frequency otherwise the result will not be resolvable. The solution then is to decimate the row first which will make the normalized frequency lower.

Alternatively, if the rows being compared have the same frequency and just a delay offset, then a cross-correlation in time would be a straight forward approach to finding the delay since it will have a maximum value at the delay offset between the two sets.

So in summary for $$x \times y$$ matrix where $$x$$ is time and $$y$$ is frequency:

1) Average all $$y$$ values versus $$x$$ to get $$v$$ which is a $$1 \times y$$ array.

2) Find the max values over a n sample window that are above a threshold $$k$$ where $$n$$ is the max width of what would be considered a unique frequency. (given as $$y=a,b,c...$$)

3) Use either the phase values of the FFT or a cross-correlation function to determined the delay between the pulses versus time on each of the relevant rows.

• I want to keep it simple. I don't want to involve more algorithmic delay in the computation as I want to find the temporal displacement in real-time while communication is still going on between the users. Commented Apr 14, 2020 at 12:19
• I see, for real time you can do 2 directly for each time t and then have a counter as t increments, and get the count for t2-t1 for when the next threshold is surpassed. Commented Apr 14, 2020 at 12:26
• @KhubaibAhmad For the simplest approach with least processing and delay would you consider an alternate approach with two simple FIR bandpass resonators rather than an FFT which is computing all the bins (This will work assuming the frequency of each unit is fixed) Commented Apr 14, 2020 at 12:59
• well thats a good approach. I'll let you know after trying it. Thanks Alot Dan Commented Apr 14, 2020 at 14:13
• Please take a look dsp.stackexchange.com/questions/66477/…. I'll update you soon once I am done with the algorithm. Commented Apr 15, 2020 at 12:43