Computing a DCT using an FFT, Signal analysis 
Computing a DCT using an FFT, Signal analysis 
Dec 9 2005, 06:40
Post
#1


Group: Members Posts: 35 Joined: 15April 04 From: Australia Member No.: 13509 
I am trying to make sense of the output of the FFT when used to calculate a type2 DCT (DCTII). I managed to scrape together enough documentation that infomed me that, if I wanted to perform a DCTII using an FFT algorithm, then the data should be arranged to be real even symmetric, with every even indexed data element set to zero.
No problem there. In addition, the outputs are correct, but half the expected value (I can see the output needs scaling by 2 to compensate for the fact that every second sample of the input is zero). However, I am tyring to figure out how the output comes to be what it is? In particular, where the zero comes from (term 4 & term 12) [counting from 0]. Input Real : 0 x0 0 x1 0 x2 0 x3 0 x3 0 x2 0 x1 0 x0 Imag : 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Output Real : ya yb yc yd 0 yd yc yb ya yb yc yd 0 yd yc yb Imag : 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [grumble: looks good in fixed width font, but spaces are not turning out here] In my real application, I am using a 16 point DCT, which equates to a 64 point FFT, but the above is representative. Any hints? PS: I know there are better ways of doing DCTs, but in my implementation, I only have access to an FFT routine. Cheers,Owen. 


Dec 9 2005, 12:57
Post
#2


Mad Scientist Group: Admin Posts: 4890 Joined: 24September 01 Member No.: 13 
Use code or codebox to get it looking right.
Do you mean you cannot do a rotation or similar before the FFT? This method looks very wasteful, there are much better ways of doing it via an FFT but they need a some pre and postprocessing. 


Dec 9 2005, 21:32
Post
#3


Group: Members Posts: 35 Joined: 15April 04 From: Australia Member No.: 13509 
QUOTE (Garf @ Dec 9 2005, 03:57 AM) Use code or codebox to get it looking right. Thanks Garf, I am not so worried about appearances, more wanting to understand why it is so. I am not that familiar with the mathematics, only how to apply the end results. What I expected, and I should have put in my original post, is... I expected... Real : ya yb yc yd yd yc yb ya ya yb yc yd yd yc yb ya Imag : (all zero) That is, no zero terms in the real output. QUOTE (Garf @ Dec 9 2005, 03:57 AM) Do you mean you cannot do a rotation or similar before the FFT? This method looks very wasteful, there are much better ways of doing it via an FFT but they need a some pre and postprocessing. Yikes!! Not familiar with 'rotations' (but I have come across the term a lot). My application is targetted to hardware (so I can't use FP), and I am doing integer maths. I am trying to get away with 16 bits of precision (very convenient for the hardware). I agree, the method is very wasteful, and I intend to use the radix4 method of preprocessing to reduce the FFT size to n (rather than 4n). At the moment, I am trying out different things  but the output did give me a surprise. Nonmathematical documentation on doing a DCT using an FFT is sparse. I only found info on how the input is arranged [first post], but nothing on how the output is arranged. I did try a 16bit implementation of BinDCT, (internally, the precision extended to 22 bits) but the inaccuracies killed my application. I guess accuracy was the tradeoff for blistering speed! I suppose my essential question is : is the output correct, or have I stuffed up? As mentioned, documentation was sparse, so I could be loading the FFT indicies wrong? Thanks for the reply Garf. Owen. 


Dec 12 2005, 01:13
Post
#4


Group: Developer Posts: 1245 Joined: 16December 02 From: Australia Member No.: 4097 
I agree that there is often not enough documentation on how to do a simple DCT using an FFT. Though it is not the best, computationalwise, it is however interesting to see the relationship between the two, esp. why the DCT is used more often than the FFT in compression (the evensymmetry creates smoothness that results in less high frequency energy).
Anyway, I was interested in this problem a while ago as well and found this link to be the best I've found. Not sure if this is any help or not, but here is some MATLAB code I wrote during that time, which calculates the DCT using FFT. CODE clear all; close all; x=[1 2 3 4]; N=length(x); d = dct(x); x2=[4 3 2 1 1 2 3 4]; f=fft(x2); n = 0:2*N1; shift=exp(j*2*pi*(N0.5)*n/(2*N)); f2 = real(f./shift); d2 = f2(1:N)/sqrt(2*N); d2(1)=d2(1)/sqrt(2); disp(d); disp(d2); Notice how I don't do any zerosetting (no idea why you need to anyway). I just take my original signal, mirror it and then take the fft. Now the thing to note is that the Fourier transform of an even symmetry signal is real but we can't make true even symmetry signals on a computer, so we have to compensate for the phase shift caused by the timeshift of 3.5 samples, hence the division by the FFT of the variable 'shift'. Then I scale the coefficients appropriately (to make them orthonormal), and voila. Hope that helps in a little way EDIT: fixed a spelling mistake This post has been edited by QuantumKnot: Dec 12 2005, 04:55 


Dec 13 2005, 23:25
Post
#5


Group: Members Posts: 35 Joined: 15April 04 From: Australia Member No.: 13509 
Hi QuantumKnot,
Thanks for your reply. Indeed I had read the same page on using an FFT to do a DCT, and understood the reasoning behind having to mirror the coefficents before placing them into a FFT. I did a bit of experimenting, and created a "brute force" DCT using floating point maths on a PC. I had some code for a DCTII, and rewrote it to perform a DCTI. The results were quite different. I also examined the documentation on FFTW, in particular about the different types of DCT. They refer to them as 00, 01, 10, and 11. The 0 stands for no shift, the 1 stands for a halfsample shift. The first figure refers to the input, the last to the output. The FFTW routines REDFT00 = DCTI, REDFT10 = DCTII, REDFT01 = DCTIII, REDFT11 = DCTIV. Apparently, filling an FFT without padding but with mirroring (as you have done) is equivilant to a DCTI. Filling with padding and mirroring (as described in my first post) is supposed to be equivilant to a DCTII (which is what I need). However, this takes the FFT to 4n. Further documentation (wikipedia in fact) hints at this "For example, a typeII DCT is equivalent to a DFT of size 4n with realeven symmetry whose evenindexed elements are zero. One of the most common methods for computing this via an FFT (e.g. the method used in FFTPACK and FFTW)..." It goes on a bit further to say computations can be shorthened to O(nlogn) + O(n) if a preprocessing step of radix4 is used (O(n) steps) and the FFT is then reduced to size n. A decent saving. As I said before, my application is power impoverished hardware, where all I have is a 16bit integer FFT routine (which works surprisingly well I must say). Additionally BinDCT may be very fast, but it is also considerably less accurate. Again, thanks for your reply. A cursory look at your sample code has been enlightening, and I am about to print it and examine it in detail. Cheers and thanks, Owen. 


Dec 14 2005, 00:08
Post
#6


Group: Developer Posts: 1245 Joined: 16December 02 From: Australia Member No.: 4097 
QUOTE (Mustardman @ Dec 14 2005, 08:25 AM) Apparently, filling an FFT without padding but with mirroring (as you have done) is equivilant to a DCTI. That is weird since I've also compared the results of the above MATLAB code with FFTW's REDFT10 r2r mode and got the same results (with the appropriate scaling). If you check FFTW's user manual (p. 11), it says that a DCTII is symmetric and even around j = 0.5, which is exactly how I'm mirroring it above. That is, a sequence of abcd is mirrored as dcbaabcd, where if we align the right 'a' as 0, then the axis of symmetry (represented as a '') is at 0.5. The DCTI is symmetric around j = 0, so the sequence abcd is mirrored as bcdabcd, where the axis of symmetry is on the 'a'. This post has been edited by QuantumKnot: Dec 14 2005, 00:10 


Dec 15 2005, 00:51
Post
#7


Group: Members Posts: 35 Joined: 15April 04 From: Australia Member No.: 13509 
QUOTE (QuantumKnot @ Dec 13 2005, 03:08 PM) That is weird since I've also compared the results of the above MATLAB code with FFTW's REDFT10 r2r mode and got the same results (with the appropriate scaling). If you check FFTW's user manual (p. 11), it says that a DCTII is symmetric and even around j = 0.5, which is exactly how I'm mirroring it above. That is, a sequence of abcd is mirrored as dcbaabcd, where if we align the right 'a' as 0, then the axis of symmetry (represented as a '') is at 0.5. The DCTI is symmetric around j = 0, so the sequence abcd is mirrored as bcdabcd, where the axis of symmetry is on the 'a'. That may explain where I am going astray. I was unsure how to perform the mirroring, so I just filled the FFT entries (excluding the zero padding, another issue) in the order abcddcba. An assumption I made (quite reasonable I think) is that since a DFT is periodic, this would affect either [1] the order of the output (when reading the same indexes), or [2] which outputs (indexes) of the FFT are relevant. I did get my 4n size integer FFT working without any grief, and the results are (within 0.05% error) consistant with the original floatingpoint DCT. I have yet to code the radix4 step, but from my reading so far, that seems fairly straight forward (just I wait!!). It could well be that I can fill the entries properly (thanks very much for the explanation you have given above!!) and do a 2n FFT. I will give that a try... Thanks, Owen. edit: OK, an update. I tried inserting values into the FFT (size 2n) as suggested. I was unhappy to find out that the outputs were nowhere near what was expected. Even worse, and more telling, the imaginary components of the output were not zero (as they should be if we shove real values in and want a DCT). Do you know if the "2n" method requires magnitude calculation of the real & imag outputs? (you don't mention anything about it in your post... I would imagine though that this would [should?] not be the case?). Delving into the FFTW manual........ Cheers, Owen. This post has been edited by Mustardman: Dec 15 2005, 04:08 


Jan 13 2006, 06:49
Post
#8


Group: Members Posts: 1593 Joined: 24March 02 From: Revere, MA Member No.: 1607 
http://www.ee.columbia.edu/~marios/courses...%20(Thesis).pdf
Here is another paper for reference describing faster way's to compute a DCT using an FFT.  College student/IT Assistant



Jan 13 2006, 10:03
Post
#9


Mad Scientist Group: Admin Posts: 4890 Joined: 24September 01 Member No.: 13 
QUOTE (HotshotGG @ Jan 13 2006, 07:49 AM) http://www.ee.columbia.edu/~marios/courses...%20(Thesis).pdf Here is another paper for reference describing faster way's to compute a DCT using an FFT. The paper is about the MDCT, not the DCT. And the MDCT is based on a type 4 DCT, while he needs a type 2 DCT. This post has been edited by Garf: Jan 13 2006, 10:04 


Jan 16 2006, 06:59
Post
#10


Group: Developer Posts: 1245 Joined: 16December 02 From: Australia Member No.: 4097 
QUOTE (Mustardman @ Dec 15 2005, 09:51 AM) edit: OK, an update. I tried inserting values into the FFT (size 2n) as suggested. I was unhappy to find out that the outputs were nowhere near what was expected. Even worse, and more telling, the imaginary components of the output were not zero (as they should be if we shove real values in and want a DCT). Do you know if the "2n" method requires magnitude calculation of the real & imag outputs? (you don't mention anything about it in your post... I would imagine though that this would [should?] not be the case?). Delving into the FFTW manual........ Cheers, Owen. Looks like I forgot to make a reply to this No, there is no need to do a magnitude calculation on the output. If you mirror the sequence (size 2n) and do an FFT, you will NOT get pure real values, since there is a phase shift. You need to correct the phase by multiplying the output with exp(j*2*pi*(N0.5)*n/(2*N)), and you will find that the output is now real (no imaginary component). 


Mar 30 2006, 21:04
Post
#11


Group: Members Posts: 1 Joined: 30March 06 Member No.: 29012 
% use FFT to commute DCT and then invert DCT using inverse FFT
x=1:6; disp(x); N=length(x); d = dct (x); x2=[fliplr(x) x]; f=fft(x2); n = 0:2*N1; shift=exp(j*2*pi*(N0.5)*n/(2*N)); f2 = real(f./shift); d2 = f2(1:N)/sqrt(2*N); d2(1)=d2(1)/sqrt(2); disp(d) disp(d2) % now invert d2(1)=d2(1)*sqrt(2); d2=d2*sqrt(2*N); f=[d2 0 fliplr(d2(2:N))].*shift; x2=ifft(f); x=real(fliplr(x2(1:N))); disp(x); 


Jun 20 2012, 15:31
Post
#12


Group: Members Posts: 1 Joined: 20June 12 Member No.: 100848 
I think the solution which is mentioned in below address in more optimized than your solution because it uses Npoint FFT to calculate Npoint DCT , but you use 2*N point FFT: http://fourier.eng.hmc.edu/e161/lectures/dct/node2.html I have checked it in MATLAB it was okey except it has forgotten to use sqrt(2) and sqrt(2*N) coeficients as you use them. clc; clear all; close all; x = [ 1, 2, 3, 4, 3, 2, 1, 5, 1, 2, 3, 4, 3, 2, 1, 5, 1, 2, 3, 4, 3, 2, 1, 5, 1, 2, 3, 4, 3, 2, 1, 5 ]; matlab_idct = idct(x); N=length(x); n = 0:N1; shift = exp( (n*pi*j) / (2*N) ); shift = shift * sqrt(2*N); shift(1) = shift(1) / sqrt(2); y = shift .* x ; y = real(ifft(y)); for i=0:((N/2)1) my_idct(2*i + 1) = y(i + 1) ; my_idct(2*i + 2) = y(N  i) ; end 


LoFi Version  Time is now: 4th July 2015  10:20 