Factor of 1/N for halfcomplex-to-real

[Lagrang3]

In this PR halfcomplex-to-real is not the inverse linear map of real-to-halfcomplex, we get instead a factor of N. This is compatible with GSL and FFTW philosophy and the fact also hold for the forward and the backward FFT.

[cosurgi]

So now all the backends behave in the same way, in all the kinds of transforms?

[Lagrang3]

So now all the backends behave in the same way, in all the kinds of transforms?

They have always behaved the same way. The real question is about the behavior of our API. Do we want halfcomplex-to-real to be the inverse of real-to-halfcomplex? In that case DO NOT merge this PR.

If we merge this, then the halfcomplex-to-real is N times the inverse of real-to-halfcomplex, so the user that wants the inverse has to take care of multiplying by 1/N himself.

[cosurgi]

@Lagrang3 ah I see. So this change affects fftw, gsl and bsl all of them in the same way. And we have to decide which way it is. For start it should be the same as complex transform. Next: how about we make it configurable via bool template parameter? This bool would have the same meaning in all parts of the code: real transforms and complex transforms. And depending on it different code would be executed.

We don't have if constexpr in C++11, but it its analog can be constructed via a helper template.

I think a template bool option would be a perfect solution. Because from what I saw in the FFT usage around sometimes one data format is needed, sometimes the other.

CC @ckormanyos

[Lagrang3]

We should be asking why isn't the dft backwards implemented as the inverse of the dft forward. The reason is that the dft backwards is another dft transform but using the opposite sign in the exponential factor. The programming interface is only providing two ways of performing a dft on a complex array. The user should know that the composition of dft backward and dft forward yields the identity times a factor N.

Why is the real case different? If we dft a real sequence a[] we obtain a complex sequence b[] with the halfcomplex symmetry, ie. b*[i] = b[N-i]. Also it can be shown that if you take the dft of a complex sequence that satisfies the halfcomplex symmtry, you obtain a real sequence. These facts are independent of the sign of the exponential factor, but it doesn't mean that the results are independent of the sign. That means we can start with a real sequence and we have two choices, either to dft backward or dft forward, in any case we will obtain a halfcomplex sequence but one is the complex conjugate of the other. On the other hand, if we start with a halfcomplex sequence we can apply either the dft forward or the dft backward, and we will obtain always a real sequence but one will be the reverse indexing of the other, ie. a'[i] = a[N-i].

What FFTW and GSL does, is to provide a real-to-halfcomplex transform that it is always a forward dft, and a halfcomplex-to-real that is always a backward dft. Clearly, if I wanted a real-to-halfcomplex computed with a backward dft, I could just do the complex conjugate of the result; and if I wanted a halfcomplex-to-real computed with a forward dft, I could just do the reverse indexing of the result. Hence they provide two functions instead of 4. In FFTW and GSL, that means that the halfcomplex-to-real is the inverse of the real-to-halfcomplex up to the factor N. But in general, it is not necessarily the case, because of what I said previously, ie. twice a forward dft (instead of one forward and a further backward) starting from a real sequence will obtain the original sequence with reversed indexes, and twice a forward dft starting from a halfcomplex sequence will result in the complex conjugate of the original sequence.

Therefore, I think we should follow this forward/backward approach, not the inverse. And hence proceed with merging this PR.

[Lagrang3]

I was thinking around the lines of abstraction. That we too could present the user with two functions: forward and backward, just like in the complex case. And it will be the difference in the type, to determine if I am performing a real-to-halfcomplex or a halfcomplex-to-real. Notice that from a real sequence, you can only go to the halfcomplex space, and from a halfcomplex you can only obtain a real sequence. We were anyways already discussing the possibility of abstracting the halfcomplex. But we will also need to abstract the real counterpart. But how can we abstract a real sequence, without selecting its container?

[ckormanyos]

follow this forward/backward approach, not the inverse. And hence proceed with merging this PR

Done. Thank you Eduardo!

[cosurgi]

and we will obtain always a real sequence but one will be the reverse indexing of the other, ie. a'[i] = a[N-i].

I imagine this as a graph with four nodes. Two nodes are halfcomplex sequences, two nodes are real sequences, the graph edges connecting the nodes are the transforms.

always a forward dft, and a halfcomplex-to-real that is always a backward dft.

So the FFTW and GSL authors have chosen one directional connection between the nodes.

But in general, it is not necessarily the case,

yes. I see it now. Thanks for discussing these various scenarios. Maybe it makes sense to discuss them in the documentation (maybe even with a simple graph of four connected nodes and the names of the transform function to call written along each graph edge).

But how can we abstract a real sequence, without selecting its container?

I think I saw you once used std::iterator_traits and std::random_access_iterator_tag. But not anymore. Writing in documentation some requirements could help here. So we might not select a container, but assume it is one of those that we need, and the code won't compile if it doesn't satisfy the needs. So that is enforcing container on the user to have some required properties. We could even do a static_assert on a decltype(operator[]) or something along these lines to verify that required methods exist in the container. Such static_assert with a corresponding message will send a clear message to users and would be better than simple compilation error when container doesn't fit the requirements.

[Lagrang3]

Yes. It is a graph with four nodes, as you have said: 2 on the real space connected by the index inversion map and 2 on the halfcomplex space connected by the complex conjugate map. Then each DFT connect a node from the real to the halfcomplex or the halfcomplex to the real.