I'd like to start by noting that for some fixed natural $N$ basis functions for my system will be generated by $f(k,x)$ as defined and explained here or in numerous other sources:

$$f(k,x) = \sqrt2 \cos \left( \frac{2 \pi k }{N}x - \frac \pi 4 \right)$$

So once I let my $k$ be equal to $0,1,2,\ldots,N-2,N-1$ it is obvious that my basis functions are sorted by frequency (in a sense that cosine periods are getting shorter).

However when I use discretization and transfer my functions ($N$ of them) into vectors consisting of $f(k,x)$ values at $x=\overline{0,N-1}$ all the sorting by frequency is lost (in a sense of positive-negative and negative-positive transitions in vector components)!

What is interesting though is that numerical experiments reveal that vector sorting (by frequency) may be obtained by quite elegant rearrangement of initial basis functions i.e. letting $k$ be equal to $0,1,N-1,2,N-2,3,N-3,\ldots$. This observation is nearly worthless without mathematical proof though. So is this a known and already documented property of discrete Hartley transform?

Here is a part of $N=8$ example.

As one can clearly see above $k=6$ function has a higher frequency than $k=4$ function. Yet the opposite can be said about discretized equivalent.

Below is a full representation of $k=\overline{0,7}$ functions. To get the desired sorting by frequency it is needed to let $k$ be $0,1,7,2,6,3,5,4$ just as proposed above.