On Sep 27, 1:31�am, Raeldor <rael...@gmail.com> wrote:> > So, I tried Robert's suggestions,"guaranteed or your money back..."> but the output still looks a little > inconsistent (and I think incorrect).you might apply for a refund.> �I think it may be because the > FFT library I am using may expect the data in a slightly different > format. �To do a basic FFT on a set of real data 0...n, the library > requires the data to be transformed into an 'even-odd' array. �This is > done by a function vDSP_ctoz, which takes an array of real numbers and > transforms it into a COMPLEX_SPLIT type, which basically contains 2 > arrays realp[] and imagp[]. �The result of the transformation call > puts the even values into the realp[] array, and the odd values into > the imagp[] array. �I don't think these are actually complex numbers, > because the FFT call is a real-valued FFT function, I think it's just > using these arrays to order the data so that the even numbers come > first followed by the odd. �What isn't clear to me is if I still need > to do this after I have created the mirrored data as Robert suggested. >there *are* fast DCT routines that do not require this mirroring and do not use the DFT as an intermediate step. i was trying only to relate the DCT to the DFT which is what your subject header was about. the DFT of a purely real input will have this "hermitian" complex- conjugate symmetry (real part is even symmetry, imag part is odd). but if the real input has even symmetry, that is x[-n] = x[n] which is the same as x[N-n] = x[n] then the DFT is not only hermitian symmetric, the imaginary part is zero so it's real and symmetric. if you start out with N/2 points, you can reflect them, but for the DCT-II, the reflection is symmetric about the n=-1/2 axis, not about the n=0 axis. but if you fix the effect of that in the frequency domain (by multiplying X[k] by e^(j*pi*k/N)) that should make it real in the frequency domain.> Also, another thing is that Robert's explanation of having to mirror > the data tied in with what I was reading in the FFTW documentation, > but the Wikipedia article for DCT-II says that the data must also be > interlaced with zero values, I quote... > > This transform is exactly equivalent (up to an overall scale factor of > 2) to a DFT of 4N real inputs of even symmetry where the even-indexed > elements are zero. That is, it is half of the DFT of the 4N inputs yn, > where y2n = 0, y2n + 1 = xn for , and y4N - n = yn for 0 < n < 2N. >what's the difference between "identical" and "equivalent"? if you do it that way (so you are starting out with N/2 points, reflecting it (getting you N points), and then zero inserting (getting you 2N points, my N is twice that of Wiki's N in the DCT article), then you are able to perfectly center this about the n=0 axis and the adjustment of e^(j*pi*k/N)) is not necessary.> Which wasn't mentioned in Robert's explanation. �So, I'm not sure > whether I should be passing in 2N values or 4N values.disregard anything specific that i'm saying. pay close attention to the Wikipedia definition and try to implement that. as i've said, i've never actually used a DCT for anything. r b-j
FFT vs DCT
Started by ●September 26, 2011
Reply by ●September 27, 20112011-09-27
Reply by ●September 27, 20112011-09-27
On 9/27/2011 1:31 AM, Raeldor wrote:> On Sep 26, 6:10 pm, Jerry Avins<j...@ieee.org> wrote: >> On 9/26/2011 6:24 PM, Raeldor wrote: >> >> ... >> >> >> >> >> >> >> >> >> >>> So, it might seem to YOU that it's a question that can be easily >>> answered by looking in somewhere like Wikipedia, but I can tell you >>> that I've read the articles there several times, and it's not easily >>> answered when you don't have a math or DSP heavy background. A lot of >>> these articles seem to say one thing, and then contradict it again a >>> few sentences later. Either that, or they use terminology which, >>> while very essential, can often be overlooked or misinterpreted by >>> someone from a different background (for me, I am a programmer). I'm >>> not asking others to do my work, believe me, I've read this stuff, >>> implemented, tested, experimented, and have been for several months, >>> but sometimes I need to go back and confirm that even my BASIC >>> understanding of some of this stuff has not sent me on the wrong >>> track. I do not ask these questions intending to be disrespectful or >>> lazy, I ask them because I hope to get a more helpful laymens-terms >>> answer from people who frequent this board because they may actually >>> care to help others understand. >> >> When you dig into an article that confuses you, cite the confusing parts >> instead of asking what could seem like a dumb or lazy question. What you >> wanted to know is actually fairly deep. It is about using the algorithm >> of an FFT to actually perform a different transform by some pre- and >> post processing. That doesn't makes the transforms the same. >> >> Jerry >> -- >> Engineering is the art of making what you want from things you can get. > > Thanks for all the positive replies. You're right, I probably should > have been more specific, but I was at a point where I was starting to > even doubt my limited understanding. > > So, I tried Robert's suggestions, but the output still looks a little > inconsistent (and I think incorrect). I think it may be because the > FFT library I am using may expect the data in a slightly different > format. To do a basic FFT on a set of real data 0...n, the library > requires the data to be transformed into an 'even-odd' array. This is > done by a function vDSP_ctoz, which takes an array of real numbers and > transforms it into a COMPLEX_SPLIT type, which basically contains 2 > arrays realp[] and imagp[]. The result of the transformation call > puts the even values into the realp[] array, and the odd values into > the imagp[] array. I don't think these are actually complex numbers, > because the FFT call is a real-valued FFT function, I think it's just > using these arrays to order the data so that the even numbers come > first followed by the odd. What isn't clear to me is if I still need > to do this after I have created the mirrored data as Robert suggested.This is another example of using the complex FFT algorithm for another purpose -- here, an efficient real-only FFT -- with some pre- and post processing. There are other examples. Just as DCT and DST can be computed with purpose-built algorithms, so can a real FFT, if one cares.> Also, another thing is that Robert's explanation of having to mirror > the data tied in with what I was reading in the FFTW documentation, > but the Wikipedia article for DCT-II says that the data must also be > interlaced with zero values, I quote... > > This transform is exactly equivalent (up to an overall scale factor of > 2) to a DFT of 4N real inputs of even symmetry where the even-indexed > elements are zero. That is, it is half of the DFT of the 4N inputs yn, > where y2n = 0, y2n + 1 = xn for , and y4N - n = yn for 0< n< 2N. > > Which wasn't mentioned in Robert's explanation. So, I'm not sure > whether I should be passing in 2N values or 4N values.Go with the documentation, and remember that Robert warned you that he was a bit fuzzy. It's only 4N values if you count the obligatory zeros as data rather than as spacers. Jerry -- Engineering is the art of making what you want from things you can get. -- Engineering is the art of making what you want from things you can get.
Reply by ●September 27, 20112011-09-27
On 9/27/11 11:35 AM, robert bristow-johnson wrote:> On Sep 27, 1:31 am, Raeldor<rael...@gmail.com> wrote: >> >> So, I tried Robert's suggestions, > > "guaranteed or your money back..." > >> but the output still looks a little >> inconsistent (and I think incorrect). > > you might apply for a refund. >okay, just for my curiousity, i decided to look at this closely. i will not bother to type in every intermediate equation, since it's basically an exercise of mathematical masturbation. i made a couple of mistakes, but am not certain why. referring to the DCT-II definition (which Wikipedia says is the most common definition of the DCT), i am changing the N in Wikipedia's definition to N/2 (to emphasize that N must be even and because i want the DFT length to be N) and i am changing their x[n] and X[k] to y[n] and Y[k] because i want the x[] and X[] to be the DFT in and out. DCT-II ("type 2 DCT"): N/2-1 Y[k] = SUM{ y[n] * cos(2*pi*(n+1/2)*k/N) } for 0 <= k < N/2 n=0 so, if you mathematically jerk off a bit (replace cos() with complex exp, turn one of the summations around, etc.) you get: Y[k] = X[k] * e^(-j*pi*k/N) for 0 <= k < N/2 where N-1 X[k] = SUM{ x[n] * e^(-j*2*pi*n*k/N) } n=0 and { (1/2) * y[n] 0 <= n < N/2 x[n] = { { (1/2) * y[N-1-n] N/2 <= n < N so, instead of keeping the evens, we keep the first half. and instead of multiplying the output of the DFT by e^(+j*pi*k/N) we multiply by e^(-j*pi*k/N) . even though i have never used the DCT for anything in any project i have ever done, i can guarantee that if the Wikipedia article is accurate about the DCT-II, the above math is also correct. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●October 2, 20112011-10-02
Raeldor <raeldor@gmail.com> writes:> 2. I keep reading about 'real-odd' and 'real-even'. What does this > mean?"Real-odd" and "real-even" is the terminology used in the FFTW documentation. FFTW takes the position that, among the zoo of Fourier-related transforms, the complex transform is the only fundamental one and all the other ones are just special cases of the (complex) FFT. The fact that some authors have chosen to call "DCT-I" a certain special case of FFT is a historical accident that caused infinite confusion. The whole DCT/DST terminology is better avoided altogether. Starting from the Fourier transform of a complex signal, one can specialize the FFT to the case where the input is real (as opposed to complex). This is an important case. (However, nobody seems to care about the case of imaginary input, which is equally interesting.) Another specialization of the FFT is for the case where the input has certain symmetries. For example, the input can be an "even" function, that is f(-x) = f(x). Or maybe the input is symmetric around -1/2 instead of around 0. Or maybe the symmetry is "odd", that is, f(x) = -f(x). Beware that "even/odd" are trickier for complex functions. The proper definition of "even" for complex functions should be f(x) = conj(f(-x)). (Somehow nobody seems to care about the case of symmetric imaginary or complex functions, which seems equally important.) FFTW implements some of these special cases, as described in manual. The meaning of these special cases can always be inferred from the meaning of the corresponding nonspecialized Fourier transform. In particular, there are no surreptitious factors of 2 or sqrt(2) or whatever to remember. What you get out of a real-even transform is the same as what you would get if you make the input symmetric and call the complex FFT routine. What you get may also be the same of what somebody else calls a DCT-I (or it may not be), but in any case it is a waste of time to try to understand the DCT as an independent concept. FFTW uses this ontology also as an implementation principle. Specifically, FFTW comes with a special-purpose compiler that generates C code for complex FFTs, real FFTs, even/odd/whetever FFTs. This compiler does not know anything about the special cases. It only knows about the complex FFT and a few algebraic rules of the form 0+x=x and similar high-school algebra. To generate C code for a real-even FFT, it generates an internal representation of the corresponding complex FFT, it notes that certain inputs are the same, it eliminates common subexpressions, it applies the high-school rules, and it outputs the result. In all cases it produces an algorithm that is as good as (or better then) the published algorithms for the special-case FFTs. I have yet to find a single "DCT" algorithm that is not a trivial specialization of a complex FFT algorithm (trivial in the sense that our compiler can do it with high-school algebra). So the real algorithmic depth of the matter lies in the various complex FFT algorithms (Cooley-Tukey, PFA, Rader, etc.) and the specializations add no further insight. I think the whole DCT/DST mess makes much more sense if you look at it in the way I just described. Regards, Matteo Frigo
Reply by ●October 2, 20112011-10-02
Matteo Frigo <athena@fftw.org> wrote:> Raeldor <raeldor@gmail.com> writes:>> 2. I keep reading about 'real-odd' and 'real-even'. What does this >> mean?> "Real-odd" and "real-even" is the terminology used in the FFTW > documentation.> FFTW takes the position that, among the zoo of Fourier-related > transforms, the complex transform is the only fundamental one and all > the other ones are just special cases of the (complex) FFT. The fact > that some authors have chosen to call "DCT-I" a certain special case of > FFT is a historical accident that caused infinite confusion. The whole > DCT/DST terminology is better avoided altogether.Well, it usually goes back to boundary conditions, and that usually goes back to the physics. In actual problems, it is usually the boundary conditions that come first. For analytical solutions, it isn't so hard to just through out the solutions that don't fit, but for numerical solutions you don't want to waste time solving for terms where the answer is known to be zero.> Starting from the Fourier transform of a complex signal, one can > specialize the FFT to the case where the input is real (as opposed to > complex). This is an important case. (However, nobody seems to care > about the case of imaginary input, which is equally interesting.)The Fourier series has periodic boundary conditions, useful for periodic functions. Other common boundary conditions are that the function or its derivative goes to zero at the boundary. DST takes care of the former, DCT the latter. You might also have a case where the function is zero at one end, the derivative at the other end. (This commonly comes out for the acoustic modes of a tube open at one end, closed at the other end.) Even more generally, the boundary condition can be some combination of the function and derivative, though that is somewhat rare.> Another specialization of the FFT is for the case where the input has > certain symmetries. For example, the input can be an "even" function, > that is f(-x) = f(x). Or maybe the input is symmetric around -1/2 > instead of around 0. Or maybe the symmetry is "odd", that is, f(x) = > -f(x). Beware that "even/odd" are trickier for complex functions. The > proper definition of "even" for complex functions should be f(x) = > conj(f(-x)). (Somehow nobody seems to care about the case of symmetric > imaginary or complex functions, which seems equally important.)If you use 0 for one of the end points, then even/odd seems a little strange. If you use a region symmetric around zero, then even or odd makes some sense, which gets you cosines and sines, respectively. But the basis functions for both DST and DCT are not all even or all odd around their midpoint.> FFTW implements some of these special cases, as described in manual. > The meaning of these special cases can always be inferred from the > meaning of the corresponding nonspecialized Fourier transform. In > particular, there are no surreptitious factors of 2 or sqrt(2) or > whatever to remember. What you get out of a real-even transform is the > same as what you would get if you make the input symmetric and call the > complex FFT routine. What you get may also be the same of what somebody > else calls a DCT-I (or it may not be), but in any case it is a waste of > time to try to understand the DCT as an independent concept.As for DCT-I, there comes the question of where the boundary, and therefor boundary condition, relative to the sample points. For the DCT-I through IV cases, the boundary is either on, or half way between, two sample points. In the general case, it could be anywhere between two points.> FFTW uses this ontology also as an implementation principle. > Specifically, FFTW comes with a special-purpose compiler that generates > C code for complex FFTs, real FFTs, even/odd/whetever FFTs. This > compiler does not know anything about the special cases. It only knows > about the complex FFT and a few algebraic rules of the form 0+x=x and > similar high-school algebra.Can it do mixed boundary conditions? Or where the boundary is other than on or half way between sample points?> To generate C code for a real-even FFT, it > generates an internal representation of the corresponding complex FFT, > it notes that certain inputs are the same, it eliminates common > subexpressions, it applies the high-school rules, and it outputs the > result. In all cases it produces an algorithm that is as good as (or > better then) the published algorithms for the special-case FFTs.That is the way it is usually taught to physics students. Consider the usual violin string problem. Assume that the solution is sums of sines and cosines. The string is fixed at one end, f(0)=0, so through out the cosine terms. The string is fixed at the other end, f(l)=0. sin(al)=0, al=2npi, a=2npi/l.> I have > yet to find a single "DCT" algorithm that is not a trivial > specialization of a complex FFT algorithm (trivial in the sense that our > compiler can do it with high-school algebra). So the real algorithmic > depth of the matter lies in the various complex FFT algorithms > (Cooley-Tukey, PFA, Rader, etc.) and the specializations add no further > insight.-- glen
Reply by ●October 3, 20112011-10-03
glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:> Well, it usually goes back to boundary conditions, and that usually > goes back to the physics. In actual problems, it is usually the > boundary conditions that come first. For analytical solutions, it > isn't so hard to just through out the solutions that don't fit, but > for numerical solutions you don't want to waste time solving for > terms where the answer is known to be zero.Correct. The special cases (real-FFT, real-{even,odd} FFT, etc.) are obviously important in practice. As a proof that I really believe this, we spent time implementing all of them in FFTW (and debugging the code). My point is that they don't have an independent *definition*, and that trying to view them as separate transforms only leads to unnecessary confusion.> The Fourier series has periodic boundary conditions, useful for > periodic functions. Other common boundary conditions are that the > function or its derivative goes to zero at the boundary. DST takes > care of the former, DCT the latter.This is incorrect. For the "DST" to be applicable the function must be zero at the boundary, as you say, and in addition the function must be odd around the boundary. If the function is not odd then you cannot apply the "DST". So the odd-ness of the function is the important property, and the fact that it is zero at the boundary is a (trivial) corollary of its odd-ness.> But the basis functions for both DST and DCT are not all even or > all odd around their midpoint.This is not what I am saying. I am saying this: consider e.g. a function f defined on [-n..n]. Let f be real and even around 0. Then its Fourier transform F, because 2 cos(x) = exp(ix)+exp(-ix), can be written in terms of cos() only, without using sin() and i=sqrt(-1). Moreover F is also real and even. So far everything is simple. Now comes the confusing part. Because half the input and half the output are redundant, people talk about a transform called the "DCT-I" that maps f([0..n]) into F([0..n]). If you look at [0..n] only, of course neither the function nor the cosine basis are symmetric. However, by doing this step you have just introduced a lot of unnecessary complexity. That is, you now have to remember whether the period is n, n+1, n-1, or 2n-1, or whatever (the period happens to be 2n+1 in this case, and the length of the DCT-I is n+1). You have to remember that some terms have a factor of 2 in front of them (coming from 2cos(x) = exp(ix)+exp(-ix)) and some don't (for x=0). You have a PI instead of 2*PI at the denominator, which makes people ask the question "why isn't this DCT-I the real part of a FFT". You have to remember whether the denominator has n+1, n-1, n, or whatever. None of this is necessary if you just remember that this is the standard FFT of a real-even function on [-n..n]. Your very post shows an instance of this confusion, because you were under the impression that f'=0 at the boundary is a sufficient condition for the DCT-I (it is not; the function must be even or else the sine components do not cancel out). All this confusion can be avoided by viewing these transforms for what they are, i.e., special cases of the FFT. This view avoids the n/n+1/2n mess and the spurious factors of 2. Regards, Matteo Frigo
Reply by ●October 3, 20112011-10-03
On 10/3/11 9:31 AM, Matteo Frigo wrote:> > My point is that they don't have an independent *definition*, and that > trying to view them as separate transforms only leads to unnecessary > confusion.A-friggin'-men! i have felt that way since the very, very first time i had ever read anything regarding the DCT. it's just another (among many) app of the DFT. Matteo, please give my best to Steve. i ain't living in the Boston area at the moment (renting the condo in Waltham), but i had always wanted to pay FFTW headquarters a visit, when i was down there. (also want to visit the office of Alan O sometime, and hadn't done that either.) -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●October 3, 20112011-10-03
Matteo Frigo <athena@fftw.org> wrote: (snip)> Your very post shows an instance of this confusion, because you were > under the impression that f'=0 at the boundary is a sufficient condition > for the DCT-I (it is not; the function must be even or else the sine > components do not cancel out).The derivative of sine at zero is never zero, so that is a sufficient condition to remove the sine terms. Or, more generally, an odd function must be zero at zero, and an even function must have a derivative zero at zero. I suppose in the general data processing sense it makes some sense to do a generalized transform routine. In the DSP case, where everything is very specific, such as the N=8 DCT commonly used in image transforms, there is no need for the generalization. -- glen
Reply by ●October 3, 20112011-10-03
glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:> The derivative of sine at zero is never zero, so that is a > sufficient condition to remove the sine terms.The premise is true. The conclusion is false. A bunch of nonzero sine coefficients can add up to zero in multiple ways---the coefficients don't have to be all zero. E.g., the DFT of f(x)=x^3 (an odd function such that f'(0)=0) has no cosine terms (since the function is odd), and therefore some sine term must be nonzero (or else the DFT of f is identically zero, which is absurd).> Or, more generally, an odd function must be zero at zero, > and an even function must have a derivative zero at zero.Yes, this is what precisely I am saying: "f is even" implies "f'(0)=0", but the converse is not true. To cancel the sine terms you need "f is even". "f'(0)=0" is not sufficient. That is, the even/odd symmetry is the important thing, and the boundary condition by itself is irrelevant. Saying that "the DCT can be used with any function f such that f'(0)=0" is wrong, as the f(x)=x^3 example demonstrates. The DCT can only be used for even functions (which then necessarily have f'(0)=0). This is why I am saying that talking about the DCT in the first place is a bad idea---it just confuses the essence of a pretty simple matter.
Reply by ●October 3, 20112011-10-03
Matteo Frigo <athena@fftw.org> wrote:>FFTW takes the position that, among the zoo of Fourier-related >transforms, the complex transform is the only fundamental one and all >the other ones are just special cases of the (complex) FFT.Hopefully there's a word "discrete" missing from the above. "Among the zoo of Fourier-related distrete transforms..." Steve
