0
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ Does your 12 point work on its own? $\endgroup$
    – user28715
    Apr 7, 2018 at 15:01
  • $\begingroup$ 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. $\endgroup$
    – pajczur
    Apr 7, 2018 at 15:16
  • $\begingroup$ 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. $\endgroup$
    – pajczur
    Apr 7, 2018 at 15:22
  • $\begingroup$ 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. $\endgroup$
    – pajczur
    Apr 7, 2018 at 15:23
  • $\begingroup$ 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. $\endgroup$
    – user28715
    Apr 7, 2018 at 15:27

2 Answers 2

0
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.