Rune Allnor skrev: Found a couple of bugs. See comments in the code.> Actually, I found this to be an interesting algorithm, so I > jotted down the sketch of a C++ MEX program appended > below. > > Check the matlab documentation about how to compile > and link it. Do note that I have NOT tested or debugged it > and run-time error checking is absent, so you shouldn't > really regard it as more than an advanced pseudocode > to elaborate on what was said above. > > If you have a go at compiling and testing this code, > it might interest you to know that I often spend lots > of time debugging index limits around loops. That is, > I often get the non-trivial index limits in loops wrong. > > Have fun. I'd be very interested in seeing how things > turn ot if you get this algorithm to work, regardless > of whether you use my program or something else. > > Rune > > //////////////////////////////////////////////////////////////////////////////////////////////// > #include "mex.h" > #include "matrix.h" > > extern void __main(); > > void mexFunction(int nlhs, mxArray *plhs[], > int nrhs, const mxArray *prhs[]) > /* > * [rxy] = NlCorrMex(x,y) > * > * x Data series #1. N by 1 column vector. > * y Data series #2 M by 1 column vector. > * rxy Correlation sequence between x and y, rxy = [-1,1]. > * max(N,M) by 1 column vector. > * > * The function computes the "nonlinear covariance" > * between sequences x and y. The mean is removed from > * the shorter sequence, and the residual is normalized > * so the L2 norm equals 1. Frames of equal length > * as the shorter sequence are extracted from the longer > * sequence. The frames have their means subtracted > * and are normalized on a frame by frame basis. > * > * An effort is made to make a fast algorithm. > */ > > { > /* Check number of input arguments */ > if (nrhs~=2) > { > mexErrMsgTxt("Need exactly two input arguments!"); > } > > for (n=0;n<2;n++) > { > if((mxIsNumeric(prhs[n])==false)||(mxIsComplex(prhs[n])==true)) > { > mexErrMsgTxt("Arguments must be real-valued numbers."); > } > if (mxGetN(prhs[n])~=1) > { > mexErrMsgTxt("Arguments must be column vectors."); > } > } > > int N = mxGetM(prhs[0]); > int M = mxGetM(prhs[1]); > int n,m; > > double *xx = mxGetPr(prhs[0]); > double *yy = mxGetPr(prhs[1]); > double *tp; > > if (N>M) > { > n = M; > M = N; > N = n; > tp = xx; > xx=yy; > yy=tp; > } > > /* > * The pointer xx now points to the shorter data sequence > * and yy to the longer. N contains the number of samples > * in xx and M the number of samples in yy. > */ > > /* > * Allocate internal work space. > */ > > double *x = new double[N]; > > /* > * Allocate return variable > */ > > plhs[0]= mxCreateDoubleMatrix(M,1,mxREAL); > tp = mxGetPr(plhs[0]); > > /* > * Initialize work spaces. > * The shorter time series. > */ > > double mean = 0; > double variance = 0; > double t0; > > for (n=0;n<N;n++) > { > mean += xx[n]; > } > > for (n=0;n<N;n++) > { > variance += (xx[n]-mean)*(xx[n]-mean); > } > > t0=1/variance; > > for (n=0;n<N;n++) > { > x[n]=(xx[n]-mean)*t0; > } > /* > * The shorter data series is now prepared. > */ > > /* > * Initialize work spaces. > * The longer time series. > */ > for (n=0;n<N-1;n++) > { > mean += yy[n]; > } > > for (n=0;n<N-1;n++) > { > variance += (yy[n]-mean)*(yy[n]-mean); > } > > t0=1/variance; > /* > * The longer data series is now prepared. > */ > > /* > * The main loop. At entry, the variables 'mean' and 'variance' > * contain data for the sequence of length N-1 ending at y[n+N-2]. > * Update the 'mean' and 'variance' variables before computing > * the clocal covariance. Remove the contribution of y[n] before > * exiting the loop. > */ > > t1 = 1.0/((double)N); > > for (n=0;n<M-N+1;n++) > { > mean = mean + t1*yy[n+N-1]; > variance = variance + t1*(yy[n+N-1]-mean)*(yy[n+N-1]-mean);The variance must be computed as a loop over the frame in yy. The same mean must be subtracted from all the samples in the current frame.> t0 = 1/variance; > tp[n]=0; > for (m=0;m<N;m++) > { > tp[n] += t0*x[n]*(yy[n+m]-mean);/////////////////////////////////////////////////////////////////// Running index should be m: tp[n] += t0*x[m]*(yy[n+m]-mean); /////////////////////////////////////////////////////////////////> } > > variance = variance - t1*(yy[n]-mean)*(yy[n]-mean); > mean = mean - t1*yy[n]; > } > > /* > * Handle trailing N-1 samples in the correlation sequence. > */ > for (n=M-N+1;n<M;n++) > { > t0 = 1/variance; > tp[n]=0; > for (m=0;m<N-(M-N+1)-n;m++) > { > tp[n] += t0*x[n]*(yy[n+m]-mean); > } > > variance = variance - t1*(yy[n]-mean)*(yy[n]-mean); > mean = mean - t1*yy[n]; > } > > /* > * Release allocated memory > */ > delete [] x; > > return; > }
Cross Correlating Road Grade information
Started by ●April 20, 2006
Reply by ●April 21, 20062006-04-21
Reply by ●April 25, 20062006-04-25
>x : Sample (blue) data, N points long. >y : Template (red) data, M points long, M>>N. > >1) Subtract the mean from x once and for all. >2) Extract an N-point frame from y, call it y'. >3) Compute the mean of y', and subtract it from y'. >4) Compute the 0'th lag correlation coefficient Rxy' > between x and y', normalized to |Rxy'| = [0,1]. >5) Store Rxy' >6) Advance k points in the main sequence y >7) Repeat from 2) >this is exactly what my first 'working' algorithm was...however I found that I had to keep the step "k" very small. Which leads to a high computation time
Reply by ●April 25, 20062006-04-25
martini skrev:> >x : Sample (blue) data, N points long. > >y : Template (red) data, M points long, M>>N. > > > >1) Subtract the mean from x once and for all. > >2) Extract an N-point frame from y, call it y'. > >3) Compute the mean of y', and subtract it from y'. > >4) Compute the 0'th lag correlation coefficient Rxy' > > between x and y', normalized to |Rxy'| = [0,1]. > >5) Store Rxy' > >6) Advance k points in the main sequence y > >7) Repeat from 2) > > > > > this is exactly what my first 'working' algorithm was...however I found > that I had to keep the step "k" very small. > > Which leads to a high computation timeIt does, if you do things directly in matlab.With C or C++ in a MEX file, on the other hand... I posted a framework for hoow to do that, complete with bugs. Rune
Reply by ●April 25, 20062006-04-25
"Rune Allnor" <allnor@tele.ntnu.no> wrote in message news:1145994396.736127.49630@v46g2000cwv.googlegroups.com...> > martini skrev: >> >x : Sample (blue) data, N points long. >> >y : Template (red) data, M points long, M>>N. >> > >> >1) Subtract the mean from x once and for all. >> >2) Extract an N-point frame from y, call it y'. >> >3) Compute the mean of y', and subtract it from y'. >> >4) Compute the 0'th lag correlation coefficient Rxy' >> > between x and y', normalized to |Rxy'| = [0,1]. >> >5) Store Rxy' >> >6) Advance k points in the main sequence y >> >7) Repeat from 2) >> > >> >> >> this is exactly what my first 'working' algorithm was...however I found >> that I had to keep the step "k" very small. >> >> Which leads to a high computation time > > It does, if you do things directly in matlab.With C or > C++ in a MEX file, on the other hand... > > I posted a framework for hoow to do that, complete with bugs. > > Rune >I still don't understand why one wouldn't do this:>x : Sample (blue) data, N points long. >y : Template (red) data, M points long, M>>N. > >1) Subtract the mean from x once and for all - zero pad out to M+N-1 points >if need be. >3) Subtract the mean from y once and for all - zero pad out to M+N-1 points >if need be. >4) Compute the correlation of x and y in their entirety.Why bother with taking segments of y at all? Thus involving steps of some size, etc? It can be done in the frequency domain by zero padding x out to M+N-1 points and y out to M+N-1 points. The bigger the array the better as far as the savings in computation goes. Then there's no overlap of computations either. The results should be free from artifacts or algorithmic concerns caused by taking segments. What am I missing? Fred
Reply by ●April 26, 20062006-04-26
Fred Marshall skrev:> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message > news:1145994396.736127.49630@v46g2000cwv.googlegroups.com... > > > > martini skrev: > >> >x : Sample (blue) data, N points long. > >> >y : Template (red) data, M points long, M>>N. > >> > > >> >1) Subtract the mean from x once and for all. > >> >2) Extract an N-point frame from y, call it y'. > >> >3) Compute the mean of y', and subtract it from y'. > >> >4) Compute the 0'th lag correlation coefficient Rxy' > >> > between x and y', normalized to |Rxy'| = [0,1]. > >> >5) Store Rxy' > >> >6) Advance k points in the main sequence y > >> >7) Repeat from 2) > >> > > >> > >> > >> this is exactly what my first 'working' algorithm was...however I found > >> that I had to keep the step "k" very small. > >> > >> Which leads to a high computation time > > > > It does, if you do things directly in matlab.With C or > > C++ in a MEX file, on the other hand... > > > > I posted a framework for hoow to do that, complete with bugs. > > > > Rune > > > > I still don't understand why one wouldn't do this: > > >x : Sample (blue) data, N points long. > >y : Template (red) data, M points long, M>>N. > > > >1) Subtract the mean from x once and for all - zero pad out to M+N-1 points > >if need be. > >3) Subtract the mean from y once and for all - zero pad out to M+N-1 points > >if need be. > >4) Compute the correlation of x and y in their entirety. > > Why bother with taking segments of y at all? Thus involving steps of some > size, etc? > It can be done in the frequency domain by zero padding x out to M+N-1 points > and y out to M+N-1 points. The bigger the array the better as far as the > savings in computation goes. Then there's no overlap of computations > either. The results should be free from artifacts or algorithmic concerns > caused by taking segments. > > What am I missing?The OP posted a figure showing the type of data he is working with. The data were clearly NOT zero mean, nor stationary, which is what I am used to work with. The general problem is to detect a general signal shape that is shorter than the "typical variation" in the background mean. Take an example: x(n) = cos(pi*n/2)+10 y(n)= sin(pi*n/2)+ 3 The cos and sin terms (with zero mean) are orthogonal. Adding the mean terms destroy this orthiogonality, as is easily seen: One period, cos and sin terms only, zero mean: x1(n) = 1 0 -1 0 y1(n) = 0 1 0 -1 <x1,y1> = 0 Adding the means: x2(n) = 11 10 9 10 y2(n) = 3 4 3 2 <x2,y2> = 33+40+27+20 = 110 So two signal shapes that ought to be orthogonal end up with a non-zero correlation by adding a non-zero mean. The "spurious" correlation <x2,y2> above can easily mask the "true" match <y2,y2>, as can be seen <y2,y2> = 9 + 16 + 9 +4 = 38 These properties suggest to me that one ought to normalize the reference and the data on a frame-by-frame basis. I don't say this is impossible to implement as an FFT, but I can't see how it can be done. Again, my credo is the do things right before one tries do things fast. Rune
Reply by ●April 26, 20062006-04-26
"Rune Allnor" <allnor@tele.ntnu.no> wrote in message news:1146031107.226883.102990@g10g2000cwb.googlegroups.com...> > Fred Marshall skrev: >> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message >> news:1145994396.736127.49630@v46g2000cwv.googlegroups.com... >> > >> > martini skrev: >> >> >x : Sample (blue) data, N points long. >> >> >y : Template (red) data, M points long, M>>N. >> >> > >> >> >1) Subtract the mean from x once and for all. >> >> >2) Extract an N-point frame from y, call it y'. >> >> >3) Compute the mean of y', and subtract it from y'. >> >> >4) Compute the 0'th lag correlation coefficient Rxy' >> >> > between x and y', normalized to |Rxy'| = [0,1]. >> >> >5) Store Rxy' >> >> >6) Advance k points in the main sequence y >> >> >7) Repeat from 2) >> >> > >> >> >> >> >> >> this is exactly what my first 'working' algorithm was...however I >> >> found >> >> that I had to keep the step "k" very small. >> >> >> >> Which leads to a high computation time >> > >> > It does, if you do things directly in matlab.With C or >> > C++ in a MEX file, on the other hand... >> > >> > I posted a framework for hoow to do that, complete with bugs. >> > >> > Rune >> > >> >> I still don't understand why one wouldn't do this: >> >> >x : Sample (blue) data, N points long. >> >y : Template (red) data, M points long, M>>N. >> > >> >1) Subtract the mean from x once and for all - zero pad out to M+N-1 >> >points >> >if need be. >> >3) Subtract the mean from y once and for all - zero pad out to M+N-1 >> >points >> >if need be. >> >4) Compute the correlation of x and y in their entirety. >> >> Why bother with taking segments of y at all? Thus involving steps of >> some >> size, etc? >> It can be done in the frequency domain by zero padding x out to M+N-1 >> points >> and y out to M+N-1 points. The bigger the array the better as far as the >> savings in computation goes. Then there's no overlap of computations >> either. The results should be free from artifacts or algorithmic >> concerns >> caused by taking segments. >> >> What am I missing? > > The OP posted a figure showing the type of data he is working with. > The data were clearly NOT zero mean, nor stationary, which is what > I am used to work with. The general problem is to detect a general > signal shape that is shorter than the "typical variation" in the > background mean. > > Take an example: > > x(n) = cos(pi*n/2)+10 > y(n)= sin(pi*n/2)+ 3 > > The cos and sin terms (with zero mean) are orthogonal. > Adding the mean terms destroy this orthiogonality, as is > easily seen: > > One period, cos and sin terms only, zero mean: > > x1(n) = 1 0 -1 0 > y1(n) = 0 1 0 -1 > > <x1,y1> = 0 > > Adding the means: > > x2(n) = 11 10 9 10 > y2(n) = 3 4 3 2 > > <x2,y2> = 33+40+27+20 = 110 > > So two signal shapes that ought to be orthogonal > end up with a non-zero correlation by adding a > non-zero mean. > > The "spurious" correlation <x2,y2> above can > easily mask the "true" match <y2,y2>, as can > be seen > > <y2,y2> = 9 + 16 + 9 +4 = 38 > > These properties suggest to me that one ought to > normalize the reference and the data on a > frame-by-frame basis. I don't say this is impossible > to implement as an FFT, but I can't see how it can > be done. > > Again, my credo is the do things right before one > tries do things fast. > > RuneRune, It seems I failed to make the point I had in mind. I have no problem with removing the means - in fact I included that step. I only suggested not breaking things up into segments. I fail to see how that helps at all. Fred
Reply by ●April 26, 20062006-04-26
Fred Marshall skrev:> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message > news:1146031107.226883.102990@g10g2000cwb.googlegroups.com... > > > > Fred Marshall skrev: > >> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message > >> news:1145994396.736127.49630@v46g2000cwv.googlegroups.com... > >> > > >> > martini skrev: > >> >> >x : Sample (blue) data, N points long. > >> >> >y : Template (red) data, M points long, M>>N. > >> >> > > >> >> >1) Subtract the mean from x once and for all. > >> >> >2) Extract an N-point frame from y, call it y'. > >> >> >3) Compute the mean of y', and subtract it from y'. > >> >> >4) Compute the 0'th lag correlation coefficient Rxy' > >> >> > between x and y', normalized to |Rxy'| = [0,1]. > >> >> >5) Store Rxy' > >> >> >6) Advance k points in the main sequence y > >> >> >7) Repeat from 2) > >> >> > > >> >> > >> >> > >> >> this is exactly what my first 'working' algorithm was...however I > >> >> found > >> >> that I had to keep the step "k" very small. > >> >> > >> >> Which leads to a high computation time > >> > > >> > It does, if you do things directly in matlab.With C or > >> > C++ in a MEX file, on the other hand... > >> > > >> > I posted a framework for hoow to do that, complete with bugs. > >> > > >> > Rune > >> > > >> > >> I still don't understand why one wouldn't do this: > >> > >> >x : Sample (blue) data, N points long. > >> >y : Template (red) data, M points long, M>>N. > >> > > >> >1) Subtract the mean from x once and for all - zero pad out to M+N-1 > >> >points > >> >if need be. > >> >3) Subtract the mean from y once and for all - zero pad out to M+N-1 > >> >points > >> >if need be. > >> >4) Compute the correlation of x and y in their entirety. > >> > >> Why bother with taking segments of y at all? Thus involving steps of > >> some > >> size, etc? > >> It can be done in the frequency domain by zero padding x out to M+N-1 > >> points > >> and y out to M+N-1 points. The bigger the array the better as far as the > >> savings in computation goes. Then there's no overlap of computations > >> either. The results should be free from artifacts or algorithmic > >> concerns > >> caused by taking segments. > >> > >> What am I missing? > > > > The OP posted a figure showing the type of data he is working with. > > The data were clearly NOT zero mean, nor stationary, which is what > > I am used to work with. The general problem is to detect a general > > signal shape that is shorter than the "typical variation" in the > > background mean. > > > > Take an example: > > > > x(n) = cos(pi*n/2)+10 > > y(n)= sin(pi*n/2)+ 3 > > > > The cos and sin terms (with zero mean) are orthogonal. > > Adding the mean terms destroy this orthiogonality, as is > > easily seen: > > > > One period, cos and sin terms only, zero mean: > > > > x1(n) = 1 0 -1 0 > > y1(n) = 0 1 0 -1 > > > > <x1,y1> = 0 > > > > Adding the means: > > > > x2(n) = 11 10 9 10 > > y2(n) = 3 4 3 2 > > > > <x2,y2> = 33+40+27+20 = 110 > > > > So two signal shapes that ought to be orthogonal > > end up with a non-zero correlation by adding a > > non-zero mean. > > > > The "spurious" correlation <x2,y2> above can > > easily mask the "true" match <y2,y2>, as can > > be seen > > > > <y2,y2> = 9 + 16 + 9 +4 = 38 > > > > These properties suggest to me that one ought to > > normalize the reference and the data on a > > frame-by-frame basis. I don't say this is impossible > > to implement as an FFT, but I can't see how it can > > be done. > > > > Again, my credo is the do things right before one > > tries do things fast. > > > > Rune > > Rune, > > It seems I failed to make the point I had in mind. > > I have no problem with removing the means - in fact I included that step. > > I only suggested not breaking things up into segments. I fail to see how > that helps at all.Sorry, I didn't state it very clearly. I think removing the mean once and for all will work, provided the mean is stationary, or m(n,N) == m(n+k,N) where m(n,N) means "the mean of the N-sample frame starting at sample n." However, I don't think the premise of a the mean being stationary is valid for the data in question. I can't quantify it, maybe it would be interesting to plot the m(n,N) for some N, and see the variation and compare it to the max deviation from the mean inside the frames, max{|m(n,N)-x(k,n,N)|} where x(k,n,N) is the k'th sample from the N-length frame in the sequence x that starts at n. Rune
Reply by ●April 26, 20062006-04-26
> >Sorry, I didn't state it very clearly. I think removing the mean once >and >for all will work, provided the mean is stationary, or > >m(n,N) == m(n+k,N) > >where m(n,N) means "the mean of the N-sample frame starting at >sample n." > >However, I don't think the premise of a the mean being stationary >is valid for the data in question. I can't quantify it, maybe it would >be interesting to plot the m(n,N) for some N, and see the variation >and compare it to the max deviation from the mean inside the >frames, max{|m(n,N)-x(k,n,N)|} where x(k,n,N) is the k'th sample >from the N-length frame in the sequence x that starts at n. > >Rune > >for an infinetly long data set it should a stationary mean at zero, but as you have mentioned, for what I have shown (not a very long set) it will not be stationary at all. I will try out your mex file right now...
Reply by ●April 26, 20062006-04-26
martini skrev:> > > >Sorry, I didn't state it very clearly. I think removing the mean once > >and > >for all will work, provided the mean is stationary, or > > > >m(n,N) == m(n+k,N) > > > >where m(n,N) means "the mean of the N-sample frame starting at > >sample n." > > > >However, I don't think the premise of a the mean being stationary > >is valid for the data in question. I can't quantify it, maybe it would > >be interesting to plot the m(n,N) for some N, and see the variation > >and compare it to the max deviation from the mean inside the > >frames, max{|m(n,N)-x(k,n,N)|} where x(k,n,N) is the k'th sample > >from the N-length frame in the sequence x that starts at n. > > > >Rune > > > > > > > for an infinetly long data set it should a stationary mean at zero, but as > you have mentioned, for what I have shown (not a very long set) it will not > be stationary at all. > > > I will try out your mex file right now...Just make sure you read the follow-up where I commented on a couple of bugs I found. Rune
Reply by ●April 26, 20062006-04-26
what complier did you use to complie it? when i use the built in matlab one, or my visual C one, I am getting a lot of errors. just going in matlab and doing mex N1CorrMex.c (which is what i named it) right away on line 27 i get an error of: N1CorrMex.c(27) : error C2143: syntax error : missing ')' before '~' which doesnt make any sense because the line is if(nrhs ~= 2) plus I am getting all sorts of other errors... let me know, martini
