tl;dr: What is a fast way to get the equation of this line?

I have to extract the "line of best fit" from this data (shown below), and many others like it (so I need a robust method). It is a a spectrogram of a linear chirp. The x axis is frequency and the y-axis is time. The linear chirp came from a truck outputting 4-140Hz linearly over exactly 8 seconds, but attenuation in the ground and some other physical phenomena mean that we are seeing a signal for around 10 seconds.

I looked into using the Hough transform to find the gradient of the line, and then cross-correlating the data with a line of this gradient. But the Hough transform is too slow to perform en mass, which is what I need to do.

I thought about finding the maxima for frequency ranges (5-10, 10-15, etc.) and plotting a line of best fit through these, but in this particular example there are only two or three areas of strong intensity.

The other spectrograms look like this but the high intensity regions will move around to different points on the line, and some may have more noise.

Whats a good way of doing this algorithmically?

My current idea is:

  1. Look for highest intesity regions.
  2. Look at many pixels in their neighbourhood and estimate the "center of mass".
  3. Draw line between them and find equation of line.

Is there a better way?

enter image description here

  • $\begingroup$ Are you generating the chirp, so you have knowledge of the chirp rate and/or time of arrival? If so, I would go back to the source data and try to dechirp it, such that the result is a near-vertical line in the spectrogram. That may be easier to pluck out and/or integrate further to suppress noise. $\endgroup$ – Jason R Dec 20 '12 at 13:28
  • $\begingroup$ The chirp is being generated by a seismic thumping truck ( geomore.com/wp-content/uploads/2012/04/… ). It should be giving a chirp of 4-120Hz over 8seconds but in the data we see (recorded from underground) it lasts 10 seconds. I cant fully explain why, am waiting for reply from someone more knowledge w.r.t. that. Sometimes it will last more than 10 seconds, sometimes less, so I do know that it will most certainly be 8-12 seconds. This does alter the gradient to observe by 1.5x, so I cant just generate the data I expect to see and 2d cross correlate it. $\endgroup$ – Mike Davies Dec 20 '12 at 14:03
  • 1
    $\begingroup$ I know nothing about seismic signal processing, but in general, it's not unusual for a signal to be extended in length by passing it through a system (in this case the earth). If the seismic channel is modeled as linear, then if that linear system has a lowpass frequency response, the observed output will be "smeared" across a longer time period. Not saying that's definitely what's happening here, but it could be an explanation. $\endgroup$ – Jason R Dec 20 '12 at 14:34

I think you can use standart method of fitting data to Straight Line. See for example explanation in 15.2 of "Numerical Recipes in C" http://apps.nrbook.com/c/index.html (page 661-662)

  1. To prepare arrays of $x_i$, $y_i$, $\sigma_i$, you must for each point on your image (frequency-time):
    $x_i$ := frequency
    $y_i$ := time
    $\sigma_i$ := 1/intensity (or 1/intensity$^2$)
    size of arrays will be N*M (N, M - x and y size of spectrogram in points)

  2. Calculate a and b (y = a + b*x) by using formula (15.2.6) from http://apps.nrbook.com/c/index.html

  3. If it will be some problem with outliers, use some robust fitting method. See for example 15.7 of http://apps.nrbook.com/c/index.html. But I think there will not such problem with your data.

This method is faster then Hough transform

It is very important to choose good expression for $\sigma_i$. May be better results will be with $\sigma_i$ := $1/\sqrt{intensity}$. May be will be better to use $\sigma_i$ := 1/(intensity-background) for (intensity > background) and exclude all points with (intensity < background). If there are some limitations on parameres of straight line (a,b), you can exclude some x-y points from $x_i$-$y_i$ arrays. But this modifications depend on actual data. It will be very interesting to see the your results.

  • $\begingroup$ SergV, is this a weighted linear regression? Are you considering all $N*M$ points in the image, but each with a weight of $\sigma_{xy}$? $\endgroup$ – Spacey Dec 21 '12 at 17:01
  • $\begingroup$ @Mohammad. Yes, it is correct. May be for best result you will need to modify expression for σ. For example, not σ=1/intensity, but σ = 1/(intensity-background) for intensity>background and exclude all points with intensity<background. But this modifications depend on actual data. It will be very interesting to see the results from Mike Davies. $\endgroup$ – SergV Dec 22 '12 at 5:47

Though the question is classic signal processing, you can use a computer vision technique called Hough line detector. The idea is that each pixel votes for all the possible lines that pass through it. Eventually, you select the line with the maximal amount of votes. In your case, you can weight each pixel vote by its intensity.

This step will give you a rough estimation of the line. Afterwards, you should throw away all outliers, find per each x coordinate center of mass (in Y) and do a least squares line fit on the points. I described something similar to the last part in this question.

Edit: After a second look on your data it seems that there are places where finding the center of mass might only confuse you. You should either remove this locations beforehand, or use some kind of RANSAC procedure.

P.S - if you post the original image, I can demonstrate the process in Matlab.

  • 2
    $\begingroup$ Great answer, thank you. I've opted for SergV's as it will run faster and I can implement it faster, but I got a lot of value from this. Thanks for taking the time. $\endgroup$ – Mike Davies Dec 21 '12 at 10:53

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.