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