https://lcamtuf.substack.com/p/is-the-frequency-domain-a-real-place [https] lcamtuf's thing SubscribeSign in Share this post [https] Is the frequency domain a real place? lcamtuf.substack.com Copy link Facebook Email Note Other Is the frequency domain a real place? The discrete Fourier transform is vital to electronics, signal processing, and radio -- but we're reading too much into what it means. Apr 07, 2024 9 Share this post [https] Is the frequency domain a real place? lcamtuf.substack.com Copy link Facebook Email Note Other 7 Share In an earlier article on the Fourier transform, I talked about the frequency domain -- a wondrous mathematical place where complex signals are transmuted into the amplitudes and phases of sine waveforms. The frequency domain allows us to perform all kinds of signal processing tricks that seem nearly impossible to pull off when we stare at the data in its most straightforward form -- that is, in the time domain. [https] A waveform (top) and a frequency view (bottom) of "Girl in Blue". At the end of that deep dive, I left one question unanswered: how real is this frequency place, anyway? The discrete Fourier transform (DFT) plays a central role in communications and signal processing -- but does it reveal some deeper, unseen truth about the universe? For example, do square waves exist at all? After all, the transform turns them into a series of odd-numbered sine harmonics -- and this model somehow predicts the behavior of electronic circuits in real life. Today, I'd like to knock the Fourier transform off the pedestal. To be sure, sine waves are ubiquitous in nature and they are well-suited for a number of tasks when employed as an analytical tool. That said, it's eminently possible to construct other well-behaved frequency domains that play by different rules -- including one where only square waves are real, and everything else is just harmonics. Revisiting discrete cosine transform (DCT) To get started, let's circle back to discrete cosine transform -- a simplified, real-numbers-only version of DFT. From the earlier article, you might recall the following DCT-II formula: \(F_k = \sum\limits_{n=0}^{N-1} s_n \cdot cos ( \pi k { n + \frac12 \ over N} )\) The operation boils down to taking a series of input values (s[n][ ]-- say, audio samples), multiplying each by the value of a particular cosine expression, and then summing the result to get the magnitude reading for a particular frequency bin (F[k]). The construct at the heart of this algorithm is the cos() expression that generates a sine wave with a frequency corresponding to the number of the current DCT bin. This is known as the basis function; we can abstract it away and rewrite the formula as: \(F_k = \sum\limits_{n=0}^{N-1} s_n \cdot B(k,n)\) In this generalized notation for a frequency-domain transform, B(k, n) returns some sort of a multiplier based on the values of k and n. Software engineers might find it intuitive to think about B(k, n) as a lookup array. In fact, let's calculate that array -- a matrix in the math parlance -- for N=16: [https] DCT-II basis function plot for N=16. If you remember the operation of DCT, this visual should be easy to parse. The first row (k=0) corresponds to the DC component -- i.e., a cosine "running" at 0 Hz, perpetually stuck at +1.00. The next row contains a cosine completing one half of a cycle, going from +1.00 to -1.00; this is followed by a full cycle at k=2, a cycle and a half at k=3, and so on. Into the square-verse! So, how does one go about building a new basis function that splits signals not into sine frequencies, but into square waves? Foreshadowing a bit, the answer appears almost ridiculously simple: [https] The Walsh matrix for N=16. This is known as the Walsh matrix. It essentially consists of square waves running at different speeds, although with some complex symmetries thrown in. And yes: every multiplier is just a +1 or a -1, so the computation boils down to flipping some signs in the input data and then summing the results. The matrix looks fairly trivial, but its design is involved. To capture all frequency and phase information, the rows have increasing sequency -- that is, each next row has just one more sign flip than the one before. Further, the pattern is carefully engineered to ensure orthogonality -- the fragile input-output symmetry that allows seamless conversions back and forth between the frequency representation and the original time-series data. Because of its complex structure, the most practical way to construct the Walsh matrix is to start by generating something known as the Hadamard matrix. For N=16, this intermediate matrix looks the following way: [https] The Hadamard matrix for N=16. Huh? Getting to know Mr. Hadamard At a glance, the plot looks chaotic, but it's simply a reordering of the Walsh layout. For example, rows #1 and #8 are swapped; the same goes for #1 and #15. And unlike Walsh, this fractal-esque pattern is actually fairly easy to create from scratch. To get the ball rolling with Hadamard, we start with the following trivial 1x1 array: \(H_0 = \begin{bmatrix} +1 \end{bmatrix}\) From there, we iteratively extend it by taking the array generated in the previous step (H[n-1]) and tiling four copies of it on a grid with twice the original dimensions. The first three copies are verbatim, and the final one -- bottom right -- has all the signs flipped. The mathematical notation for this extension is: \(H_{n} = N_{n-1} \otimes \begin{bmatrix} +1 & +1 \\ +1 & -1 \end {bmatrix}\) The fancy operator ([?]) is known as the Kronecker product, but it's really just glorified copy-and-paste. The first extension is: \(H_1 = \begin{bmatrix} \begin{array}{c | c } +1 & +1 \\ \hline +1 & -1 \end{array}\end{bmatrix} \) Another iteration gives us: \(H_2 = \begin{bmatrix} \begin{array}{c c | c c} +1 & +1 & +1 & +1 \\ +1 & -1 & +1 & -1 \\ \hline +1 & +1 & -1 & -1 \\ +1 & -1 & -1 & +1 \ end{array} \end{bmatrix}\) ...and so on. On a computer, the matrix can be computed by following this tiling algorithm, but there's a cute bitwise arithmetic trick we can employ instead: as it turns out, the value of the Hadamard function at a particular cell can be divined by calculating x & y and then checking if the number of bits set in the result is even or odd. The following C code does just that, and then displays an 8x8 Hadamard matrix on the screen: #include #include #define HD_LEN 8 int8_t hadamard(uint32_t x, uint32_t y) { return (__builtin_popcount(x & y) % 2) ? -1 : 1; } int main() { for (uint32_t y = 0; y < HD_LEN; y++) { for (uint32_t x = 0; x < HD_LEN; x++) printf("%+d ", hadamard(x, y)); putchar('\n'); } } In principle, this matrix is sufficient to construct a frequency-domain transform. That said, its ordering of the frequency bins is unintuitive, so let's try to fix that before we leave. Enter Mr. Walsh To turn the Hadamard matrix in the nicely-ordered flavor showcased earlier, we need to sort the rows based on their sequency. I'm not aware of an algorithm more elegant than counting the number of zero crossings, but in any case, this implementation is pretty fast -- and prints a sorted 8x8 array on the screen: #include #include #include #define HD_LEN 8 int8_t hadamard(uint32_t x, uint32_t y) { return (__builtin_popcount(x & y) % 2) ? -1 : 1; } int8_t walsh_array[HD_LEN][HD_LEN]; void precompute_walsh() { int8_t tmp[HD_LEN]; for (uint32_t y = 0; y < HD_LEN; y++) { uint32_t sign_cnt = 0; int8_t prev_val = 1; for (int x = 0; x < HD_LEN; x++) { tmp[x] = hadamard(x, y); if (tmp[x] != prev_val) { sign_cnt++; prev_val = tmp[x]; } } memcpy(walsh_array[sign_cnt], tmp, HD_LEN); } } int main() { precompute_walsh(); for (uint32_t y = 0; y < HD_LEN; y++) { for (uint32_t x = 0; x < HD_LEN; x++) printf("%+d ", walsh_array[y][x]); putchar('\n'); } } Equipped with this, we can gut the DCT implementation and make a "discrete square transform" and its inverse: ...previous code here, sans main()... void sqft(double* out_buf, double* in_buf, uint32_t len) { for (uint32_t bin_no = 0; bin_no < len; bin_no++) { double sum = 0; for (uint32_t s_no = 0; s_no < len; s_no++) sum += in_buf[s_no] * walsh_array[s_no][bin_no]; out_buf[bin_no] = sum; } } void isqft(double* out_buf, double* in_buf, uint32_t len) { for (int s_no = 0; s_no < len; s_no++) { double sum = 0; for (int bin_no = 0; bin_no < len; bin_no++) sum += in_buf[bin_no] * walsh_array[bin_no][s_no]; out_buf[s_no] = sum / len; } } Technically, this is called the Walsh-Hadamard transform (WHT), but never mind. Let's confirm that it works: ...continuing from previous snippet... void print_buf(const char* prefix, double* buf, uint32_t len) { printf("%s", prefix); for (uint32_t i = 0; i < len; i++) printf("%+.2f ", buf[i]); putchar('\n'); } int main() { double in[HD_LEN] = { 1, 1, 1, 1, 5, 5, 5, 5 }; double sq_freq[HD_LEN], out[HD_LEN]; precompute_walsh(); sqft(sq_freq, in, HD_LEN); isqft(out, sq_freq, HD_LEN); print_buf("Input : ", in, HD_LEN); print_buf("SQFT : ", sq_freq, HD_LEN); print_buf("ISQFT : ", out, HD_LEN); return 0; } The input is a square-wave-ish sequence of numbers: 1 1 1 1 5 5 5 5. The output from DFT or DCT would be a bunch of harmonics across multiple frequency bins: DCT : +24.00 -10.25 -0.00 +3.60 +0.00 -2.41 -0.00 +2.04 In contrast, the frequency-domain representation generated by the program shows non-zero components only in F[0] (DC) and F[1], confirming that we have an algorithm that deconstructs the signal into pure square waves: SQFT : +24.00 -16.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 Finally, we can verify that the inverse function -- isqft() -- transforms the frequency domain back to what we started with: ISQFT : +1.00 +1.00 +1.00 +1.00 +5.00 +5.00 +5.00 +5.00 The Walsh-Hadamard transform, being computationally efficient and well-suited for certain types of data, finds uses here and there. The point isn't that it needs to be used more; it's just that discrete Fourier doesn't have a monopoly on truth. [ ] Subscribe For more articles on electronics, click here. 9 Share this post [https] Is the frequency domain a real place? lcamtuf.substack.com Copy link Facebook Email Note Other 7 Share 7 Comments [https] [ ] Megan Jones 3 hrs ago This is a great article, but for future ones, could you use [https] Rust instead of C for any code? Rust is now the preferred programming language for code examples. Expand full comment Reply Share 3 replies by lcamtuf and others lcamtuf 15 hrs agoAuthor PS #1: Apologies to all the folks who signed up for this after the xz article and didn't realize what's coming! author PS #2: Here's a person experimenting with Walsh-Hadamard for image compression: http://rotormind.com/blog/2019/ hadamard-days-night/ Expand full comment Reply Share 5 more comments... Top Latest Discussions No posts Ready for more? [ ] Subscribe (c) 2024 lcamtuf Privacy [?] Terms [?] Collection notice Start WritingGet the app Substack is the home for great writing Share Copy link Facebook Email Note Other This site requires JavaScript to run correctly. Please turn on JavaScript or unblock scripts