GLSL FFT Issue

[eli-jordan]

I'm getting some weird results when attempting to port glsl-fft to processing

• fragment shader source
• vertex shader source this is auto generated by the library I'm using.
• driver source i.e. the thing that generates the configuration of the passes.
• test cases

The only notable difference I can see between the two implementations is the way the vertexes are defined. In glsl-fft a triangle is drawn, whereas in my implementation a quad is drawn. I don't see how this could make a difference.

One thing that I thought may have an impact on this, is the different in the coordinate systems for gl_FragCoord and texture coordinates. AFAIK gl_FragCoord has the origin at the bottom left corner, whereas texture coordinates have the origin at the upper left corner. Perhaps regl is doing some adjustment to the coordinate system that is not being done in my solution?

Example 1

Input (all real components):

0, 0, 0, 1
0, 0, 0, 0
0, 0, 0, 0
0, 0, 0, 0

Expected:

 + 1.0 + 0.0i,  + 0.0 + 1.0i,  - 1.0 + 0.0i,  + 0.0 - 1.0i, 
 + 1.0 + 0.0i,  + 0.0 + 1.0i,  - 1.0 + 0.0i,  + 0.0 - 1.0i, 
 + 1.0 + 0.0i,  + 0.0 + 1.0i,  - 1.0 + 0.0i,  + 0.0 - 1.0i, 
 + 1.0 + 0.0i,  + 0.0 + 1.0i,  - 1.0 + 0.0i,  + 0.0 - 1.0i,

Actual:

 + 1.0 + 0.0i,  + 0.0 - 1.0i,  - 1.0 + 0.0i,  + 0.0 + 1.0i, 
 + 1.0 + 0.0i,  + 0.0 - 1.0i,  - 1.0 + 0.0i,  + 0.0 + 1.0i, 
 + 1.0 + 0.0i,  + 0.0 - 1.0i,  - 1.0 + 0.0i,  + 0.0 + 1.0i, 
 + 1.0 + 0.0i,  + 0.0 - 1.0i,  - 1.0 + 0.0i,  + 0.0 + 1.0i,

Notice that column 1 and 3 (0 indexed) are swapped in the actual output as compared to the expected output.

Example 2

Input (all real components):

1, 0, 0, 0
0, 0, 0, 0
0, 0, 0, 0
0, 0, 0, 0

Expected:

 + 1.0 + 0.0i,  + 1.0 + 0.0i,  + 1.0 + 0.0i,  + 1.0 + 0.0i, 
 + 0.0 + 1.0i,  + 0.0 + 1.0i,  + 0.0 + 1.0i,  + 0.0 + 1.0i, 
 - 1.0 + 0.0i,  - 1.0 + 0.0i,  - 1.0 + 0.0i,  - 1.0 + 0.0i, 
 + 0.0 - 1.0i,  + 0.0 - 1.0i,  + 0.0 - 1.0i,  + 0.0 - 1.0i,

Actual:

 + 1.0 + 0.0i,  + 1.0 + 0.0i,  + 1.0 + 0.0i,  + 1.0 + 0.0i, 
 + 0.0 - 1.0i,  + 0.0 - 1.0i,  + 0.0 - 1.0i,  + 0.0 - 1.0i, 
 - 1.0 + 0.0i,  - 1.0 + 0.0i,  - 1.0 + 0.0i,  - 1.0 + 0.0i, 
 + 0.0 + 1.0i,  + 0.0 + 1.0i,  + 0.0 + 1.0i,  + 0.0 + 1.0i,

Notice that row 1 and 3 (0 indexed) are swapped in the actual output as compared to the expected output.

[rreusser]

Interesting… it looks… to me… like both are exactly the complex conjugate. I wonder if this generalizes. If so, I wouldn't necessarily know the cause, the the solution would sure be easy…

Also, regarding the triangle, the motivation is simply that it (might) have better cache performance than a quad: https://michaldrobot.com/2014/04/01/gcn-execution-patterns-in-full-screen-passes/ It should otherwise be exactly identical.

[rreusser]

Ah, and regarding coordinate systems, I'm not sure I'd say the texture coordinates origin is in the top-left. I visualize normalized device coords and texture coords the same, except one is in [-1, 1] x [-1, 1] while the other is in [0, 1] x [0, 1].

[eli-jordan]

Thanks for taking a look!

Regarding the coordinate systems, I was going off this image

which is from here https://gamedevelopment.tutsplus.com/tutorials/a-beginners-guide-to-coding-graphics-shaders--cms-23313

and that when I have a very simple shader that does something like

void main(void) {
   // resolution is a vec2 with the screen width and height
   vec2 pos = gl_FragCoord.xy / resolution;
   gl_FragColor = texture(lena, pos);
}

I get an upside down image.

---

Regarding the complex conjugate idea. I don't think thats right, because there are several entries in the result matrix that actually match the expected result, so taking the conjugate of those would result in an incorrect answer.

For example in Example 1 above the expected and actual result at (0, 0) is + 1.0 + 0.0i along with the rest of the entries in the 0th and 3rd columns.

---

I'll keep digging. Let me know if you have other ideas :)

[rreusser]

Ah, I guess that argument works. I've always instead thought of textures as upright and if something accidentally gets flipped, I prefer to flip the input on texture creation (regl has a flipY flag for this purpose).

Also, you are probably right about complex conjugate not being the right answer, though to be clear, in the above examples, I think the complex conjugate would solve it—though it probably doesn't generalize.

[rreusser]

And just in case anything was lost in translation, my implementation was reworked from this code, and then fairly carefully tested against ndarray-fft: https://david.li/filtering/

[eli-jordan]

Also, you are probably right about complex conjugate not being the right answer, though to be clear, in the above examples, I think the complex conjugate would solve it—though it probably doesn't generalize.

Ah, yes you are right. I missed that -0 and +0 are the same thing :-S

Also, now that I think about it, the geometric interpretation of the conjugate is flipping the i axis. Which makes me think there is something up with the coordinates again.

I'll try calculating the conjugate, and running it through a bunch more examples and see if that works.

[eli-jordan]

It seems to work for any real valued input with dimensions <= 8 i.e. 2,4,8
As soon as I add imaginary values in the input or make the width or height >= 16 it fails.

Thats an odd result for sure.

[rreusser]

FWIW I find that 4x4 is a very effective size for testing, but think 0/1 might be too degenerate. Sometimes even just a serial counter e.g. [[0, 1, 2,3], [4, 5, 6, 7], ...] is sufficient to break the symmetry and ensure a meaningful signal.

[eli-jordan]

Yeah I ran a test with a few thousand randomly generated inputs with various combinations of dimension 2,4 and 8 with the real values between -100 and 100 and 0 imaginary value. Applying the complex conjugate worked every time.

[rreusser]

Isn't the complex conjugate equivalent to time reversal (in 1d)? A flipped input dimension would be consistent with that.

[eli-jordan]

Right, thats what I was thinking when I said

Also, now that I think about it, the geometric interpretation of the conjugate is flipping the i axis. Which makes me think there is something up with the coordinates again.

Also noticing that the size issue was just a matter of precision. My tolerance was 0.001 which was too small apparently. Making that larger means that my tests pass at larger sizes up to 128x128. Though the precision issues get larger as the size goes up.

I'll try flipping the y dimension on the input, and see if that works out.

[eli-jordan]

No dice, when trying various transformations of the input I still get wrong results.

Below is the result from inverting the y axis.
Note, I didn't change anything in the shader, just inverted the input data at the beginning of the passes.

Input:

 1,  2,  3,  4
 5,  6,  7,  8
 9, 10, 11, 12
13, 14, 15, 16

Inverted to:

13, 14, 15, 16 
 9, 10, 11, 12
 5,  6,  7,  8
 1,  2,  3,  4

Expected:

 + 136.0 + 0.0i,       - 7.9999995 + 8.0i,  - 8.0 + 0.0i,  - 8.0 - 8.0i, 
 - 31.999998 + 32.0i,  + 0.0 + 0.0i,        + 0.0 + 0.0i,  + 0.0 + 0.0i, 
 - 32.0 + 0.0i,        + 0.0 + 0.0i,        + 0.0 + 0.0i,  + 0.0 + 0.0i, 
 - 32.0 - 32.0i,       + 0.0 + 0.0i,        + 0.0 + 0.0i,  + 0.0 + 0.0i, 

Actual:

 + 136.0 + 0.0i,   - 8.0 - 8.0i,  - 8.0 + 0.0i,  - 8.0 + 8.0i, 
 + 32.0  + 32.0i,  + 0.0 + 0.0i,  + 0.0 + 0.0i,  + 0.0 + 0.0i, 
 + 32.0  + 0.0i,   + 0.0 + 0.0i,  + 0.0 + 0.0i,  + 0.0 + 0.0i, 
 + 32.0  - 32.0i,  + 0.0 + 0.0i,  + 0.0 + 0.0i,  + 0.0 + 0.0i, 

---

I think there is something else going on, because I start getting the wrong answer when the input has non-zero imaginary components even when applying the conjugate after reading the data.

[eli-jordan]

I figured it out 🥳🕺
After all that, it was such a dumb mistake. I had to convert from boolean to int as interpreted in glsl. 0 == true and 1 == false I had it the other way around https://github.com/eli-jordan/generative-art/blob/master/processing-app/src/main/java/glslfft/GlslFft.java#L149

[eli-jordan]

Ok, next issue is non-square matrices. Its never easy 😭

[rreusser]

I didn't see your previous comment at the time, but glad it worked! FWIW apart from initial pass setup, I didn't need any particular changes for non-square matrices as long as they were power-of-two, but ymmv 😄

[eli-jordan]

Thanks!
Yeah, I'm guessing I mixed up width and height variables somewhere. I just have to go through and find where!

[eli-jordan]

I've been looking for the error that cases failure on non-square input, but with no success. If you have any ideas where to look, I'd very much welcome any input!

I'm running on a 4x2 input. Below is the input and output for each pass.

Input

1, 2, 3, 4
5, 6, 7, 8

Pass configurations:

[forward=true, input=glslfft.GlslFft$FftBuffer@38e2fda7, output=glslfft.GlslFft$FftBuffer@337bd16c, normalization=1.0, horizontal=true, subtransformSize=2.0, resolution=(0.25, 0.5)]
[forward=true, input=glslfft.GlslFft$FftBuffer@337bd16c, output=glslfft.GlslFft$FftBuffer@7b230bdc, normalization=1.0, horizontal=true, subtransformSize=4.0, resolution=(0.25, 0.5)]
[forward=true, input=glslfft.GlslFft$FftBuffer@7b230bdc, output=glslfft.GlslFft$FftBuffer@573bb2eb, normalization=1.0, horizontal=false, subtransformSize=2.0, resolution=(0.25, 0.5)]

------------- Pass 0: Inputs -------------
 + 1.0 + 0.0i,  + 2.0 + 0.0i,  + 3.0 + 0.0i,  + 4.0 + 0.0i, 
 + 5.0 + 0.0i,  + 6.0 + 0.0i,  + 7.0 + 0.0i,  + 8.0 + 0.0i, 

------------- Pass 0: Outputs -------------
 + 6.0 + 0.0i,  + 8.0 + 0.0i,  + 10.0 + 0.0i,  + 12.0 + 0.0i, 
 - 4.0 + 0.0i,  - 4.0 + 0.0i,  - 4.0 + 0.0i,  - 4.0 + 0.0i, 

------------- Pass 1: Inputs -------------
 + 6.0 + 0.0i,  + 8.0 + 0.0i,  + 10.0 + 0.0i,  + 12.0 + 0.0i, 
 - 4.0 + 0.0i,  - 4.0 + 0.0i,  - 4.0 + 0.0i,  - 4.0 + 0.0i, 

------------- Pass 1: Outputs -------------
 + 2.0 + 0.0i,  + 4.0 + 0.0i,  + 6.0 + 0.0i,  + 8.0 + 0.0i, 
 - 4.0 + 4.0i,  - 4.0 + 4.0i,  - 4.0 + 4.0i,  - 4.0 + 4.0i, 

------------- Pass 2: Inputs -------------
 + 2.0 + 0.0i,  + 4.0 + 0.0i,  + 6.0 + 0.0i,  + 8.0 + 0.0i, 
 - 4.0 + 4.0i,  - 4.0 + 4.0i,  - 4.0 + 4.0i,  - 4.0 + 4.0i, 

------------- Pass 2: Outputs -------------
 + 8.0 + 0.0i,  - 4.0 + 0.0i,  + 12.0 + 0.0i,  - 4.0 + 0.0i, 
 - 8.0 + 8.0i,  + 0.0 + 0.0i,  - 8.0 + 8.0i,  + 0.0 + 0.0i,

Expected result

 + 36.0 + 0.0i,  - 3.9999998 + 4.0i,  - 4.0 + 0.0i,  - 4.0 - 4.0i, 
 - 16.0 + 0.0i,  + 0.0 + 0.0i,  + 0.0 + 0.0i,  + 0.0 + 0.0i,