General > General Technical Chat
FFT window function calculation detail
rhb:
Compute the FFT of the "flat topped" triangle of length N, insert N zeros in the center and back transform. Now you can see the half sample shift.
Padding a zero gives you an N-1 length operator stored in an N length array. It really doesn't matter what you do as long as you understand what you are doing in both domains. If you use N in your calculations rather than N-1 you will get a subtley wrong answer.
@Nominal Animal
The DFT is defined over the semi-closed interval from 0 to 2*Pi. It is sampled at increments of i/N for i=0 to N. For N prime and i=N, N/N = 1. There are no integers between 0 and 1. QED
I have an FFT code I wrote based on one that was circulating at UT Austin in my research group. It is an arbitrary prime factor code in 2 pages of FORTRAN. It will handle *any* length, but run time is N**2 for prime lengths.
It's a *very* good code and I strongly recommend it for general use if you don't want the footprint of FFTW.
https://www.jjj.de/fft/fftpage.html
There's an interesting note from someone who used the algorithm in the 70's while working at IBM. I've got the Ferguson paper, but not the Glassman one.
A bit of trivia from Jorge's FFT page:
I happened upon your FFT page and noticed that you have a generalized radix
FFT algorithm attributed to Glassman (listed as unknown). I don't know if it
is of any interest to you, but I can cite a reference by him which describes
the algorithm. It appeared in the IEEE Transactions on Computers in the
February, 1970 issue. It was written by J. A. Glassman who is listed as
working at Hughes Aircraft in Canoga Park California. The title of the paper
is "A Generalization of the Fast Fourier Transform" and describes the
calculation of a n-point FFT where n is not restricted to be a power of 2.
He provides a Fortran listing of the algorithm as well as the mathematical
development of it.
I used this algorithm and a modification of the code to implement a piece of
test equipment back in the early 1970's when I was working for IBM. I have
never seen another paper by Glassman, but I recently came across a reference
to his generalized radix algorithm which provides an alternate (and perhaps
simpler) derivation. Here is the reference:
W. E. Ferguson, Jr., "A simple derivation of Glassman general-n fast
Fourier transform," Comput. and Math. with Appls., vol. 8, no. 6, pp.
401-411, 1982.
Also, in Report AD-A083 811, NTIS, Dec. 1979.
Wayne Torrey
rhb:
I just realized an important issue may be getting overlooked.
You *must* pad the series to at *least* 2*N to prevent wraparound. The other thing you *must* do is to remove any linear trend so that the first and last sample have the same value so there is no discontinuity. You then *may* want to add back in the linear trend when you return to the original domain. But if there is a discontinuity between the first and last sample it has the effect of a spike in time and creates a uniform bias to all samples in the frequency domain. In all cases you must pay close attention to avoiding unidentified phase shifts. They are lurking everywhere.
If you are resampling by Fourier transform, the input and output series must be the same length . Failure to ensure that will result in phase errors. To make it work you vary the padding of the two series so that N1*dT1 = N2*dT2. This is because of the semi-closed interval over which the DFT is defined.
In general, the speed of an FFT is dependent upon how many factors N has. Dave Hale's prime factor FFT uses benchmarked data to create the "wisdom" file. The test exhaustively tests the performance of all factorizations up to some limit. It then uses those times to decide how much to pad the input. It's in FFTW. Long ago, before it was included the FFTW folks benchmarked it against FFTW and it was faster. Not by much, but still faster. N = 2**n is the most factors you can have at those sizes.
Edit: A few caveats about the DFT. There are 3 scaling factors used, 1/N,N & 1/sqrt(N). Moving data between software packages using different scaling can be a problem if you're not paying attention. There are also 2 sign conventions for the exponent and 2 intervals commonly used, -Pi to Pi and 0 to 2*Pi. The former is the mathematical definition, but the implementations are typically 0 to 2*Pi. F(i*Pi-x) = F((i+2)*Pi-x) which is why the "negative" frequencies are at the end of the transform in descending order. So there are 12 possible implementations of the same concept. This is a common source of confusion and error.
Hopefully this has taken a bit of the mystery out of the FFT. The bottom line is that choosing N so that it is highly composite with lots of small factors provides considerable potential performance gains at very low cost.
Each factor of N is a multiply that can be replaced by an add.
Nominal Animal:
--- Quote from: rhb on December 10, 2022, 04:02:35 pm ---The DFT is defined over the semi-closed interval from 0 to 2*Pi. It is sampled at increments of i/N for i=0 to N.
--- End quote ---
For i=0 to N-1, inclusive, to be exact.
--- Quote from: rhb on December 10, 2022, 04:02:35 pm ---For N prime and i=N, N/N = 1. There are no integers between 0 and 1. QED
--- End quote ---
That has nothing to do with sampling a waveform, and completely irrelevant here. It's like saying "The sky is blue, therefore you are wrong."
What really matters, for prime N, is that \$e^{-2 \pi i / N}\$ is the \$N\$'th complex root of unity. Until you understand this, and why this is so, you really do not grasp Fourier transforms.
You currently surely understand Cooley-Tukey prime factorization based fast discrete Fourier transform, but that is just a small part of all possible fast discrete Fourier transforms. Please, open your mind, before you lead more people astray with your limited understanding and thus misleading advice.
--- Quote from: rhb on December 10, 2022, 04:02:35 pm ---I have an FFT code I wrote based on one that was circulating at UT Austin in my research group. It is an arbitrary prime factor code in 2 pages of FORTRAN. It will handle *any* length, but run time is N**2 for prime lengths.
--- End quote ---
Again, just because you think Cooley-Tukey is the only possible approach, does not make it so.
You already have FFTW installed and are familiar with it, so comparing the run time (after gathering wisdom) for say 1048573 (prime), 1048575 (3×5×5×11×31×41), 1048576 (2**20), 1048583 (prime) DFTs would tell you something important about your own understanding.
Because you didn't bother doing even this, shows that you are one of those who relies on their beliefs, and are incapable of considering they might be wrong, and therefore prefer to lead others astray, instead of checking their own understanding.
Despicable behaviour, rhb.
rhb:
Please explain how the FFT works for prime length series. Where is the computational savings?
As a reminder here's the definition from Brigham.
dft.pdf (12.44 kB - downloaded 56 times.)
What you describe is the behavior of Dave Hale's FFT code. It's likely not computing *any* of the specified series lengths, but a longer series which has been zero padded. Hale is generally considered the best mathematical programmer of his peers from the Stanford Exploration Project. The FFTW guys were quite stunned when Hale pretty consistently beat them on difficult lengths in a major benchmark they published.
I spent most of my career in the oil industry fixing errors in complex numerical codes for which an FFT was just one of many operations and there were often many FFTs, as in millions of 3D FFTs in one code I worked on. In that environment you *really* care about every nuance of the DFT and fast algorithms for computing it.
For a brief period my job was tuning a 3D FFT algorithm on Intel i386 & i860 Hypercubes.
electr_peter:
--- Quote from: Nominal Animal on November 24, 2022, 02:12:16 pm ---... with \$a_0 = 0.5\$ being a Hann window, simplifying (via \$\cos(2 x) = 1 - 2 \sin^2(x)\$) to
\$w_n = \sin^2\left(\frac{\pi n}{N}\right), \quad 0 \le n \le N\$
--- End quote ---
Most window function are defined as above, with \$w_0 = 0\$. Considering what was said above (w.r.t. window length \$L\$ and DFT length \$N\$), zeros are not counted as points inside window. So with \$N\$ data points maximum Hann window length is \$L = N-1\$ and we have one extra zero padded as well (my interpratation). As padded zeros essentially interpolates FFT spectrum points with sinc and do not add any extra info, this leaves me puzzled :-//
Whats is the point of having zeros in window filter definition in the first place? Is this just a "quirk" of going from integral over interval to finite sum? For big \$N\$ it converges to same result, for small \$N\$ window filters probably should not have any zeros to preserve all raw data points. Padding any amount of zeros (or none at all) is up to practical consideration.
--- Quote from: rhb on December 10, 2022, 04:02:35 pm ---Compute the FFT of the "flat topped" triangle of length N, insert N zeros in the center and back transform. Now you can see the half sample shift.
--- End quote ---
Not sure I completely understand idea behind this. Can you show an example calculation?
--- Quote from: rhb on December 10, 2022, 05:49:56 pm ---A few caveats about the DFT. There are 3 scaling factors used, 1/N,N & 1/sqrt(N). Moving data between software packages using different scaling can be a problem if you're not paying attention. There are also 2 sign conventions for the exponent and 2 intervals commonly used, -Pi to Pi and 0 to 2*Pi. The former is the mathematical definition, but the implementations are typically 0 to 2*Pi. F(i*Pi-x) = F((i+2)*Pi-x) which is why the "negative" frequencies are at the end of the transform in descending order. So there are 12 possible implementations of the same concept. This is a common source of confusion and error.
--- End quote ---
Good point, I crosscheck implementations to be sure. FFT(*) functions have implicit assumptions that should be checked before use.
Just for reference, fft() function help in R states
--- Quote ---The FFT is fastest when the length of the series being transformed is highly composite (i.e., has many factors). If this is not the case, the transform may take a long time to compute and will use a large amount of memory.
--- End quote ---
Not sure which implementation it uses internally. Simple test (replicated many times to get stable result) with 1024 (2^10), 1155 (3*5*7*11), 1021 (prime) shows ~15x longer calculation for prime and ~1.7x longer for composite vs power of 2.
--- Code: --- relative duration
fft(seq_len(1024)) 1.000
fft(seq_len(1155)) 1.683
fft(seq_len(1021)) 14.663
--- End code ---
Navigation
[0] Message Index
[#] Next page
[*] Previous page
Go to full version