# Connexions

You are here: Home » Content » Automatic Generation of Prime Length FFT Programs » Appendix: A Matlab Program for Generating Prime Length FFT Programs

### Lenses

What is a lens?

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

#### Affiliated with (What does "Affiliated with" mean?)

This content is either by members of the organizations listed or about topics related to the organizations listed. Click each link to see a list of all content affiliated with the organization.
• Rice Digital Scholarship

This collection is included in aLens by: Digital Scholarship at Rice University

Click the "Rice Digital Scholarship" link to see all content affiliated with them.

• Featured Content

This collection is included inLens: Connexions Featured Content
By: Connexions

Click the "Featured Content" link to see all content affiliated with them.

Click the tag icon to display tags associated with this content.

#### Also in these lenses

• UniqU content

This collection is included inLens: UniqU's lens
By: UniqU, LLC

Click the "UniqU content" link to see all content selected in this lens.

• Lens for Engineering

This module and collection are included inLens: Lens for Engineering
By: Sidney Burrus

Click the "Lens for Engineering" link to see all content selected in this lens.

### Recently Viewed

This feature requires Javascript to be enabled.

### Tags

(What is a tag?)

These tags come from the endorsement, affiliation, and other lenses that include this content.

Inside Collection (Technical report):

Report by: C. Sidney Burrus. E-mail the author

# Appendix: A Matlab Program for Generating Prime Length FFT Programs

Module by: Ivan Selesnick, C. Sidney Burrus. E-mail the authors


function [u,ip,op,ADDS,MULTS] = ff(p,e);
% [u,ip,op,ADDS,MULTS] = ff(p,e);
% u  : multiplicative constants
% ip : input permutation
% op : output permutation

K = length(p);
N = prod(p.^e);
P = N + 1;
[pr, ipr] = primitive_root(P);
Red_Adds = 2 * N * (K - sum(1./(p.^e)) );

FS = sprintf('fft%d.m',P);
fid = fopen(FS,'w');

fprintf(fid,'function y = fft%d(x,u,ip,op)\n',P);
fprintf(fid,'%% y = fft%d(x,u,ip,op)\n',P);
fprintf(fid,'%% y  : the %d point DFT of x \n',P);
fprintf(fid,'%% u  : a vector of precomputed multiplicative constants\n');
fprintf(fid,'%% ip : input permutation\n');
fprintf(fid,'%% op : ouput permutation\n');

Pstr = sprintf('[%d',p(1));
for k = 2:K, Pstr = [Pstr, sprintf(',%d',p(k))]; end
Pstr = [Pstr,']'];
Estr = sprintf('[%d',e(1));
for k = 2:K, Estr = [Estr, sprintf(',%d',e(k))]; end
Estr = [Estr,']'];
PEstr = sprintf('[%d',p(1)^e(1));
for k = 2:K, PEstr = [PEstr, sprintf(',%d',p(k)^e(k))]; end
PEstr = [PEstr,']'];

fprintf(fid,'\n');
S = sprintf('y = zeros(%d,1);\n',P);
fprintf(fid,S);
S1 = sprintf('x = x(ip);');
S2 = sprintf('%% input permutation\n');
fprintf(fid,'%-50s%s',S1,S2);
S1 = sprintf(['x(2:%d) = KRED(',Pstr,',',Estr,',%d,x(2:%d));'],P,K,P);
S2 = sprintf('%% reduction operations\n');
fprintf(fid,'%-50s%s',S1,S2);

e_table = [0:e(1)]';
a = e(1)+1;
for i = 2:K
e_table = [kron(ones(e(i)+1,1),e_table), kron([0:e(i)]',ones(a,1))];
a = a * (e(i)+1);
end
R = prod(e+1);

% ------------------------ MULTIPLICATIVE CONSTANTS ------------------------
k = rp(P,ipr,0:N);
I = sqrt(-1);
W = exp(-I*2*pi*k/P);
h = W(2:P);
h = h(N:-1:1);
h = pfp(p.^e,K,h);
h = itKRED(p,e,K,h);
u = h(1);

S1 = sprintf('y(1) = x(1)+x(2);');
S2 = sprintf('%% DC term calculation\n');
fprintf(fid,'%-50s%s',S1,S2);

SLINE = '--------------------------------------------------------------------------------';
SB = ' block : 1 ';
SC = SLINE;
BL = 21;
SC(BL:BL-1+length(SB)) = SB;
fprintf(fid,'%% %s\n',SC);
S1 = sprintf('y(2) = x(2)*u(1);');
fprintf(fid,'%-40s\n',S1);
a = 1;
MULTS = 1;
for i = 2:R
v = e_table(i,:);
f = find(v>0);
q = p(f);
t = v(f);
L = prod(q-1)*prod(q.^(t-1));

B = prod(q.^t);
bs = sprintf('%d',q(1)^t(1));
for k = 2:length(q), bs = [bs, sprintf(' * %d',q(k)^t(k))]; end
if length(q) > 1
SB = sprintf(' block : %d = %s ',B,bs);
SC = SLINE;
SC(BL:BL-1+length(SB)) = SB;
fprintf(fid,'%% %s\n',SC);
else
SB = sprintf(' block : %d ',B);
SC = SLINE;
SC(BL:BL-1+length(SB)) = SB;
fprintf(fid,'%% %s\n',SC);
end
if prod(q.^t) == 2
S1 = sprintf('y(%d) = x(%d)*u(%d);',a+2,a+2,MULTS+1);
fprintf(fid,'%-40s\n',S1);
Mk = 1;
else
d = []; r = []; c = []; Q = []; Qt = [];
for j = 1:length(q)
[dk,rk,ck,Qk,Qtk] = A_data(q(j)^t(j));
if dk > 1
d = [d dk]; r = [r rk]; c = [c ck]; Q = [Q Qk]; Qt = [Qt Qtk];
end
end
[g,C1] = cgc(Q,r,c,length(Q));
Mk = prod(r);
BEG = int2str(a+2); FIN = int2str(a+1+L);
XX = ['x(',BEG,':',FIN,')']; YY = 'v';
kpi(d,g,r,c,length(Q),YY,XX,fid);
S1 = ['v = v.*u(',int2str(MULTS+1),':',int2str(MULTS+Mk),');'];
fprintf(fid,'%-40s\n',S1);
[g,C2] = cgc(Qt,c,r,length(Q));
XX = 'v'; YY = ['y(',BEG,':',FIN,')'];
kpit(d,g,c,r,length(Q),YY,XX,fid);
end

c = [];
r = [];
lq = length(q);
for j = 1:lq
[fk,rk,ck] = C_data(q(j),t(j));
r = [r rk]; c = [c ck];
end
f = (q-1).*(q.^(t-1));
temp = Kcrot(q,t,lq,h(a+1:a+L));
temp = KFt(f,r,c,temp);
u = [u; temp(:)];
a = a + L;
MULTS = MULTS + Mk;
end

u(1) = u(1)-1;
fprintf(fid,'%% %s\n',SLINE);
S1 = sprintf('y(2) = y(1)+y(2);');
S2 = sprintf('%% DC term calculation\n');
fprintf(fid,'%-50s%s',S1,S2);
S1 = sprintf(['y(2:%d) = tKRED(',Pstr,',',Estr,',%d,y(2:%d));'],P,K,P);
S2 = sprintf('%% transpose reduction operations\n');
fprintf(fid,'%-50s%s',S1,S2);
S1 = sprintf('y = y(op);');
S2 = sprintf('%% output permutation\n');
fprintf(fid,'%-50s%s',S1,S2);
fprintf(fid,'\n');

MULTS = 2 * MULTS;
fprintf(fid,'%% For complex data - \n');
fprintf(fid,'%% Total Number of Real Multiplications : %d\n',MULTS);
fprintf(fid,'%% Total Number of Real Additions: %d\n\n',ADDS);
fclose(fid);

%%%%%%%%%%%%%%%%%%%% COMPUTE INPUT AND OUTPUT PERMUTATIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%

id = 1:P;   % identity permutation
ip = rp(P,pr,id);
ip(2:P) = pfp(p.^e,K,ip(2:P));

op = id;
op(2:P) = pfpt(p.^e,K,op(2:P));
op(2:P) = op(P:-1:2);
op = rpt(P,ipr,op);

%%%%%%%%%%%%%%%%% PUT MULTIPLICATIVE CONSTANTS AND PERMUTATIONS IN A FILE %%%%%%%%%%%%%%

CFS = sprintf('cap%d.m',P);
fid = fopen(CFS,'w');
fprintf(fid,'\n%% The multiplicative constants for the %d point FFT\n\n',P);
fprintf(fid,'I = sqrt(-1);\n');
fprintf(fid,'u = [\n');
for k = 1:MULTS/2
if abs(real(u(k))) < 0.000001
fprintf(fid,'%25.15f*I\n',imag(u(k)));
elseif abs(imag(u(k))) < 0.00001
fprintf(fid,'%25.15f\n',real(u(k)));
else
fprintf(fid,'%25.15f + %25.15f*I\n',real(u(k)),imag(u(k)));
end
end
fprintf(fid,'];\n\n');
fprintf(fid,'\n%% The input permutation for the %d point FFT\n\n',P);
fprintf(fid,'ip = [\n');
for k = 1:P
fprintf(fid,'   %d\n',ip(k));
end
fprintf(fid,'];\n\n');
fprintf(fid,'\n%% The output permutation for the %d point FFT\n\n',P);
fprintf(fid,'op = [\n');
for k = 1:P
fprintf(fid,'   %d\n',op(k));
end
fprintf(fid,'];\n\n');
fclose(fid);


The following programs print the program statements that carry out the operation IDkIIDkI and IDktIIDktI. They are modeled after kpi in the text.

function kpi(d,g,r,c,n,Y,X,fid)
% kpi(d,g,r,c,n,Y,X,fid);
% Kronecker Product : A(d(1)) kron ... kron A(d(n))
% g : permutation of 1,...,n
% r : [r(1),...,r(n)]
% c : [c(1),..,c(n)]
% r(i) : rows of A(d(i))
% c(i) : columns of A(d(i))
% n : number of terms
for i = 1:n
a = 1;
for k = 1:(g(i)-1)
if i > find(g==k)
a = a * r(k);
else
a = a * c(k);
end
end
b = 1;
for k = (g(i)+1):n
if i > find(g==k)
b = b * r(k);
else
b = b * c(k);
end
end
% Y = (I(a) kron A(d(g(i))) kron I(b)) * X;
if i == 1
S1 = sprintf([Y,' = ID%dI(%d,%d,',X,');    '],d(g(i)),a,b);
S2 = sprintf(['%% ',Y,' = (I(%d) kron D%d kron I(%d)) * ',X],a,d(g(i)),b);
fprintf(fid,'%-35s%s\n',S1,S2);
elseif d(g(i)) ~= 1
S1 = sprintf([Y,' = ID%dI(%d,%d,',Y,');    '],d(g(i)),a,b);
S2 = sprintf(['%% ',Y,' = (I(%d) kron D%d kron I(%d)) * ',Y],a,d(g(i)),b);
fprintf(fid,'%-35s%s\n',S1,S2);
end
end

function kpit(d,g,r,c,n,Y,X,fid)
% kpit(g,r,c,n,Y,X,fid);
% (transpose)
% Kronecker Product : A(d(1))' kron ... kron A(d(n))'
% g : permutation of 1,...,n
% r : [r(1),...,r(n)]
% c : [c(1),..,c(n)]
% r(i) : rows of A(d(i))'
% c(i) : columns of A(d(i))'
% n : number of terms

for i = 1:n
a = 1;
for k = 1:(g(i)-1)
if i > find(g==k)
a = a * r(k);
else
a = a * c(k);
end
end
b = 1;
for k = (g(i)+1):n
if i > find(g==k)
b = b * r(k);
else
b = b * c(k);
end
end
% x = (I(a) kron A(d(g(i)))'' kron I(b)) * x;
if i == n
S1 = sprintf([Y,' = ID%dtI(%d,%d,',X,');    '],d(g(i)),a,b);
S2 = sprintf(['%% ',Y,' = (I(%d) kron D%d'' kron I(%d)) * ',X],a,d(g(i)),b);
fprintf(fid,'%-35s%s\n',S1,S2);
elseif d(g(i)) ~= 1
S1 = sprintf([X,' = ID%dtI(%d,%d,',X,');    '],d(g(i)),a,b);
S2 = sprintf(['%% ',X,' = (I(%d) kron D%d'' kron I(%d)) * ',X],a,d(g(i)),b);
fprintf(fid,'%-35s%s\n',S1,S2);
end
end


## Programs for Computing Multiplicative Constants

The following programs carry out the operation of Fd1FdKFd1FdK where FF is the reconstruction matrix in a linear convolution algorithm. See the appendix, Bilinear Forms for Linear Convolution.'

function u = KFt(f,r,c,u)
% u = (F^t kron ... kron F^t)*u
% (transpose)
% f = [f(1),...,f(K)]
% r : r(i) = rows of F(i)
% c : c(i) = columns of F(i)
% u : length(u) = prod(c);
K = length(f);
for i = 1:K
m = prod(c(1:i-1));
n = prod(r(i+1:K));
u = IFtI(f(i),r(i),c(i),m,n,u);
end

function y = IFtI(s,r,c,m,n,x);
% y = (I(m) kron F(s)^t kron I(n))*x
% (transpose)
% r : rows of F(s)
% c : columns of F(s)
v = 0:n:n*(c-1);
u = 0:n:n*(r-1);
for i = 0:m-1
for j = 0:n-1
y(v+i*c*n+j+1) = Ftop(s,x(u+i*r*n+j+1));
end
end

function y = Ftop(k,x)
if k == 1, y = x;
elseif k == 2, y = F2t(x);
elseif k == 3, y = F3t(x);
elseif k == 4, y = F4t(x);
elseif k == 6, y = F6t(x);
elseif k == 8, y = F8t(x);
elseif k == 18, y = F18t(x);
end


The following programs carry out the operation of Gp1e1GpKeKGp1e1GpKeK were GG is given by Equation 13 and Equation 14 from Bilinear Forms for Circular Convolution.

function x = Kcrot(p,e,K,x)
% Kronecker product of Cyclotomic Reduction Operations.
% x = (G(p(1)^e(1)) kron ... kron G(p(K)^(K)))^t*x
% (transpose)
% p : p = [p(1),...,p(K)];
% e : e = [e(1),...,e(K)];
a = (p-1).*((p).^(e-1));
r = a;		% r(i) = number of rows of G(i)
c = 2*a-1;	% c(i) = number of columns of G(i)
m = 1;
n = prod(r);
for i = 1:K
n = n / r(i);
x = IcrotI(p(i),e(i),m,n,x);
m = m * c(i);
end

function y = IcrotI(p,e,m,n,x)
%  y = (eye(m) kron G(p^e)^t kron eye(n))*x
%  (transpose)
a = (p-1)*(p^(e-1));
c = a;
r = 2*a-1;
y = zeros(r*m*n,1);
v = 0:n:(r-1)*n;
u = 0:n:(c-1)*n;
for i = 0:m-1
for j = 0:n-1
y(v+i*r*n+j+1) = crot(p,e,x(u+i*c*n+j+1));
end
end

function y = crot(p,e,x)
% y = crot(p,x)
% cyclotomic reduction matrix (transpose)
% length(x) == 2*n-1
% length(y) == n
% where n = (p-1)*(p^(e-1))
n = (p-1)*(p^(e-1));
y = zeros(2*n-1,1);
if p == 2
n = p^(e-1);
y(1:n) = x;
y(n+1:2*n-1) = -x(1:n-1);
else
y(1:n) = x;
L = p^(e-1);
y(n+1:n+L) = -x(1:L);
a = L;
for k = 2:p-1
y(n+1:n+L) = y(n+1:n+L) - x(a+1:a+L);
a = a + L;
end
b = 2*n-1 - p*(p^(e-1));
y(p*L+1:p*L+b) = x(1:b);
end


The following programs tell the programs for code generation relevant information about the bilinear forms for cyclotomic convolution. Specifically, they indicates the linear convolution out of which these cyclotomic convolution are composed, and the dimensions of the corresponding matrices. See the appendix Bilinear Forms for Linear Convolution.

function [d,r,c,Q,Qt] = A_data(n)
% A : A matrix in bilinear form for cyclotomic convolution
% d : linear convolution modules used
% r : rows
% c : columns
% Q : Q(i) = cost associated with D(d(i))
% Qt : Qt(i) = cost associated with D(d(i))'
if n == 2, d = [1];
elseif n == 4, d = [2];
elseif n == 8, d = [2 2];
elseif n == 16, d = [2 2 2];
elseif n == 3, d = [2];
elseif n == 9, d = [2 3];
elseif n == 27, d = [2 3 3];
elseif n == 5, d = [2 2];
elseif n == 7, d = [2 3];
end
r = []; c = []; Q = []; Qt = [];
for k = 1:length(d)
[rk, ck, Qk, Qtk] = D_data(d(k));
r = [r rk]; c = [c ck]; Q = [Q Qk]; Qt = [Qt Qtk];
end

function [r,c,Q,Qt] = D_data(d);
% D : D matrix in bilinear form for linear convolution
% r : rows
% c : columns
% Q : cost associated with D(d)
% Qt : cost associated with D(d)'
if d == 1, r = 1; c = 1; Q = 0; Qt = 0;
elseif d == 2, r = 3; c = 2; Q = 1; Qt = 2;
elseif d == 3, r = 5; c = 3; Q = 7; Qt = 9;
end

function [f,r,c] = C_data(p,e)
% f : length of linear convolution
% r : rows
% c : columns
f = prod((p-1).*(p.^(e-1)));
% (Euler Totient Function)
r = 2*f-1;
c = F_data(f);

function c = F_data(n)
% c : columns of F matrix
if n == 1, c = 1;
elseif n == 2, c = 3;
elseif n == 4, c = 9;
elseif n == 8, c = 27;
elseif n == 3, c = 5;
elseif n == 6, c = 15;
elseif n == 18, c = 75;
end


### Programs for Inverse Transpose Reduction Operations

function x = itKRED(P,E,K,x)
% x = itKRED(P,E,K,x);
% (inverse transpose)
% P : P = [P(1),...,P(K)];
% E : E = [E(K),...,E(K)];
for i = 1:K
a = prod(P(1:i-1).^E(1:i-1));
c = prod(P(i+1:K).^E(i+1:K));
p = P(i);
e = E(i);
for j = e-1:-1:0
x(1:a*c*(p^(j+1))) = itRED(p,a,c*(p^j),x(1:a*c*(p^(j+1))));
end
end

function y = itRED(p,a,c,x)
% y = itRED(p,a,c,x);
% (inverse transpose)
y = zeros(a*c*p,1);
for i = 0:c:(a-1)*c
for j = 0:c-1
A = x(i*p+j+1);
for k = 0:c:c*(p-2)
A = A + x(i*p+j+k+c+1);
end
y(i+j+1) = A;
for k = 0:c:c*(p-2)
y(i*(p-1)+j+k+a*c+1) = p*x(i*p+j+k+1) - A;
end
end
end
y = y/p;


## Programs for Permutations

The permutation of Equation 18 from Preliminaries is implemented by pfp . It calls the function pfp2I . The transpose is implemented by pfpt and it calls pfpt2I .

function x = pfp(n,K,x)
% x = P(n(1),...,n(K)) * x
% n = [n(1),...,n(K)];
% length(x) = prod(n(1),...,n(K))
a = prod(n);
s = 1;
for i = K:-1:2
a = a / n(i);
x = pfp2I(a,n(i),s,x);
s = s * n(i);
end

function y = pfp2I(a,b,s,x)
% y = kron(P(a,b),I(s)) * x;
% length(x) = a*b*s
n = a * b;
y = zeros(n*s,1);
k1 = 0;
k2 = 0;
for k = 0:n-1
i1 = s * (k1 + b * k2);
i2 = s * k;
for i = 1:s
y(i1 + i) = x(i2 + i);
end
k1 = k1 + 1;
k2 = k2 + 1;
if k1 >= b
k1 = k1 - b;
end
if k2 >= a
k2 = k2 - a;
end
end

function x = pfpt(n,K,x)
% x = P(n(1),...,n(K))' * x
% (tanspose)
% n = [n(1),...,n(K)];
% length(x) = prod(n(1),...,n(K))
% a = prod(n);
a = n(1);
s = prod(n(2:K));
for i = 2:K
s = s / n(i);
x = pfpt2I(a,n(i),s,x);
a = a * n(i);
end

function y = pfpt2I(a,b,s,x)
% y = P(a,b)' kron I(s) * x;
% (transpose)
% length(x) = a*b*s
n = a * b;
y = zeros(n*s,1);
k1 = 0;
k2 = 0;
for k = 0:n-1
i1 = s * (k1 + b * k2);
i2 = s * k;
for i = 1:s
y(i2 + i) = x(i1 + i);
end
k1 = k1 + 1;
k2 = k2 + 1;
if k1 >= b
k1 = k1 - b;
end
if k2 >= a
k2 = k2 - a;
end
end


The following Matlab programs implement Rader's permutation and its transpose. They require the primitive root to be passed to them as an argument.

function y = rp(p,r,x)
% p : prime
% r : a primitive root of p
% x : length(x) == p
a = 1;
y = zeros(p,1);
y(1) = x(1);
for k = 2:p
y(k) = x(a+1);
a = rem(a*r,p);
end

function y = rpt(p,r,x)
% (transpose)
% p : prime
% r : a primitive root of p
% x : length(x) == p
a = 1;
y = zeros(p,1);
y(1) = x(1);
for k = 2:p
y(a+1) = x(k);
a = rem(a*r,p);
end

function [R, R_inv] = primitive_root(N)

%  function [R, R_inv] = primitive_root(N)
%  Ivan Selesnick
%  N is assumed to be prime.  This function returns R,
%  the smallest	primitive root of N, and R_inv, the
%  inverse of R modulo N.

m = 0:(N-2);
for x = 1:(N-1)
if ( 1:(N-1) == sort(rem2(x,m,N)) )
R = x;
break
end
end
for x = 1:N
if rem(x*R,N) == 1
R_inv = x;
break
end
end
`

## Content actions

### Give feedback:

#### Collection as:

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

For detailed instructions on how to download this content's EPUB to your specific device, click the "(?)" link.

#### Module as:

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

For detailed instructions on how to download this content's EPUB to your specific device, click the "(?)" link.

#### Collection to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

#### Module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks