6 ms·
There’s an alternate way to save space: the discrete Hartley transform. https://news.ycombinator.com/item?id=27386319 https://news.ycombinator.com/item?id=2738
by sillysaurusx 4y ago
There’s an alternate way to save space: the discrete Hartley transform.
https://news.ycombinator.com/item?id=27386319 https://news.ycombinator.com/item?id=27386319
https://twitter.com/theshawwn/status/1400383554673065984 https://twitter.com/theshawwn/status/1400383554673065984
dht(x) = fft( (1 + 1j)*x ).real / area(x)**0.5
area(x) = np.prod(np.shape(x))
The dht is its own inverse, which is a neat trick. No need for an inverse fft. In other words, dht(dht(x)) == x.
Notice that the output is the real component of the fft. That means the storage required is exactly identical to the input signal. And since the input signal is copied to the real and imaginary components (1 + 1j)*x, the intermediate representation size is also equal to the input signal size. (Just take sin and cos of the input signal.)
The dht has almost all the same properties of the fft. You can blur an image by taking a hamming window and multiplying. You can resize an image by fftshift + pad with zeros + fftshift. Etc.
In practice the dht seems to have all the benefits of fft with very few downsides. I’m amazed it isn’t more popular.
- enriquto 4y agoThe convolution theorem is much uglier for the hartley transform. I guess this is the main reason. You can see the Fourier transform as a simple variant of the Hartley transform where the linear filtering can be done directly in the frequency domain, without all the fuss of the hartley transform convolutions. Another advantage of Fourier with respect to Harley is that derivatives become products with polynomials. Thus the Fourier transform transforms linear pde into trivial algebraic equations. For the Hartley transform, you get more complicated functional equations, so that it's not as useful. > You can blur an image by taking a hamming window and multiplying. It doesn't seem like you can? (Unless your image has some symmetries).
- sillysaurusx 4y agoHere's an image: https://i.imgur.com/hNCwisy.png https://i.imgur.com/hNCwisy.png Here's fftshift(dht2(img)): https://i.imgur.com/EB6yIb4.png https://i.imgur.com/EB6yIb4.png Here's a hamming window: https://i.imgur.com/NsGKb0X.png https://i.imgur.com/NsGKb0X.png Here's a hamming window raised to power of 20: https://i.imgur.com/ejCviD7.png https://i.imgur.com/ejCviD7.png Here's the result of multiplying that with the transformed image: https://i.imgur.com/trHdWZS.png https://i.imgur.com/trHdWZS.png And here's dht(fftshift(that)): https://i.imgur.com/of4bvhv.png https://i.imgur.com/of4bvhv.png Presto, blurry image. No symmetries required. We can increase the blur factor by changing the power from 20 to 200: https://i.imgur.com/yqkR3qb.png https://i.imgur.com/yqkR3qb.png And 2000 (megablur): https://i.imgur.com/7h9patv.png https://i.imgur.com/7h9patv.png The code I used to blur: def blur(img, factor): x = dht2(img) x = fftshift(x) x = x * hamming(img) ** factor x = fftshift(x) x = dht2(x) return x There doesn't seem to be any limitations -- anything you can do with FFT, you can do here. And since the DHT isn't complex, you don't need to deal with phase/amplitude. It's kind of like working with wavelet coefficients in a jpg transform, except wavelets are a lot harder to work with compared to fft signals (which is why people use fft instead of wavelets everywhere). But dht is the best of both worlds -- all the power of fft, none of the hassle of complex values.
- zeristor 4y agoThat’s impressive thanks. I’m we did Fourier transforms as lab exercise using lasers for my Masters, that was one of quite a few amazing and memorable experiments. Another was using impedance to pulse lasers, set me to thinking about impedance as a metaphor in software development. Chemistry A-Level had rates of reaction in Physical chemistry, a few certain key parts being critical in a reaction.
- enriquto 4y ago> There doesn't seem to be any limitations -- anything you can do with FFT, you can do here. There's certainly no limitations if all you want to do is to convolve your image with symmetric kernels. But you asked for advantages of the Fourier transform, if any. Convolution of two non-symmetric images is one of them ;) There's also the fact that derivatives have an easier relationship with the DFT than the DHT. If you are interested in discrete linear pde, that's a big advantage. For signal processing, the main advantage of the DFT is that the amplitudes and phases are cleanly separated. For the applications where you only want the amplitude, it is already there in the DFT, but you have to untangle it from the DHT, because the DHT intermingles the phases and amplitudes. An example: imagine that you have an image affected by periodic noise, and you want to remove this noise. This periodic noise appears in the DFT as a frequency with very large amplitude (and an arbitrary phase, corresponding to the phase of the periodic noise, that you don't really care about). Thus, it is easy to locate and filter-out this periodic noise using the DFT: just find the frequency with large amplitude and set it to zero. Now, for the DHT, the frequency that corresponds to the noise period is not necessarily large (in magnitude). It may even be zero, depending on the phase of the noise! To find that frequency, you have to essentially recover the DFT and compute the amplitudes. Of course the DHT and the DFT are equivalent and there's an easy conversion between both, so you can do exactly the same things with either. The difference is just a matter of convenience. It's certainly useful to know both. If you are scared of complex numbers for some reason (e.g., if your programming language makes it difficult to work with them) then you may prefer the DHT. Otherwise, the DFT neatly separates the phases and amplitudes and it's slightly easier to interpret.
- goldenkey 4y agoTo explain further, HT is unable to compose complex valued functions - it is only a complete system for real valued functions, that is: R -> R. So naturally, it's simpler, because it's less expressive. Let me wow you since HT is not even on the map of exotic or alternative wavelet decompositions. For example, Hermite transform: https://en.wikipedia.org/wiki/Hermite_transform https://en.wikipedia.org/wiki/Hermite_transform or if looking in the discrete space: https://en.wikipedia.org/wiki/Hadamard_transform https://en.wikipedia.org/wiki/Hadamard_transform See: https://en.wikipedia.org/wiki/Orthogonal_functions https://en.wikipedia.org/wiki/Orthogonal_functions I'm amazed that folks think Fourier series is so special when it's just one of many in a deluge of complete function systems.
- kragen 4y agothe fourier transform is very special indeed because of what its basis functions are the eigenvectors of
- goldenkey 4y agoSo you use one special thing, the Laplacian, to connote the specialty of another? I can go ahead and do the same duality for those other transforms and their respective basis functions. Specialness in mathematics is a fruitless and fealty affair - G.H. Hardy
- kragen 4y agothat's not it
- tails4e 4y agoGenuinely curious - How does this save space? If the input was real, the first operation (inside the fft call) makes it complex, which doubles the storage space. Sure the output is real again, but the intermediate storage was 2x the input and output.
- sillysaurusx 4y agoSince the input is copied to the real and imaginary components, you don’t actually need to copy it at all. Think of it like having a pointer (in the C sense) to the real and imaginary image buffers. Set both pointers to the input signal, then run your FFT algo on that. The key is that it minimizes memory bandwidth, not computation. It’s highly cache coherent.
- tails4e 4y agoThanks, so the real and imag components always map to the same value during the processing, i.e. After every butterfly the real and iamg will be the same value? Cause again if the real and imag components ever differ it means more storage needed.
- sillysaurusx 4y agoCorrect, they are identical by definition.
- wdkrnls 4y agoWhat is `area(x)` here? Some kind of numerical integration like the trapezoid rule? Can someone please provide example code which can actually be tested?