| << Chapter < Page | Chapter >> Page > |
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)) );
ADDS = 2 * Red_Adds;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:Ke_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);DC_ADDS = 2;
ADDS = ADDS + DC_ADDS;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))]; endif length(q)>1
SB = sprintf(' block : %d = %s ',B,bs);SC = SLINE;
SC(BL:BL-1+length(SB)) = SB;fprintf(fid,'%% %s\n',SC);
elseSB = sprintf(' block : %d ',B);
SC = SLINE;SC(BL:BL-1+length(SB)) = SB;
fprintf(fid,'%% %s\n',SC);end
if prod(q.^t) == 2S1 = sprintf('y(%d) = x(%d)*u(%d);',a+2,a+2,MULTS+1);
fprintf(fid,'%-40s\n',S1);Mk = 1;
elsed = []; 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];
endend
[g,C1]= cgc(Q,r,c,length(Q));
ADDS = ADDS + C1;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));ADDS = ADDS + C2;
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;
ADDS = 2* ADDS;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 permutationip = 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/2if 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
endfprintf(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:Pfprintf(fid,' %d\n',op(k));
endfprintf(fid,'];\n\n');fclose(fid);
Notification Switch
Would you like to follow the 'Automatic generation of prime length fft programs' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|