P: n/a

Hi everyone,
I'm debugging an FFT routine and it has thrown up a perplexing error
which I cannot trace. Cannot anyone help? I've included the relevant
code below with comments:
rlft3(data_in1, speq_out1,1,512,256,1); // Perform FFT on 1st set of
input data
rlft3(data_in2, speq_out2,1,512,256,1); // Now 2nd set
// NB: data_in1/2, speq_out1/2 are indexed from 1
// data_in* are 3D arrays, speq_out* are 2D arrays
// data_in1[1][1][1] has "good data"
// Now do phase correlation
float *sp1, *sp2, r, im;
sp1 = &data_in1[1][1][1];
sp2 = &data_in2[1][1][1];
// data_in1[1][1][1] has "good data"
// Loop through elements multiplying and normalizing
for (int j=1;j<=(1*512*256)/2; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
// data_in1[1][1][1] has "good data"
// data_in1[1][1][1] has "good data"
// Do the same for the spectrum elements
sp1 = &speq_out1[1][1];
sp2 = &speq_out2[1][1];
for (j=1;j<=1024; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
int ab=1; // Does nothing  just put here for breakpoint purposes
// data_in1[1][1][1] is junk  "Expression cannot be evaluated" .
rlft3(data_in1, speq_out1, 1, 512, 256 , 1); // Do inverse FFT
Somewhere, the array data_in1 is being replaced by junk, although the
pointer value/memory address itself does not change. Whats happening?
TIA
Paul  
Share this Question
P: n/a
 pa**@paullee.com wrote:
Hi everyone,
I'm debugging an FFT routine and it has thrown up a perplexing error
which I cannot trace. Cannot anyone help? I've included the relevant
code below with comments:
rlft3(data_in1, speq_out1,1,512,256,1); // Perform FFT on 1st set of
input data
rlft3(data_in2, speq_out2,1,512,256,1); // Now 2nd set
You've NOT included the relevant code. Where is the declaration
of these functions. What is data_in1 declared as? What is Speq_out?
What are the definitions both in the caller and the claledd function?
You seem to pass in the size 512 256 1 but you don't use it at all in
the function.
What is fac?
>
// NB: data_in1/2, speq_out1/2 are indexed from 1
// data_in* are 3D arrays, speq_out* are 2D arrays
// data_in1[1][1][1] has "good data"
// Now do phase correlation
float *sp1, *sp2, r, im;
sp1 = &data_in1[1][1][1];
sp2 = &data_in2[1][1][1];
// data_in1[1][1][1] has "good data"
Why 1's here?
>
// Loop through elements multiplying and normalizing
for (int j=1;j<=(1*512*256)/2; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
Impossible to say, but my guess is you write off the end of the
allocation of the data_in1 and 2 arrays here.
>
// data_in1[1][1][1] has "good data"
// data_in1[1][1][1] has "good data"
// Do the same for the spectrum elements
sp1 = &speq_out1[1][1];
sp2 = &speq_out2[1][1];
for (j=1;j<=1024; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
int ab=1; // Does nothing  just put here for breakpoint purposes
// data_in1[1][1][1] is junk  "Expression cannot be evaluated" .
Possibly because there's a bug in this loop... again you possibly
write outside the allocation.
>
rlft3(data_in1, speq_out1, 1, 512, 256 , 1); // Do inverse FFT  
P: n/a

The code in those functions isn't relevant, its just there because I
anticipated someone saying "what are you trying to do here?"
what I'm interested in is that somewhere between:
for (j=1;j<=1024; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
and:
rlft3(data_in1, speq_out1, 1, 512, 256 , 1);
the data_in1 array gets filled with duff data. Only theres nothing
there to change it. I've put breakdpoints in (I'm using VC++ 6) and
examined the quanities in the arrays
and they seem fine until just before the final rlft3(...) call.  
P: n/a
 pa**@paullee.com wrote:
The code in those functions isn't relevant, its just there because I
anticipated someone saying "what are you trying to do here?"
what I'm interested in is that somewhere between:
for (j=1;j<=1024; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
and:
rlft3(data_in1, speq_out1, 1, 512, 256 , 1);
the data_in1 array gets filled with duff data. Only theres nothing
there to change it. I've put breakdpoints in (I'm using VC++ 6) and
examined the quanities in the arrays
and they seem fine until just before the final rlft3(...) call.
Possibly you're writing past the end of an array or corrupting the stack
somewhere (possibly before this code is even called). You're the only
one who can see enough of the code to make a judgment  you haven't
posted nearly enough for anyone to have the first clue. For example, we
can't see declarations for data_in1 etc., so we can't work out whether
you're writing off the end of any of the arrays.
One general comment is that array indexing from 1 is a bad idea in C++
(and in programming in general IMHO).
Tom  
P: n/a

I agree that starting arrays from 1 is bad, but unfortunately this is
something of a legacy project inherited from formula in a book....which
were translated from fortran, and its 1based arrays.
I'm still stumped.
The various arrays are defined as:
float *** data_in1=f3tensor(1,1,1,512,1,256); // width = 1...512,
height =1...256
float *** data_in2=f3tensor(1,1,1,512,1,256);
float **speq_out1 = matrix(1,1,1,1024);
float **speq_out2 = matrix(1,1,1,1024);
with the functions defined as:
#define NR_END 1
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl,
long ndh)
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh]
*/
{
long i,j,nrow=nrhnrl+1,ncol=nchncl+1,ndep=ndhndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
//if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t = nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)) );
//if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] = ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *)
malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(flo at)));
//if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] = ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i1]+ncol;
t[i][ncl]=t[i1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch]
*/
{
long i, nrow=nrhnrl+1,ncol=nchncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
//if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m = nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))) ;
//if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] = ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}  
P: n/a
 pa**@paullee.com wrote:
I agree that starting arrays from 1 is bad, but unfortunately this is
something of a legacy project inherited from formula in a book....which
were translated from fortran, and its 1based arrays.
I'm still stumped.
The various arrays are defined as:
float *** data_in1=f3tensor(1,1,1,512,1,256); // width = 1...512,
height =1...256
float *** data_in2=f3tensor(1,1,1,512,1,256);
float **speq_out1 = matrix(1,1,1,1024);
float **speq_out2 = matrix(1,1,1,1024);
with the functions defined as:
#define NR_END 1
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl,
long ndh)
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh]
*/
{
long i,j,nrow=nrhnrl+1,ncol=nchncl+1,ndep=ndhndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
//if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t = nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)) );
//if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] = ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *)
malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(flo at)));
//if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] = ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i1]+ncol;
t[i][ncl]=t[i1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
Yikes, all that just so you can use the syntax a[i][j][k]. But there's
an easier (or at least safer) way  C++ operator overloading. Obviously,
verifying that the above code is correct is not easy by inspection
(though it looks correct), because a simple out by one error anywhere is
likely to spell disaster. In general, such code should be avoided like
the plague.
Going back to your original code, one thing jumps out:
for (j=1;j<=1024; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
That appears to be writing off the end of your matrices (it's writing to
&speq_out1[1][1] + 2048 I think, while +1024 is actually the last
value), and presumably hosing the value of data_in1 in the process.
Tom  
P: n/a

>
Going back to your original code, one thing jumps out:
for (j=1;j<=1024; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
That appears to be writing off the end of your matrices (it's writing to
&speq_out1[1][1] + 2048 I think, while +1024 is actually the last
value), and presumably hosing the value of data_in1 in the process.
I've changed it to:
for (j=1;j<1024 + 1; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
 and it still doesn't work!  
P: n/a
 pa**@paullee.com wrote:
>Going back to your original code, one thing jumps out:
for (j=1;j<=1024; j++) { r = sp1[0]*sp2[0] + sp1[1]*sp2[1]; im = sp1[1]*sp2[0]  sp1[0]*sp2[1]; sp1[0] = fac*r/sqrt(r*r + im*im); sp2[1] = fac*im/sqrt(r*r + im*im); sp1+=2; sp2+=2; }
That appears to be writing off the end of your matrices (it's writing to &speq_out1[1][1] + 2048 I think, while +1024 is actually the last value), and presumably hosing the value of data_in1 in the process.
I've changed it to:
for (j=1;j<1024 + 1; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
sp2+=2;
}
 and it still doesn't work!
That's not changed anything! How is <1024+1 any different to <=1024?
It's still writing out to &speq_out1[1][1] + 2048. That's 1024 float
elements off the end of the allocated storage AFAICT, hosing 4K of heap.
Perhaps the loop should have <=512 as the condition? Only you can say,
but it's definitely wrong at the moment.
Tom  
P: n/a

I've checked the code and it seems to be OK. The memory reserved for
the four arrays don't overlap so theres no chance of any memory being
overwritten.
I'm bewildered. Its almost as if you go to a light swittch, turn it on
and off dozens of times, go to a sofa, read a book and the light turns
on by itself! Perhaps theres a ghost in the machine.... :)  
P: n/a
 pa**@paullee.com wrote:
I've checked the code and it seems to be OK. The memory reserved for
the four arrays don't overlap so theres no chance of any memory being
overwritten.
Of course the memory doesn't overlap, but if you write off the end of an
array, you might hit one of the other arrays (or any other heap memory
for that matter).
I'm bewildered. Its almost as if you go to a light swittch, turn it on
and off dozens of times, go to a sofa, read a book and the light turns
on by itself! Perhaps theres a ghost in the machine.... :)
I've already told you the problem. You have this code:
float **speq_out1 = matrix(1,1,1,1024);
float **speq_out2 = matrix(1,1,1,1024);
So &speq_out1[1][1] is a pointer to the first of 1024 floats (from
looking at your matrix function).
sp1 = &speq_out1[1][1];
sp2 = &speq_out2[1][1];
So sp1 is a pointer to the first of 1024 floats.
for (j=1;j<=1024; j++)
{
r = sp1[0]*sp2[0] + sp1[1]*sp2[1];
im = sp1[1]*sp2[0]  sp1[0]*sp2[1];
sp1[0] = fac*r/sqrt(r*r + im*im);
sp2[1] = fac*im/sqrt(r*r + im*im);
sp1+=2;
Here, sp1 is incremented by *2*. This happens 1024 times.
sp2+=2;
}
Hence on the final loop, you write to (&speq_out1[1][1])[2046]. BOOM!
Tom  
P: n/a

But it isn't speq1 thats causing problems  its data_in1 !  
P: n/a
 pa**@paullee.com wrote:
But it isn't speq1 thats causing problems  its data_in1 !
Writing off the end of speq1 may well involve overwriting data_in1!
They may well be placed sequentially in memory.
Incidentally, do you ever free the memory allocated by "matrix" etc.?
If not, you should, and it may show up some other heap corruption bugs
elsewhere in the code.
Tom  
P: n/a

Tom Widmer wrote:
>
Incidentally, do you ever free the memory allocated by "matrix" etc.?
If not, you should, and it may show up some other heap corruption bugs
elsewhere in the code.
The memory was allocated by:
> m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); m += NR_END; m = nrl;
Freeing that is going to be a whole lotta fun.   This discussion thread is closed Replies have been disabled for this discussion.   Question stats  viewed: 7474
 replies: 12
 date asked: Jul 7 '06
