function y = fft61(x,u,ip,op)
% y = fft61(x,u,ip,op)
% y  : the 61 point DFT of x 
% u  : a vector of precomputed multiplicative constants
% ip : input permutation
% op : ouput permutation

y = zeros(61,1);
x = x(ip);                                        % input permutation
x(2:61) = KRED([2,3,5],[2,1,1],3,x(2:61));        % reduction operations
y(1) = x(1)+x(2);                                 % DC term calculation
% -------------------- block : 1 -------------------------------------------------
y(2) = x(2)*u(1);                       
% -------------------- block : 2 -------------------------------------------------
y(3) = x(3)*u(2);                       
% -------------------- block : 4 -------------------------------------------------
v = ID2I(1,1,x(4:5));              % v = (I(1) kron D2 kron I(1)) * x(4:5)
v = v.*u(3:5);                          
y(4:5) = ID2tI(1,1,v);             % y(4:5) = (I(1) kron D2' kron I(1)) * v
% -------------------- block : 3 -------------------------------------------------
v = ID2I(1,1,x(6:7));              % v = (I(1) kron D2 kron I(1)) * x(6:7)
v = v.*u(6:8);                          
y(6:7) = ID2tI(1,1,v);             % y(6:7) = (I(1) kron D2' kron I(1)) * v
% -------------------- block : 6 = 2 * 3 -----------------------------------------
v = ID2I(1,1,x(8:9));              % v = (I(1) kron D2 kron I(1)) * x(8:9)
v = v.*u(9:11);                         
y(8:9) = ID2tI(1,1,v);             % y(8:9) = (I(1) kron D2' kron I(1)) * v
% -------------------- block : 12 = 4 * 3 ----------------------------------------
v = ID2I(1,2,x(10:13));            % v = (I(1) kron D2 kron I(2)) * x(10:13)
v = ID2I(3,1,v);                   % v = (I(3) kron D2 kron I(1)) * v
v = v.*u(12:20);                        
v = ID2tI(1,3,v);                  % v = (I(1) kron D2' kron I(3)) * v
y(10:13) = ID2tI(2,1,v);           % y(10:13) = (I(2) kron D2' kron I(1)) * v
% -------------------- block : 5 -------------------------------------------------
v = ID2I(1,2,x(14:17));            % v = (I(1) kron D2 kron I(2)) * x(14:17)
v = ID2I(3,1,v);                   % v = (I(3) kron D2 kron I(1)) * v
v = v.*u(21:29);                        
v = ID2tI(1,3,v);                  % v = (I(1) kron D2' kron I(3)) * v
y(14:17) = ID2tI(2,1,v);           % y(14:17) = (I(2) kron D2' kron I(1)) * v
% -------------------- block : 10 = 2 * 5 ----------------------------------------
v = ID2I(1,2,x(18:21));            % v = (I(1) kron D2 kron I(2)) * x(18:21)
v = ID2I(3,1,v);                   % v = (I(3) kron D2 kron I(1)) * v
v = v.*u(30:38);                        
v = ID2tI(1,3,v);                  % v = (I(1) kron D2' kron I(3)) * v
y(18:21) = ID2tI(2,1,v);           % y(18:21) = (I(2) kron D2' kron I(1)) * v
% -------------------- block : 20 = 4 * 5 ----------------------------------------
v = ID2I(1,4,x(22:29));            % v = (I(1) kron D2 kron I(4)) * x(22:29)
v = ID2I(3,2,v);                   % v = (I(3) kron D2 kron I(2)) * v
v = ID2I(9,1,v);                   % v = (I(9) kron D2 kron I(1)) * v
v = v.*u(39:65);                        
v = ID2tI(1,9,v);                  % v = (I(1) kron D2' kron I(9)) * v
v = ID2tI(2,3,v);                  % v = (I(2) kron D2' kron I(3)) * v
y(22:29) = ID2tI(4,1,v);           % y(22:29) = (I(4) kron D2' kron I(1)) * v
% -------------------- block : 15 = 3 * 5 ----------------------------------------
v = ID2I(1,4,x(30:37));            % v = (I(1) kron D2 kron I(4)) * x(30:37)
v = ID2I(3,2,v);                   % v = (I(3) kron D2 kron I(2)) * v
v = ID2I(9,1,v);                   % v = (I(9) kron D2 kron I(1)) * v
v = v.*u(66:92);                        
v = ID2tI(1,9,v);                  % v = (I(1) kron D2' kron I(9)) * v
v = ID2tI(2,3,v);                  % v = (I(2) kron D2' kron I(3)) * v
y(30:37) = ID2tI(4,1,v);           % y(30:37) = (I(4) kron D2' kron I(1)) * v
% -------------------- block : 30 = 2 * 3 * 5 ------------------------------------
v = ID2I(1,4,x(38:45));            % v = (I(1) kron D2 kron I(4)) * x(38:45)
v = ID2I(3,2,v);                   % v = (I(3) kron D2 kron I(2)) * v
v = ID2I(9,1,v);                   % v = (I(9) kron D2 kron I(1)) * v
v = v.*u(93:119);                       
v = ID2tI(1,9,v);                  % v = (I(1) kron D2' kron I(9)) * v
v = ID2tI(2,3,v);                  % v = (I(2) kron D2' kron I(3)) * v
y(38:45) = ID2tI(4,1,v);           % y(38:45) = (I(4) kron D2' kron I(1)) * v
% -------------------- block : 60 = 4 * 3 * 5 ------------------------------------
v = ID2I(1,8,x(46:61));            % v = (I(1) kron D2 kron I(8)) * x(46:61)
v = ID2I(3,4,v);                   % v = (I(3) kron D2 kron I(4)) * v
v = ID2I(9,2,v);                   % v = (I(9) kron D2 kron I(2)) * v
v = ID2I(27,1,v);                  % v = (I(27) kron D2 kron I(1)) * v
v = v.*u(120:200);                      
v = ID2tI(1,27,v);                 % v = (I(1) kron D2' kron I(27)) * v
v = ID2tI(2,9,v);                  % v = (I(2) kron D2' kron I(9)) * v
v = ID2tI(4,3,v);                  % v = (I(4) kron D2' kron I(3)) * v
y(46:61) = ID2tI(8,1,v);           % y(46:61) = (I(8) kron D2' kron I(1)) * v
% --------------------------------------------------------------------------------
y(2) = y(1)+y(2);                                 % DC term calculation
y(2:61) = tKRED([2,3,5],[2,1,1],3,y(2:61));       % transpose reduction operations
y = y(op);                                        % output permutation

% For complex data - 
% Total Number of Real Multiplications : 400
% Total Number of Real Additions: 1908