| << Chapter < Page | Chapter >> Page > |
The reduction matrix of
Equation 44 in Preliminaries is implemented
by
KRED which calls
RED .
Its transpose and inverse transpose are implemented by
tRED ,
tRED ,
itKRED and
itRED .
function x = KRED(P,E,K,x)
% x = KRED(P,E,K,x);% 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:0x(1:a*c*(p^(j+1))) = RED(p,a,c*(p^j),x(1:a*c*(p^(j+1))));
endend
function y = RED(p,a,c,x)
% y = RED(p,a,c,x);y = zeros(a*c*p,1);
for i = 0:c:(a-1)*cfor j = 0:c-1
y(i+j+1) = x(i*p+j+1);for k = 0:c:c*(p-2)
y(i+j+1) = y(i+j+1) + x(i*p+j+k+c+1);y(i*(p-1)+j+k+a*c+1) = x(i*p+j+k+1) - x(i*p+j+c*(p-1)+1);
endend
end
function x = tKRED(P,E,K,x)
% x = tKRED(P,E,K,x);% (transpose)
% P : P = [P(1),...,P(K)];
% E : E = [E(K),...,E(K)];
for i = K:-1:1a = 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 = 0:e-1
x(1:a*c*(p^(j+1))) = tRED(p,a,c*(p^j),x(1:a*c*(p^(j+1))));end
end
function y = tRED(p,a,c,x)
% y = tRED(p,a,c,x);% (transpose)
y = zeros(a*c*p,1);for i = 0:c:(a-1)*c
for j = 0:c-1y(i*p+j+c*(p-1)+1) = x(i+j+1);
for k = 0:c:c*(p-2)y(i*p+j+k+1) = x(i+j+1) + x(i*(p-1)+j+k+a*c+1);
y(i*p+j+c*(p-1)+1) = y(i*p+j+c*(p-1)+1) - x(i*(p-1)+j+k+a*c+1);end
endend
The operations of
and
are carried out by
ID2I and
ID3I .
Their transposes by
ID2tI and
ID3tI .
The functions
D2 and
D3 are listed in the appendix,
`Bilinear Forms for Linear Convolution.'Two versions of
ID2I are listed here.
One of them calls
D2 in a loop,
while the other version puts the
D2 code in the loop instead of using a function call.
There are several ways to implement theform
. But this is a simple
and straightforward method.It is modeled after
IAI in the text.
function y = ID2I(m,n,x)
y = zeros(m*n*3,1);v = 0:n:2*n;
u = 0:n:n;for i = 0:m-1
for j = 0:n-1y(v+i*3*n+j+1) = D2(x(u+i*2*n+j+1));
endend
function y = ID2I(m,n,x)
y = zeros(m*n*3,1);for i = 0:n:n*(m-1)
i2 = 2*i;i3 = 3*i;
for j = 1:nj2 = i2 + j;
j3 = i3 + j;y(j3) = x(j2);
y(n+j3) = x(n+j2);y(2*n+j3) = x(j2) + x(n+j2);
endend
function y = ID2tI(m,n,x)
y = zeros(m*n*2,1);v = 0:n:n;
u = 0:n:2*n;for i = 0:m-1
for j = 0:n-1y(v+i*2*n+j+1) = D2t(x(u+i*3*n+j+1));
endend
function y = ID3I(m,n,x)
y = zeros(m*n*5,1);v = 0:n:4*n;
u = 0:n:2*n;for i = 0:m-1
for j = 0:n-1y(v+i*5*n+j+1) = D3(x(u+i*3*n+j+1));
endend
function y = ID3tI(m,n,x)
y = zeros(m*n*3,1);v = 0:n:2*n;
u = 0:n:4*n;for i = 0:m-1
for j = 0:n-1y(v+i*3*n+j+1) = D3t(x(u+i*5*n+j+1));
endend
Notification Switch
Would you like to follow the 'Automatic generation of prime length fft programs' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|