# 1D as a 3D FFT - what’s wrong here?

I implemented the algorithm to calculate 2D DFT derived from 1D DFT. It works great, and makes my calculations much more efficiency then regular 1D DFT.

But now I want to make 3D DFT derived from 1D DFT and it doesn’t work for me. For 3 days I’ve tried to solve it, but I can’t, so I would like to ask you for help.

I make all math transforms by following the explanations in found in book “Theory and application of digital signal processing” Lawrence R. Rabiner, Bernard Gold - chapter 6.8 titled “A Unified Approach to the FFT“.

According to book example, let’s say we have $\ N=60$ point signal. Let me remind definition of the DFT, which is:

$\ X ( k ) = \sum _ { n = 0} ^ { N - 1} x ( n ) W ^ { n k }$

$\ N=60$ points we can expres as a $60= 5\times 12$. So $\ n$ and $\ k$ can be written as:

$\ n = M l + m$

$\ k = L r + s$

Where:

$\ M = 12$

$\ L = 5$

$\ l = 0,1,..,L - 1$

$\ m = 0,1,..,M - 1$

$\ s = 0,1,..,L - 1$

$\ r = 0,1,..,M - 1$

So now we can write DFT as:

$\ X ( k ) = X ( L r + s ) = \sum _ { m = 0} ^ { M - 1} \sum _ { l = 0} ^ { L - 1} x ( M l + m ) W ^ { ( M l + m ) ( L r + s ) }$

And after some tweaking with power of W, we find that:

$\ W ^ { ( M l L r ) } = 1$

And we can write DFT as:

$\ X ( L r + s ) = \sum _ { m = 0} ^ { M - 1} W ^ {m (L r + s) } \sum _ { l = 0} ^ { L - 1} x ( M l + m ) W ^ { M s l }$

So great. Now I have 2D DFT, and it works great for me.

THE PROBLEMS STARTS HERE

Following my example with $\ N=60=5x12$, we can see that 12 can be expressed as $\ 12=3x4$.

So $\ m$ and $\ r$, which have values $\ 0, 1, …, M-1$, we can express as:

$\ m = Z c + z$

$\ r = C v + y$

Where:

$\ Z = 4$

$\ C = 3$

$\ c = 0, 1, .., C - 1$

$\ z = 0, 1, .., Z - 1$

$\ y = 0, 1 , .., C - 1$

$\ v = 0, 1, .., Z - 1$

So now the $\ n$ and $\ k$ are:

$\ n = M l + Z c + z$

$\ k = L (C v + y) + s$

And then the DFT is:

$\ X ( k ) = X ( L (C v + y) + s ) = \sum _ { z = 0} ^ { Z - 1} \sum _ { c = 0} ^ { C - 1} \sum _ { l = 0} ^ { L - 1} x ( M l + Z c + z ) W ^ { ( M l + Z c + z ) ( L (C v + y) + s ) }$

And now after some tweaking with power of W, we find that:

$\ W ^ { ( M l L C v ) } = 1$

$\ W ^ { ( M l L y ) } = 1$

$\ W ^ { ( Z c L C v ) } = 1$

And finally we can write DFT as:

$\ X ( L (C v + y) + s ) = \sum _ { z = 0} ^ { Z - 1} W ^ {z (L C v + L y + s) } \sum _ { c = 0} ^ { C - 1} W ^ {c (Z L y + Z s) } \sum _ { l = 0} ^ { L - 1}x ( M l + Z c + z ) W ^ { M s l }$

And it doesn’t work for me, I have some crazy output values. And Have no idea if my math transforms (which I described above) are wrong, or maybe I just made wrong implementation. Of course I am not asking you about my implementation (which I haven’t shown here). First I want to be sure that my math transforms are OK?

For any help great thanks in advance.

• Does your 12 point work on its own?
– user28715
Apr 7, 2018 at 15:01
• Hello, Great thanks for any reply, Stanley, you're my good spirit :) But to the point: I am not sure what you are asking. If you mean 12 point inside that 3D N=5x3x4 DFT, then I am not sure how to check it. But if you are asking about 2D DFT, so yes it works for input N=12, after I make matrix 3x4. Also works for N=60=5x12, and for any other 2D matrix. But I wrote it clearly in my question, just before the line "THE PROBLEMS STARTS HERE". So I don't think you are asking about 2D DFT. Apr 7, 2018 at 15:16
• The problem I have only with the 3D matrix FFT. I tried everything, changed and checked my code 100 times, so I am almost sure my code is OK, that is why I am asking about math transform. If anyone can assure me my math transform is good, then I know I have some mistake in the code. But it would be stupid to work hours with the code trying to implement math equation which is wrong. Apr 7, 2018 at 15:22
• Of course I can post my code here also, but I was afraid the question would be too long. So first I decide to ask clear question about math. Apr 7, 2018 at 15:23
• If I recall the book section. They don’t go to 3D. They say that you can decompose the column or row DFTs in the same manner. From a programming perspective, you can call a DFT 3 and DFT 4 from a DFT 12.
– user28715
Apr 7, 2018 at 15:27

Sometimes just writing things out brute force can shed some light. I think you have express $M$ and $L$ in terms of their factors. To be honest,I'm perfectly happy using FFTW. I haven't rolled my own DFTs for a very long time.

So I wrote a little Matlab script for the 12 point transform in terms of 4 point and 3 points

with input

intxt = {'x[0]'} {'x[1]'} {'x[2]'} {'x[3]'} {'x[4]'} {'x[5]'} {'x[6]'} {'x[7]'} {'x[8]'} {'x[9]'} {'x[10]'} {'x[11]'}

where we form the 2D matrix

  {'x[0]'}    {'x[1]' }    {'x[2]' }
{'x[3]'}    {'x[4]' }    {'x[5]' }
{'x[6]'}    {'x[7]' }    {'x[8]' }
{'x[9]'}    {'x[10]'}    {'x[11]'}

clear all
M=12;
in=[0:(M-1)];
intxt=num2cell(in);
intxt=cellfun(@(x) ['x[',int2str(x),']'],intxt,'UniformOutput',false)

zzz=DFT2D(intxt,3,4)

%%
function y=txtdft(x)  % 1d dft
y=[];
n=length(x);
for k=1:n
yy=[];
for i=1:n
if i < n
yy=[yy strcat(['W_{',int2str(n),'}^{',...
int2str(i-1),',',int2str(k-1),'}'],strcat(x{i},'+'))];
else
yy=[yy strcat(['W_{',int2str(n),'}^{',int2str(i-1),',...
',int2str(k-1),'}'],strcat(x{i}))];
end
% yy=['(',yy,')'];
end
y=[y {yy}];
end
end
%%

function [y]=twiddle(x)
[m,n]=size(x);
y=cell(m,n);
D=m*n;
for i=1:m
for k=1:n
twid= ['W_{',int2str(D),'}^{', ...
int2str(i-1),',',int2str(k-1),'}'];
y(i,k)={[ twid, '(',x{i,k},')']};
end
end
end
%%

function [y]=DFT2D( x, mm,nn )
y=cell(size(x));
length(x)
assert(length(x) == mm*nn,'not proper');
x12d=reshape(x,mm,nn)'   % mm = number colums nn=number rows
y_out=cell(size(x12d));
for i=1:mm
x11=x12d(:,i);
y_out(:,i)= txtdft(x11)';
end
zz=twiddle(y_out)';
y_out2=cell(size(zz));
for i=1:nn
x11=zz(:,i);
y_out2(:,i)= txtdft(x11)';
end
y=reshape(y_out2,length(x),1);   % read out rowszzz

end


And Pretty Printing the terms:

$$\begin{multline} Y[0]=W_{3}^{0,0}W_{12}^{0,0}(W_{4}^{0,0}x[0]+W_{4}^{1,0}x[3]+W_{4}^{2,0}x[6]+W_{4}^{3,0}x[9])+W_{3}^{1,0}W_{12}^{0,1}(W_{4}^{0,0}x[1]+W_{4}^{1,0}x[4]+W_{4}^{2,0}x[7]+W_{4}^{3,0}x[10])+W_{3}^{2,0}W_{12}^{0,2}(W_{4}^{0,0}x[2]+W_{4}^{1,0}x[5]+W_{4}^{2,0}x[8]+W_{4}^{3,0}x[11]) \end{multline}$$ $$\begin{multline} Y[1]=W_{3}^{0,1}W_{12}^{0,0}(W_{4}^{0,0}x[0]+W_{4}^{1,0}x[3]+W_{4}^{2,0}x[6]+W_{4}^{3,0}x[9])+W_{3}^{1,1}W_{12}^{0,1}(W_{4}^{0,0}x[1]+W_{4}^{1,0}x[4]+W_{4}^{2,0}x[7]+W_{4}^{3,0}x[10])+W_{3}^{2,1}W_{12}^{0,2}(W_{4}^{0,0}x[2]+W_{4}^{1,0}x[5]+W_{4}^{2,0}x[8]+W_{4}^{3,0}x[11]) \end{multline}$$ $$\begin{multline} Y[2]=W_{3}^{0,2}W_{12}^{0,0}(W_{4}^{0,0}x[0]+W_{4}^{1,0}x[3]+W_{4}^{2,0}x[6]+W_{4}^{3,0}x[9])+W_{3}^{1,2}W_{12}^{0,1}(W_{4}^{0,0}x[1]+W_{4}^{1,0}x[4]+W_{4}^{2,0}x[7]+W_{4}^{3,0}x[10])+W_{3}^{2,2}W_{12}^{0,2}(W_{4}^{0,0}x[2]+W_{4}^{1,0}x[5]+W_{4}^{2,0}x[8]+W_{4}^{3,0}x[11]) \end{multline}$$ $$\begin{multline} Y[3]= W_{3}^{0,0}W_{12}^{1,0}(W_{4}^{0,1}x[0]+W_{4}^{1,1}x[3]+W_{4}^{2,1}x[6]+W_{4}^{3,1}x[9])+W_{3}^{1,0}W_{12}^{1,1}(W_{4}^{0,1}x[1]+W_{4}^{1,1}x[4]+W_{4}^{2,1}x[7]+W_{4}^{3,1}x[10])+W_{3}^{2,0}W_{12}^{1,2}(W_{4}^{0,1}x[2]+W_{4}^{1,1}x[5]+W_{4}^{2,1}x[8]+W_{4}^{3,1}x[11]) \end{multline}$$ $$\begin{multline} Y[4]=W_{3}^{0,1}W_{12}^{1,0}(W_{4}^{0,1}x[0]+W_{4}^{1,1}x[3]+W_{4}^{2,1}x[6]+W_{4}^{3,1}x[9])+W_{3}^{1,1}W_{12}^{1,1}(W_{4}^{0,1}x[1]+W_{4}^{1,1}x[4]+W_{4}^{2,1}x[7]+W_{4}^{3,1}x[10])+W_{3}^{2,1}W_{12}^{1,2}(W_{4}^{0,1}x[2]+W_{4}^{1,1}x[5]+W_{4}^{2,1}x[8]+W_{4}^{3,1}x[11]) \end{multline}$$ $$\begin{multline} Y[5]=W_{3}^{0,2}W_{12}^{1,0}(W_{4}^{0,1}x[0]+W_{4}^{1,1}x[3]+W_{4}^{2,1}x[6]+W_{4}^{3,1}x[9])+W_{3}^{1,2}W_{12}^{1,1}(W_{4}^{0,1}x[1]+W_{4}^{1,1}x[4]+W_{4}^{2,1}x[7]+W_{4}^{3,1}x[10])+W_{3}^{2,2}W_{12}^{1,2}(W_{4}^{0,1}x[2]+W_{4}^{1,1}x[5]+W_{4}^{2,1}x[8]+W_{4}^{3,1}x[11]) \end{multline}$$ $$\begin{multline} Y[6]= W_{3}^{0,0}W_{12}^{2,0}(W_{4}^{0,2}x[0]+W_{4}^{1,2}x[3]+W_{4}^{2,2}x[6]+W_{4}^{3,2}x[9])+W_{3}^{1,0}W_{12}^{2,1}(W_{4}^{0,2}x[1]+W_{4}^{1,2}x[4]+W_{4}^{2,2}x[7]+W_{4}^{3,2}x[10])+W_{3}^{2,0}W_{12}^{2,2}(W_{4}^{0,2}x[2]+W_{4}^{1,2}x[5]+W_{4}^{2,2}x[8]+W_{4}^{3,2}x[11]) \end{multline}$$ $$\begin{multline} Y[8]= W_{3}^{0,1}W_{12}^{2,0}(W_{4}^{0,2}x[0]+W_{4}^{1,2}x[3]+W_{4}^{2,2}x[6]+W_{4}^{3,2}x[9])+W_{3}^{1,1}W_{12}^{2,1}(W_{4}^{0,2}x[1]+W_{4}^{1,2}x[4]+W_{4}^{2,2}x[7]+W_{4}^{3,2}x[10])+W_{3}^{2,1}W_{12}^{2,2}(W_{4}^{0,2}x[2]+W_{4}^{1,2}x[5]+W_{4}^{2,2}x[8]+W_{4}^{3,2}x[11]) \end{multline}$$ $$\begin{multline} Y[8]=W_{3}^{0,2}W_{12}^{2,0}(W_{4}^{0,2}x[0]+W_{4}^{1,2}x[3]+W_{4}^{2,2}x[6]+W_{4}^{3,2}x[9])+W_{3}^{1,2}W_{12}^{2,1}(W_{4}^{0,2}x[1]+W_{4}^{1,2}x[4]+W_{4}^{2,2}x[7]+W_{4}^{3,2}x[10])+W_{3}^{2,2}W_{12}^{2,2}(W_{4}^{0,2}x[2]+W_{4}^{1,2}x[5]+W_{4}^{2,2}x[8]+W_{4}^{3,2}x[11]) \end{multline}$$ $$\begin{multline} Y[9]=W_{3}^{0,0}W_{12}^{3,0}(W_{4}^{0,3}x[0]+W_{4}^{1,3}x[3]+W_{4}^{2,3}x[6]+W_{4}^{3,3}x[9])+W_{3}^{1,0}W_{12}^{3,1}(W_{4}^{0,3}x[1]+W_{4}^{1,3}x[4]+W_{4}^{2,3}x[7]+W_{4}^{3,3}x[10])+W_{3}^{2,0}W_{12}^{3,2}(W_{4}^{0,3}x[2]+W_{4}^{1,3}x[5]+W_{4}^{2,3}x[8]+W_{4}^{3,3}x[11]) \end{multline}$$ $$\begin{multline} Y[10]=W_{3}^{0,1}W_{12}^{3,0}(W_{4}^{0,3}x[0]+W_{4}^{1,3}x[3]+W_{4}^{2,3}x[6]+W_{4}^{3,3}x[9])+W_{3}^{1,1}W_{12}^{3,1}(W_{4}^{0,3}x[1]+W_{4}^{1,3}x[4]+W_{4}^{2,3}x[7]+W_{4}^{3,3}x[10])+W_{3}^{2,1}W_{12}^{3,2}(W_{4}^{0,3}x[2]+W_{4}^{1,3}x[5]+W_{4}^{2,3}x[8]+W_{4}^{3,3}x[11]) \end{multline}$$ $$\begin{multline} Y[11]= W_{3}^{0,2}W_{12}^{3,0}(W_{4}^{0,3}x[0]+W_{4}^{1,3}x[3]+W_{4}^{2,3}x[6]+W_{4}^{3,3}x[9])+W_{3}^{1,2}W_{12}^{3,1}(W_{4}^{0,3}x[1]+W_{4}^{1,3}x[4]+W_{4}^{2,3}x[7]+W_{4}^{3,3}x[10])+W_{3}^{2,2}W_{12}^{3,2}(W_{4}^{0,3}x[2]+W_{4}^{1,3}x[5]+W_{4}^{2,3}x[8]+W_{4}^{3,3}x[11])\end{multline}$$

For 60 points the 2 d matrix is

 {'x[0]' }    {'x[1]' }    {'x[2]' }    {'x[3]' }    {'x[4]' }
{'x[5]' }    {'x[6]' }    {'x[7]' }    {'x[8]' }    {'x[9]' }
{'x[10]'}    {'x[11]'}    {'x[12]'}    {'x[13]'}    {'x[14]'}
{'x[15]'}    {'x[16]'}    {'x[17]'}    {'x[18]'}    {'x[19]'}
{'x[20]'}    {'x[21]'}    {'x[22]'}    {'x[23]'}    {'x[24]'}
{'x[25]'}    {'x[26]'}    {'x[27]'}    {'x[28]'}    {'x[29]'}
{'x[30]'}    {'x[31]'}    {'x[32]'}    {'x[33]'}    {'x[34]'}
{'x[35]'}    {'x[36]'}    {'x[37]'}    {'x[38]'}    {'x[39]'}
{'x[40]'}    {'x[41]'}    {'x[42]'}    {'x[43]'}    {'x[44]'}
{'x[45]'}    {'x[46]'}    {'x[47]'}    {'x[48]'}    {'x[49]'}
{'x[50]'}    {'x[51]'}    {'x[52]'}    {'x[53]'}    {'x[54]'}
{'x[55]'}    {'x[56]'}    {'x[57]'}    {'x[58]'}    {'x[59]'}


and we shove the first column into the script. the first column will look like: $$\begin{multline} Y[0,0]=W_{3}^{0,0}W_{12}^{0,0}(W_{4}^{0,0}x[0]+W_{4}^{1,0}x[15]+W_{4}^{2,0}x[30]+W_{4}^{3,0}x[45])+W_{3}^{1,0}W_{12}^{0,1}(W_{4}^{0,0}x[5]+W_{4}^{1,0}x[20]+W_{4}^{2,0}x[35]+W_{4}^{3,0}x[50])+W_{3}^{2,0}W_{12}^{0,2}(W_{4}^{0,0}x[10]+W_{4}^{1,0}x[25]+W_{4}^{2,0}x[40]+W_{4}^{3,0}x[55]) \end{multline}$$

I'm not going to go through the tedium of going through each column DFT, doing the twiddle one the 12 by 5 matrix, and 12 row 5 point DFTs,

but your welcome to modify it yourself

OK, I need to answer my own question. The math transform I described is GOOD, the problem was like always, stupid mistake in the code. But Stanleys answer helped me much. Thanks a lot.