> just make sure that the integer lengths of the two input signals add to a > number no more than 1 sample longer than the size of the FFT buffer.what you mean? you mean to calculate the size of the output signal? btw, my function isnt working, i spotted another mistake, which is that my j loop was going up to j<blocksize instyead of j<blocksize/2, but still the output signal isn't what can be expected (what i use for test is a delta function convoluting another delta function, and i get like 8 non-zero samples in the output signal)
FFT convolution's complex multiplication problem
Started by ●August 12, 2005
Reply by ●August 14, 20052005-08-14
Reply by ●August 15, 20052005-08-15
Michel Rouzic wrote:> I think it should be alright like that : > > for (i=0; i<nblocks; i++) > for (j=0; j<blocksize/2; j++) > { > a=block[i][j]; > b=block[i][blocksize-j-1]; > c=smin[j]; > d=smin[blocksize-j-1]; > > block[i][j]=a*c - b*d; > block[i][blocksize-j-1]=b*c + a*d; > } > > but it's still giving me weird results, so there must be something > wrong still, do you see it?Without seeing the rest of the function it's difficult for me to say. However, you mentioned you are using a real-to-real FFT. As robert said, are you sure about the way data is being stored in the block[][] array? I recommend first coding it using the simple 1-dimensional complex FFT. Once you have that working, then go ahead and optimize for speed using the real-to-real FFT.
Reply by ●August 16, 20052005-08-16
BobM wrote:> Michel Rouzic wrote: > > I think it should be alright like that : > > > > for (i=0; i<nblocks; i++) > > for (j=0; j<blocksize/2; j++) > > { > > a=block[i][j]; > > b=block[i][blocksize-j-1]; > > c=smin[j]; > > d=smin[blocksize-j-1]; > > > > block[i][j]=a*c - b*d; > > block[i][blocksize-j-1]=b*c + a*d; > > } > > > > but it's still giving me weird results, so there must be something > > wrong still, do you see it? > > Without seeing the rest of the function it's difficult for me to say. > However, you mentioned you are using a real-to-real FFT. As robert > said, are you sure about the way data is being stored in the block[][] > array? > > I recommend first coding it using the simple 1-dimensional complex FFT. > Once you have that working, then go ahead and optimize for speed using > the real-to-real FFT.well i never even dealt with real-to-complex FFT, only with real-to-real, and well i thought that with the two corrections I did it would be ok (i dont see anything wrong in what i have corrected now). if you want i can post a version of the function that i reduced to the strict necessary, but that's still about 120 lines long..
Reply by ●August 16, 20052005-08-16
> > Michel Rouzic wrote: > > > I think it should be alright like that : > > > > > > for (i=0; i<nblocks; i++) > > > for (j=0; j<blocksize/2; j++) > > > { > > > a=block[i][j]; > > > b=block[i][blocksize-j-1]; > > > c=smin[j]; > > > d=smin[blocksize-j-1]; > > > > > > block[i][j]=a*c - b*d; > > > block[i][blocksize-j-1]=b*c + a*d; > > > } > > > > > > but it's still giving me weird results, so there must be something > > > wrong still, do you see it? > well i never even dealt with real-to-complex FFT, only with > real-to-real, and well i thought that with the two corrections I did it > would be ok (i dont see anything wrong in what i have corrected now). > if you want i can post a version of the function that i reduced to the > strict necessary, but that's still about 120 lines long..I do see something else that may be wrong here. It seems you have set up a symmetry between items in the array ({0,blocksize-1}, {1,blocksize-2}, etc.) in order to reference {real,imaginary} pairs. My assumption is that the item at index [0] should be DC and at [blocksize/2] should be Nyquist. This suggests that the symmetry should actually be like this: 0: DC {1,blocksize-1}: element 1 {real,imaginary} {2,blocksize-2}: element 2 {real,imaginary} . . . blocksize/2: Nyquist In other words, element [1] is paired with element [blocksize-1]. Right now you have element [0] paired with element [blocksize-1]. With this in mind I'd try substituting with this code: for (i=0; i<nblocks; i++) { a=block[i][0]; // DC c=smin[0]; block[i][0]=a*c; // no need for imaginary parts for (j=1; j<blocksize/2; j++) { a=block[i][j]; b=block[i][blocksize-j]; c=smin[j]; d=smin[blocksize-j]; block[i][j]=a*c - b*d; block[i][blocksize-j]=b*c + a*d; } a=block[i][blocksize/2]; // Nyquist c=smin[blocksize/2]; block[i][blocksize/2]=a*c; } But otherwise I recommend going with the complex FFT. There may be some dimensional issues going on in other parts of your code that you haven't shown. In my understanding you need to account for a size difference in the FFT array when using the real-to-real FFT. If you use the complex FFT, then you won't have to worry about that.
Reply by ●August 17, 20052005-08-17
BobM wrote:> > > Michel Rouzic wrote: > > > > I think it should be alright like that : > > > > > > > > for (i=0; i<nblocks; i++) > > > > for (j=0; j<blocksize/2; j++) > > > > { > > > > a=block[i][j]; > > > > b=block[i][blocksize-j-1]; > > > > c=smin[j]; > > > > d=smin[blocksize-j-1]; > > > > > > > > block[i][j]=a*c - b*d; > > > > block[i][blocksize-j-1]=b*c + a*d; > > > > } > > > > > > > > but it's still giving me weird results, so there must be something > > > > wrong still, do you see it? > > well i never even dealt with real-to-complex FFT, only with > > real-to-real, and well i thought that with the two corrections I did it > > would be ok (i dont see anything wrong in what i have corrected now). > > if you want i can post a version of the function that i reduced to the > > strict necessary, but that's still about 120 lines long.. > > I do see something else that may be wrong here. It seems you have set > up a symmetry between items in the array ({0,blocksize-1}, > {1,blocksize-2}, etc.) in order to reference {real,imaginary} pairs. My > assumption is that the item at index [0] should be DC and at > [blocksize/2] should be Nyquist. This suggests that the symmetry should > actually be like this: > > 0: DC > {1,blocksize-1}: element 1 {real,imaginary} > {2,blocksize-2}: element 2 {real,imaginary} > . > . > . > blocksize/2: Nyquist > > In other words, element [1] is paired with element [blocksize-1]. Right > now you have element [0] paired with element [blocksize-1]. > > With this in mind I'd try substituting with this code: > > for (i=0; i<nblocks; i++) > { > a=block[i][0]; // DC > c=smin[0]; > block[i][0]=a*c; // no need for imaginary parts > for (j=1; j<blocksize/2; j++) > { > a=block[i][j]; > b=block[i][blocksize-j]; > c=smin[j]; > d=smin[blocksize-j]; > > block[i][j]=a*c - b*d; > block[i][blocksize-j]=b*c + a*d; > } > a=block[i][blocksize/2]; // Nyquist > c=smin[blocksize/2]; > block[i][blocksize/2]=a*c; > } > > But otherwise I recommend going with the complex FFT. There may be some > dimensional issues going on in other parts of your code that you > haven't shown. In my understanding you need to account for a size > difference in the FFT array when using the real-to-real FFT. If you use > the complex FFT, then you won't have to worry about that.wow, thank you BobM, it works perfectly with the code you did. I think you pointed out something that I didn't understand about the structure of real-to-real FFT, it's what matches to the nyquist frequency. on a signal that is multiple of 2, if N/2 is the nyquist, what is the sample that symetrically matches to it (i guess i mean N/2 - 1)?
Reply by ●August 17, 20052005-08-17
Michel Rouzic wrote:> BobM wrote: > >>>>Michel Rouzic wrote: >>>> >>>>>I think it should be alright like that : >>>>> >>>>> for (i=0; i<nblocks; i++) >>>>> for (j=0; j<blocksize/2; j++) >>>>> { >>>>> a=block[i][j]; >>>>> b=block[i][blocksize-j-1]; >>>>> c=smin[j]; >>>>> d=smin[blocksize-j-1]; >>>>> >>>>> block[i][j]=a*c - b*d; >>>>> block[i][blocksize-j-1]=b*c + a*d; >>>>> } >>>>> >>>>>but it's still giving me weird results, so there must be something >>>>>wrong still, do you see it? >>> >>>well i never even dealt with real-to-complex FFT, only with >>>real-to-real, and well i thought that with the two corrections I did it >>>would be ok (i dont see anything wrong in what i have corrected now). >>>if you want i can post a version of the function that i reduced to the >>>strict necessary, but that's still about 120 lines long.. >> >>I do see something else that may be wrong here. It seems you have set >>up a symmetry between items in the array ({0,blocksize-1}, >>{1,blocksize-2}, etc.) in order to reference {real,imaginary} pairs. My >>assumption is that the item at index [0] should be DC and at >>[blocksize/2] should be Nyquist. This suggests that the symmetry should >>actually be like this: >> >>0: DC >>{1,blocksize-1}: element 1 {real,imaginary} >>{2,blocksize-2}: element 2 {real,imaginary} >>. >>. >>. >>blocksize/2: Nyquist >> >>In other words, element [1] is paired with element [blocksize-1]. Right >>now you have element [0] paired with element [blocksize-1]. >> >>With this in mind I'd try substituting with this code: >> >>for (i=0; i<nblocks; i++) >>{ >> a=block[i][0]; // DC >> c=smin[0]; >> block[i][0]=a*c; // no need for imaginary parts >> for (j=1; j<blocksize/2; j++) >> { >> a=block[i][j]; >> b=block[i][blocksize-j]; >> c=smin[j]; >> d=smin[blocksize-j]; >> >> block[i][j]=a*c - b*d; >> block[i][blocksize-j]=b*c + a*d; >> } >> a=block[i][blocksize/2]; // Nyquist >> c=smin[blocksize/2]; >> block[i][blocksize/2]=a*c; >>} >> >>But otherwise I recommend going with the complex FFT. There may be some >>dimensional issues going on in other parts of your code that you >>haven't shown. In my understanding you need to account for a size >>difference in the FFT array when using the real-to-real FFT. If you use >>the complex FFT, then you won't have to worry about that. > > > wow, thank you BobM, it works perfectly with the code you did. I think > you pointed out something that I didn't understand about the structure > of real-to-real FFT, it's what matches to the nyquist frequency. on a > signal that is multiple of 2, if N/2 is the nyquist, what is the sample > that symetrically matches to it (i guess i mean N/2 - 1)?DC and Fs/2 are cosine-only terms. The sine component of FS/2 disappears because it's sampled at its zero crossings. At DC sin(2*pi*0*t) is zero. Sometimes, to save space, the terms for DC and Fs/2 are stored together as if they were an imaginary pair. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●August 17, 20052005-08-17
Michel Rouzic wrote:> wow, thank you BobM, it works perfectly with the code you did. I think > you pointed out something that I didn't understand about the structure > of real-to-real FFT, it's what matches to the nyquist frequency. on a > signal that is multiple of 2, if N/2 is the nyquist, what is the sample > that symetrically matches to it (i guess i mean N/2 - 1)?Good to hear it worked! Jerry answered your question much more clearly than I could have. But how about a rundown of how an FFT result is typically stored.. Regarding the complex FFT, think of an 8-point FFT result array called "array[8][2]" which will hold the results. Here is the way this array might be arranged: REAL PARTS array[0][0]: DC array[1][0]: bin 1 array[2][0]: bin 2 array[3][0]: bin 3 array[4][0]: Nyquist array[5][0]: bin 5 (-3) array[6][0]: bin 6 (-2) array[7][0]: bin 7 (-1) IMAGINARY PARTS array[0][1]: DC array[1][1]: bin 1 array[2][1]: bin 2 array[3][1]: bin 3 array[4][1]: Nyquist array[5][1]: bin 5 (-3) array[6][1]: bin 6 (-2) array[7][1]: bin 7 (-1) Note that real parts are stored in [0] and imaginary are in [1]. Bins above Nyquist are the so-called negative frequencies. Note the symmetry around Nyquist to the respective positive frequencies. If you were to graph these out, you might considering plotting the bins above Nyquist below 0 (-3,-2,-1 shown above). For a real signal the real bins symmetrical to each other around Nyquist will be equal to each other. The imaginary bins symmetrical around Nyquist will be equal to -1 * each other (complex conjugate). DC and Nyquist will always be real for a real signal. The real-to-real FFT you're using takes advantage of these redundancies and is stored like this: array[0]: DC (real) array[1]: bin 1 (real) array[2]: bin 2 (real) array[3]: bin 3 (real) array[4]: Nyquist (real) array[5]: bin 3 (imaginary) array[6]: bin 2 (imaginary) array[7]: bin 1 (imaginary) Hope this helps, Bob
Reply by ●August 17, 20052005-08-17
yeah i understand better now. now i think i'll have to review some of my functions that do things in the frequency domain because I didn't understand perfectly the storage of the real-to-real FFT. Now, if a 8-point FFT result is like that :> array[0]: DC (real) > array[1]: bin 1 (real) > array[2]: bin 2 (real) > array[3]: bin 3 (real) > array[4]: Nyquist (real) > array[5]: bin 3 (imaginary) > array[6]: bin 2 (imaginary) > array[7]: bin 1 (imaginary)how would be a 7-point or 9-point FFT, I mean, what would disappear or appear? BobM wrote:> Michel Rouzic wrote: > > wow, thank you BobM, it works perfectly with the code you did. I think > > you pointed out something that I didn't understand about the structure > > of real-to-real FFT, it's what matches to the nyquist frequency. on a > > signal that is multiple of 2, if N/2 is the nyquist, what is the sample > > that symetrically matches to it (i guess i mean N/2 - 1)? > > Good to hear it worked! Jerry answered your question much more clearly > than I could have. > > But how about a rundown of how an FFT result is typically stored.. > > Regarding the complex FFT, think of an 8-point FFT result array called > "array[8][2]" which will hold the results. Here is the way this array > might be arranged: > > REAL PARTS > array[0][0]: DC > array[1][0]: bin 1 > array[2][0]: bin 2 > array[3][0]: bin 3 > array[4][0]: Nyquist > array[5][0]: bin 5 (-3) > array[6][0]: bin 6 (-2) > array[7][0]: bin 7 (-1) > > IMAGINARY PARTS > array[0][1]: DC > array[1][1]: bin 1 > array[2][1]: bin 2 > array[3][1]: bin 3 > array[4][1]: Nyquist > array[5][1]: bin 5 (-3) > array[6][1]: bin 6 (-2) > array[7][1]: bin 7 (-1) > > Note that real parts are stored in [0] and imaginary are in [1]. Bins > above Nyquist are the so-called negative frequencies. Note the symmetry > around Nyquist to the respective positive frequencies. If you were to > graph these out, you might considering plotting the bins above Nyquist > below 0 (-3,-2,-1 shown above). > > For a real signal the real bins symmetrical to each other around > Nyquist will be equal to each other. The imaginary bins symmetrical > around Nyquist will be equal to -1 * each other (complex conjugate). DC > and Nyquist will always be real for a real signal. > > The real-to-real FFT you're using takes advantage of these redundancies > and is stored like this: > > array[0]: DC (real) > array[1]: bin 1 (real) > array[2]: bin 2 (real) > array[3]: bin 3 (real) > array[4]: Nyquist (real) > array[5]: bin 3 (imaginary) > array[6]: bin 2 (imaginary) > array[7]: bin 1 (imaginary) > > Hope this helps, > Bob
Reply by ●August 19, 20052005-08-19
Michel Rouzic wrote:> how would be a 7-point or 9-point FFT, I mean, what would disappear or > appear?There is no Nyquist element for an odd size. (This is, by the way, documented in the FFTW manual.)
Reply by ●August 26, 20052005-08-26
BobM wrote: I tested my function further, and found out that it still doesnt work properly. When I convolve with a function that is an array of 20,000 samples set to 0, and the element 4,000 being set to 1, my output signal looks like the orginal input (with a longer length of course), unmodified, except that on every section of 32,768 samples, the first 4,000 are set to zero. I also tried with a windowed-sinc function, i used 88.2 Hz as a cutoff frequency, 8.82 Hz as the bandwidth of the rolloff, and when convolving with the input signal, not only the output signal has regions regularly set to zero, but the output signal doesnt look/sound at all like it's been low-passed filtered I made sure that my sinc function generator works properly and that my blackman function works properly too. any ideas on what can be wrong?
