Search code examples
image-processingsignal-processingfftcorrelationcross-correlation

Correlation (offset detection) issues - Signal power concentrated at edge of domain


I'm in a bit of a bind - I am in too deep to quickly apply another technique, so here goes nothing...

I'm doing line tracking by correlating each row of a matrix with the row below and taking the max of the correlation to compute the offset. It works extremely well EXCEPT when the signals are up against the edge of the domain. It simply gives a 0. I suspect this is is because it is advantageous to simply add in place rather than shift in 0's to the edge. Here are some example signals that cause the issue. These signals aren't zero-mean, but they are when I correlate (I subtract the mean). I get the correct offset for the third image, but not for the first two.

Sample Bad signal 1

Sample Bad signal 2

Good Signal 1

Here is my correlation code

x0 -= mean(x0)
x1 -= mean(x1)
x0 /= max(x0)
x1 /= max(x1)


c = signal.correlate(x1, x0, mode='full')
m = interp_peak_offset(c)
foffset =(m - len(x0) + 1) * (f[2] - f[1])

I have tried clipping one of the signals by 20 samples on each side, correlating the gradient of the signal, and some other wonky methods with no success...

Any help is greatly appreciated! Thanks so much!


Solution

  • I ended up minimizing the average absolute difference between the two vectors. For each time shift, I computed the absolute difference/number of points of overlap. Here is my function that does so

    def offset_using_diff(x0, x1, f):
    #Finds the offset of x0 from x1 such that x0(f) ~ x1(f - foffset). Does so by 
    #minimizing the average absolute difference between the two signals, with one signal 
    #shifted.
    #In other words, we minimize |x0 - x1|/N where N is the number of points overlapping 
    #between x1 and the shifted version of x0
    
    #Args:
    #    x0,x1 (vector): data
    #    f (vector): frequency vector
    
    #Returns:
    #   foffset (float): frequency offset
    
    OMAX = min(len(x0) // 2, 100) # max offset in samples
    
    dvec = zeros((2 * OMAX,))
    offsetvec = arange(-OMAX + 1, OMAX + 1)
    
    y0 = x0.copy()
    y1 = x1.copy()
    
    y0 -= min(y0)
    y1 -= min(y1)
    
    y0 = pad(y0, (100, 100), 'constant', constant_values=(0, 0))
    y1 = pad(y1, (100, 100), 'constant', constant_values=(0, 0))
    
    for i, offset in enumerate(offsetvec):
        d0 = roll(y0, offset)
        d1 = y1
    
        iinds1 = d0 != 0
        iinds2 = d1 != 0
        iinds = logical_and(iinds1, iinds2)
        d0 = d0[iinds]
        d1 = d1[iinds]
    
        diff = d0 - d1
        dvec[i] = sum(abs(diff))/len(d0)
    
    m = interp_peak_offset(-1*dvec)
    foffset = (m - OMAX + 1)*(f[2]-f[1])
    return foffset