| General > General Technical Chat |
| FFT window function calculation detail |
| << < (8/8) |
| electr_peter:
I get much better picture of weights effect on FFT due to such insightful comments. Essential point - always check actual implementation for frequency response, as it may not be what you expect. With checks it is clear which triangle window implementation works better with odd/even \$N\$. (Implicit) assumption of even \$N\$ has strong implications for window weights, most sources do not discuss evenness topic if at all. Theory and practice differ slightly, this was the root of my misunderstanding. Formula above (lacks normalization by \$2 \pi i k\$) is good for quick approximation, but is has some "wavy" amplitude error, especially near \$k = 0.5\$ - it shows lower amplitude (check rectangle window \$N = 21\$ vs padded FFT). Compare this function with padded FFT as well before use. |
| SiliconWizard:
Yep. And use a decent tool for verifying your assumptions. Can be Matlab, Scilab, Octave, Julia, whatever you're comfortable with. |
| rhb:
I suggest you read the Wikipedia page you linked to. According to it Rader is computing the DFT by padding with zeros to a length that can be computed quickly. Hale's innovation, if original to Dave, is in testing series lengths near prime lengths to find the quickest longer transform. The FFTW "wisdom" file is machine and compiler specific. "Geologists" rarely, if ever, use the FFT. Geophysicists, however, are another matter. The entire field of digital signal processing started at MIT in Norbert Weiner's Geophysical Analysis Group. And geophysicists were doing FFTs *before* Cooley-Tukey appeared. They never published because they thought the interchange of order of operation was obvious to anyone using a desk calculator and pad of paper. In seismic processing a "window" is the length of the non-zero portion of a time series. Zeros on the end don't count. You can easily verify this by comparing the transform of 99 ones followed by 925 zeros and 100 ones followed by 924 zeros. A very common task in seismic processing is to estimate the variation of the signal spectrum over a series of windows in time. Probably the most common window shape for that is a raised cosine edge boxcar. If you want to define "window" differently you are free to do so. I've worked in 7 seismic processing environments. Each one has required about 6 months to learn the local jargon well enough to know who is bullshitting and who actually knows what they are doing. I went into reflection seismology because there were no jobs for a hard rock petrologist. So I switched from light to elastic waves and from hard rocks to soft rocks. Reg |
| Nominal Animal:
In the hopes that I might post something useful that is not derailed by misleading advise from others, I'll sum up the important details here. First, the discrete windowed Fourier transform is defined as $$X_k = \sum_{n = 0}^{N - 1} w_n x_n e^{-2 i \pi k n / N} \tag{DFT}\label{DFT}$$ where \$N\$ is the size of the Fourier transform, \$w_n\$ are the window coefficients, \$x_n\$ are the (time-domain) samples, and \$X_k\$ are the frequency domain complex values, with \$\lvert X_k \rvert\$ being the amplitude, and \$\arg X_k\$ the phase. (Note that if you consider \$w_n x_n\$ and \$X_k\$ \$N\$-dimensional vectors \$\mathbf{x}\$ and \$\mathbf{X}\$, then the above can be written in linear algebra terms as \$\mathbf{X} = \mathbf{F} \mathbf{x}\$, where \$\mathbf{F}\$ is an \$N \times N\$ square matrix, \$\mathbf{F}_{j k} = \omega^{j k} / \sqrt{N}\$, where \$\omega = e^{-2 \pi i / N}\$ is the \$N\$'th complex root of unity, with \$\omega^n = \omega^{n \mod N}\$. \$\mathbf{F}\$ is often called the DFT matrix. One way to achieve a Fast Fourier transform is to exploit symmetries in this matrix.) The frequency response of a discrete window function is obtained by setting \$x_n = 1\$, i.e. $$X_k = \sum_{n = 0}^{N - 1} w_n e^{-2 i \pi k n / N}$$ This is what you want to use when evaluating the response of your discrete window, using the same machinery you use for the DFT calculation. You can evaluate it at an arbitrary frequency \$f\$, $$X(f) = \sum_{n \lvert w_n \ne 0} w_n e^{-2 i \pi f n}$$ such that \$X_k = X(k/N)\$; and similarly for the windowed Fourier transform, $$X(f) = \sum_{n \lvert w_n \ne 0} w_n x_n e^{-2 i \pi f n} \tag{WFT}\label{WFT}$$ Discrete Fourier Transforms, including all Fast Fourier Transforms, on real data calculate \$\eqref{DFT}\$ for frequencies \$f = k/N\$ for \$k = 0 \dots \lfloor N/2 \rfloor\$. \$\eqref{DFT}\$ and \$\eqref{WFT}\$ are exactly equivalent. This shows two things * \$N\$, the size of the FFT, only affects the frequency binning, the frequencies \$f = k/N\$ where the discrete transform \$\eqref{DFT}\$ samples the continuous transform \$\eqref{WFT}\$. * Only samples at nonzero window coefficients \$w_n\$ contribute to the frequency response. In other words, window coefficients \$w_n = 0\$ are useless, unless they are only used to pad the FFT to a desired larger size.When we have a symmetric window, \$w_j = 0\$ for \$j \lt 0\$ and \$j \gt J\$, and \$w_j = w_{J - j}\$, there may be a difference in the frequency response \$X(f)\$ depending on whether \$J\$ is odd or even; i.e. whether the peak \$w_j = \max(w)\$ consists of an odd or an even number of coefficients. This difference depends on the exact window function shape. For many, like the Hann window, there is no real difference at all. To capture all information in a signal using windowed DFT, we need to overlap the windows, by generating a new DFT using samples \$x_{n+L}\$. The exact amount of overlap (or the amount of "time" covered by a single DFT window) depends on the windowing function, not on the FFT size \$N\$. For a rectangular window \$w_j = 1\$ for \$j = 0 \dots J-1\$, and \$w_j = 0\$ for \$j \lt 0\$ and \$j \ge J\$, no overlap is needed. For many windowing functions, \$L < J\$. As shown earlier, the zero coefficients \$w_j\$ do not contribute to the frequency spectrum at all. Some windowing functions use \$w_0 = 0\$ because it is convenient, and define the width of the window by including this useless zero coefficient –– again, useless, because it does not contribute. However, they often define the required overlap or \$L\$ in terms of their definition of the width of the window. It is very common to define \$L = \lfloor J/2 \rfloor\$ (half the width of the window, including leading and/or trailing zeroes in the window function coefficients, with symmetry defined as \$w_j = w_{J-j}\$. Even if you were to omit all leading and trailing zero window coefficients \$w_j\$, your window overlap \$L\$ does not change. This means that in terms of the required overlap \$L\$, and the size of the FFT \$N\$, the most useful way to define the window function width \$W\$ is in terms of distance between first and last nonzero coefficients, i.e. \$W = j_{\max} - j_{\min} + 1\$, where \$w_{j_{\min}}\$ is the smallest-index nonzero coefficient, and \$w_{j_{\max}}\$ is the largest-index nonzero coefficient. Then, \$W\$ is the number of samples needed per FFT, and \$L \le W \le N\$. When \$N \gt W\$, the data is simply padded with zeroes. If your sample rate is \$S\$ samples per second, you'll need to \$W S / L\$ multiplications per second to window the incoming signal, and compute \$S/L\$ FFTs of size \$N\$ each second. The only issue with this definition of \$W\$ is that the way window function overlaps are defined in the literature, you cannot rely on \$L / W/2\$, and must instead calculate \$L\$ based on whatever definitions are commonly used for that window function. With regards to "Fast Fourier Transforms", the definition \$\eqref{DFT}\$ for \$k = 0 \dots N-1\$ involves \$(N-1)^2 \approx O(N^2)\$ multiplications. Any method that computes the same in the order of \$O(N \log N)\$ operations, qualifies as a "Fast Fourier Transform". Since the \$O\$ notation does not contain any coefficients, just the asymptotic behaviour: not how fast an approach is, but how it scales when \$N\$ increases. Fastest Fourier Transforms are based on power-of-two \$N\$, but there are almost as fast implementations for composite \$N\$ with small factors (like say \$N = 3888 = 2^4 3^5\$). Rader's and Bluestein's Fast Fourier Transforms can handle prime \$N\$, but tend to be slower than the fastest algorithms for similar \$N\$ that are composite or powers of two, by some constant factor depending on the implementation details, but typically on the order of 2× to 10×. Both also rely on convolution: Rader's calculates the exact DFT using the convolution of two DFTs of size \$N-1\$ (which must be composite if \$N\$ is prime), and Bluestein's by using the convolution of two DFTs of size \$N_2 \ge 2 N + 1\$ where \$N_2\$ is composite with small prime factors, or a power of two. Even though Bluestein's pads the two DFTs with zeroes, it still calculates the exact DFT of the original \$N\$, because it is a convolution of two DFTs of suitably defined subsequences. In the worst case, \$N\$ being a Mersenne prime, Bluestein's does two DFTs of size \$2 N + 2\$, followed by \$N+1\$ complex multiplications. It is crucial to understand that both Rader's and Bluestein's calculates the exact DFT for the given \$N\$, per the discrete Fourier transform formula, even though they use convolutions of discrete Fourier transforms at other sizes to do so. They are no trick or approximation, they are mathematically exact. So, when one needs a frequency spectrum with frequency bins at specific inverse prime intevals, \$N\$ prime, Rader's (if \$N-1\$ has only small prime factors) or Bluestein's (if \$2 N + 2\$ is a power of two, or when \$N-1\$ has large prime factors) will do that for you, scaling just like all FFTs do as \$N\$ increases, using at most a constant times the number of operations compared to if \$N\$ was a power of two (and that constant depending on the implementation, typically on the order of 2 to 10). |
| Nominal Animal:
As a proof of my advice, please take a look at my Window Function Spectrum Analysis page. It is a self-contained HTML page (so you can download it and open it locally; no network connection needed). For a practical example, put coefficients 0.25 0.50 0.75 1.00 0.75 0.50 0.25 0.00 in the first (upper/red) text area, and coefficients 0.25 0.50 0.75 1.00 1.00 0.75 0.50 0.25 int the second (lower/blue) text area, and click on the Calculate button. This shows you the difference between odd and even-width triangle filters. The page shows you the frequency spectrum you should expect if you take an FFT of those coefficients (real, imaginary format). Feel free to compare the shown values to whatever FFT implementation you use, with the FFT window size the same as the number of coefficients input. In the plot, vertical axis is normalized to the sum of the coefficients for each spectrum. The vertical lines show you where the FFT samples the spectrum (if your FFT window size were the same as the number of coefficients). If you add eight zeros to both text areas, instead of five vertical lines, you'll see nine. Padding or prepending zeros to the coefficients only changes the FFT size (where the vertical lines are shown), and you see how that affects the sampling of the continuous spectrum. The horizontal axis is linear in frequency, with DC at the left edge, and half the sample rate at the right edge. |
| Navigation |
| Message Index |
| Previous page |