Google Git
Sign in
android / platform / external / fft2d / 878be62896befc91a6268bf48287f8ea5eef8e90 / . / src / fft2d / fft2d / fftsg.f
blob: 6fb26c496119428cdbd010e6bff70f1367682c10 [file] [log] [blame]
! Fast Fourier/Cosine/Sine Transform
! dimension :one
! data length :power of 2
! decimation :frequency
! radix :split-radix
! data :inplace
! table :use
! subroutines
! cdft: Complex Discrete Fourier Transform
! rdft: Real Discrete Fourier Transform
! ddct: Discrete Cosine Transform
! ddst: Discrete Sine Transform
! dfct: Cosine Transform of RDFT (Real Symmetric DFT)
! dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
!
!
! -------- Complex DFT (Discrete Fourier Transform) --------
! [definition]
! <case1>
! X(k) = sum_j=0^n-1 x(j)*exp(2*pi*i*j*k/n), 0<=k<n
! <case2>
! X(k) = sum_j=0^n-1 x(j)*exp(-2*pi*i*j*k/n), 0<=k<n
! (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
! [usage]
! <case1>
! ip(0) = 0 ! first time only
! call cdft(2*n, 1, a, ip, w)
! <case2>
! ip(0) = 0 ! first time only
! call cdft(2*n, -1, a, ip, w)
! [parameters]
! 2*n :data length (integer)
! n >= 1, n = power of 2
! a(0:2*n-1) :input/output data (real*8)
! input data
! a(2*j) = Re(x(j)),
! a(2*j+1) = Im(x(j)), 0<=j<n
! output data
! a(2*k) = Re(X(k)),
! a(2*k+1) = Im(X(k)), 0<=k<n
! ip(0:*) :work area for bit reversal (integer)
! length of ip >= 2+sqrt(n)
! strictly,
! length of ip >=
! 2+2**(int(log(n+0.5)/log(2.0))/2).
! ip(0),ip(1) are pointers of the cos/sin table.
! w(0:n/2-1) :cos/sin table (real*8)
! w(),ip() are initialized if ip(0) = 0.
! [remark]
! Inverse of
! call cdft(2*n, -1, a, ip, w)
! is
! call cdft(2*n, 1, a, ip, w)
! do j = 0, 2 * n - 1
! a(j) = a(j) / n
! end do
! .
!
!
! -------- Real DFT / Inverse of Real DFT --------
! [definition]
! <case1> RDFT
! R(k) = sum_j=0^n-1 a(j)*cos(2*pi*j*k/n), 0<=k<=n/2
! I(k) = sum_j=0^n-1 a(j)*sin(2*pi*j*k/n), 0<k<n/2
! <case2> IRDFT (excluding scale)
! a(k) = (R(0) + R(n/2)*cos(pi*k))/2 +
! sum_j=1^n/2-1 R(j)*cos(2*pi*j*k/n) +
! sum_j=1^n/2-1 I(j)*sin(2*pi*j*k/n), 0<=k<n
! [usage]
! <case1>
! ip(0) = 0 ! first time only
! call rdft(n, 1, a, ip, w)
! <case2>
! ip(0) = 0 ! first time only
! call rdft(n, -1, a, ip, w)
! [parameters]
! n :data length (integer)
! n >= 2, n = power of 2
! a(0:n-1) :input/output data (real*8)
! <case1>
! output data
! a(2*k) = R(k), 0<=k<n/2
! a(2*k+1) = I(k), 0<k<n/2
! a(1) = R(n/2)
! <case2>
! input data
! a(2*j) = R(j), 0<=j<n/2
! a(2*j+1) = I(j), 0<j<n/2
! a(1) = R(n/2)
! ip(0:*) :work area for bit reversal (integer)
! length of ip >= 2+sqrt(n/2)
! strictly,
! length of ip >=
! 2+2**(int(log(n/2+0.5)/log(2.0))/2).
! ip(0),ip(1) are pointers of the cos/sin table.
! w(0:n/2-1) :cos/sin table (real*8)
! w(),ip() are initialized if ip(0) = 0.
! [remark]
! Inverse of
! call rdft(n, 1, a, ip, w)
! is
! call rdft(n, -1, a, ip, w)
! do j = 0, n - 1
! a(j) = a(j) * 2 / n
! end do
! .
!
!
! -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
! [definition]
! <case1> IDCT (excluding scale)
! C(k) = sum_j=0^n-1 a(j)*cos(pi*j*(k+1/2)/n), 0<=k<n
! <case2> DCT
! C(k) = sum_j=0^n-1 a(j)*cos(pi*(j+1/2)*k/n), 0<=k<n
! [usage]
! <case1>
! ip(0) = 0 ! first time only
! call ddct(n, 1, a, ip, w)
! <case2>
! ip(0) = 0 ! first time only
! call ddct(n, -1, a, ip, w)
! [parameters]
! n :data length (integer)
! n >= 2, n = power of 2
! a(0:n-1) :input/output data (real*8)
! output data
! a(k) = C(k), 0<=k<n
! ip(0:*) :work area for bit reversal (integer)
! length of ip >= 2+sqrt(n/2)
! strictly,
! length of ip >=
! 2+2**(int(log(n/2+0.5)/log(2.0))/2).
! ip(0),ip(1) are pointers of the cos/sin table.
! w(0:n*5/4-1) :cos/sin table (real*8)
! w(),ip() are initialized if ip(0) = 0.
! [remark]
! Inverse of
! call ddct(n, -1, a, ip, w)
! is
! a(0) = a(0) / 2
! call ddct(n, 1, a, ip, w)
! do j = 0, n - 1
! a(j) = a(j) * 2 / n
! end do
! .
!
!
! -------- DST (Discrete Sine Transform) / Inverse of DST --------
! [definition]
! <case1> IDST (excluding scale)
! S(k) = sum_j=1^n A(j)*sin(pi*j*(k+1/2)/n), 0<=k<n
! <case2> DST
! S(k) = sum_j=0^n-1 a(j)*sin(pi*(j+1/2)*k/n), 0<k<=n
! [usage]
! <case1>
! ip(0) = 0 ! first time only
! call ddst(n, 1, a, ip, w)
! <case2>
! ip(0) = 0 ! first time only
! call ddst(n, -1, a, ip, w)
! [parameters]
! n :data length (integer)
! n >= 2, n = power of 2
! a(0:n-1) :input/output data (real*8)
! <case1>
! input data
! a(j) = A(j), 0<j<n
! a(0) = A(n)
! output data
! a(k) = S(k), 0<=k<n
! <case2>
! output data
! a(k) = S(k), 0<k<n
! a(0) = S(n)
! ip(0:*) :work area for bit reversal (integer)
! length of ip >= 2+sqrt(n/2)
! strictly,
! length of ip >=
! 2+2**(int(log(n/2+0.5)/log(2.0))/2).
! ip(0),ip(1) are pointers of the cos/sin table.
! w(0:n*5/4-1) :cos/sin table (real*8)
! w(),ip() are initialized if ip(0) = 0.
! [remark]
! Inverse of
! call ddst(n, -1, a, ip, w)
! is
! a(0) = a(0) / 2
! call ddst(n, 1, a, ip, w)
! do j = 0, n - 1
! a(j) = a(j) * 2 / n
! end do
! .
!
!
! -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
! [definition]
! C(k) = sum_j=0^n a(j)*cos(pi*j*k/n), 0<=k<=n
! [usage]
! ip(0) = 0 ! first time only
! call dfct(n, a, t, ip, w)
! [parameters]
! n :data length - 1 (integer)
! n >= 2, n = power of 2
! a(0:n) :input/output data (real*8)
! output data
! a(k) = C(k), 0<=k<=n
! t(0:n/2) :work area (real*8)
! ip(0:*) :work area for bit reversal (integer)
! length of ip >= 2+sqrt(n/4)
! strictly,
! length of ip >=
! 2+2**(int(log(n/4+0.5)/log(2.0))/2).
! ip(0),ip(1) are pointers of the cos/sin table.
! w(0:n*5/8-1) :cos/sin table (real*8)
! w(),ip() are initialized if ip(0) = 0.
! [remark]
! Inverse of
! a(0) = a(0) / 2
! a(n) = a(n) / 2
! call dfct(n, a, t, ip, w)
! is
! a(0) = a(0) / 2
! a(n) = a(n) / 2
! call dfct(n, a, t, ip, w)
! do j = 0, n
! a(j) = a(j) * 2 / n
! end do
! .
!
!
! -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
! [definition]
! S(k) = sum_j=1^n-1 a(j)*sin(pi*j*k/n), 0<k<n
! [usage]
! ip(0) = 0 ! first time only
! call dfst(n, a, t, ip, w)
! [parameters]
! n :data length + 1 (integer)
! n >= 2, n = power of 2
! a(0:n-1) :input/output data (real*8)
! output data
! a(k) = S(k), 0<k<n
! (a(0) is used for work area)
! t(0:n/2-1) :work area (real*8)
! ip(0:*) :work area for bit reversal (integer)
! length of ip >= 2+sqrt(n/4)
! strictly,
! length of ip >=
! 2+2**(int(log(n/4+0.5)/log(2.0))/2).
! ip(0),ip(1) are pointers of the cos/sin table.
! w(0:n*5/8-1) :cos/sin table (real*8)
! w(),ip() are initialized if ip(0) = 0.
! [remark]
! Inverse of
! call dfst(n, a, t, ip, w)
! is
! call dfst(n, a, t, ip, w)
! do j = 1, n - 1
! a(j) = a(j) * 2 / n
! end do
! .
!
!
! Appendix :
! The cos/sin table is recalculated when the larger table required.
! w() and ip() are compatible with all routines.
!
!
subroutine cdft(n, isgn, a, ip, w)
integer n, isgn, ip(0 : *), nw
real*8 a(0 : n - 1), w(0 : *)
nw = ip(0)
if (n .gt. 4 * nw) then
nw = n / 4
call makewt(nw, ip, w)
end if
if (isgn .ge. 0) then
call cftfsub(n, a, ip, nw, w)
else
call cftbsub(n, a, ip, nw, w)
end if
end
!
subroutine rdft(n, isgn, a, ip, w)
integer n, isgn, ip(0 : *), nw, nc
real*8 a(0 : n - 1), w(0 : *), xi
nw = ip(0)
if (n .gt. 4 * nw) then
nw = n / 4
call makewt(nw, ip, w)
end if
nc = ip(1)
if (n .gt. 4 * nc) then
nc = n / 4
call makect(nc, ip, w(nw))
end if
if (isgn .ge. 0) then
if (n .gt. 4) then
call cftfsub(n, a, ip, nw, w)
call rftfsub(n, a, nc, w(nw))
else if (n .eq. 4) then
call cftfsub(n, a, ip, nw, w)
end if
xi = a(0) - a(1)
a(0) = a(0) + a(1)
a(1) = xi
else
a(1) = 0.5d0 * (a(0) - a(1))
a(0) = a(0) - a(1)
if (n .gt. 4) then
call rftbsub(n, a, nc, w(nw))
call cftbsub(n, a, ip, nw, w)
else if (n .eq. 4) then
call cftbsub(n, a, ip, nw, w)
end if
end if
end
!
subroutine ddct(n, isgn, a, ip, w)
integer n, isgn, ip(0 : *), j, nw, nc
real*8 a(0 : n - 1), w(0 : *), xr
nw = ip(0)
if (n .gt. 4 * nw) then
nw = n / 4
call makewt(nw, ip, w)
end if
nc = ip(1)
if (n .gt. nc) then
nc = n
call makect(nc, ip, w(nw))
end if
if (isgn .lt. 0) then
xr = a(n - 1)
do j = n - 2, 2, -2
a(j + 1) = a(j) - a(j - 1)
a(j) = a(j) + a(j - 1)
end do
a(1) = a(0) - xr
a(0) = a(0) + xr
if (n .gt. 4) then
call rftbsub(n, a, nc, w(nw))
call cftbsub(n, a, ip, nw, w)
else if (n .eq. 4) then
call cftbsub(n, a, ip, nw, w)
end if
end if
call dctsub(n, a, nc, w(nw))
if (isgn .ge. 0) then
if (n .gt. 4) then
call cftfsub(n, a, ip, nw, w)
call rftfsub(n, a, nc, w(nw))
else if (n .eq. 4) then
call cftfsub(n, a, ip, nw, w)
end if
xr = a(0) - a(1)
a(0) = a(0) + a(1)
do j = 2, n - 2, 2
a(j - 1) = a(j) - a(j + 1)
a(j) = a(j) + a(j + 1)
end do
a(n - 1) = xr
end if
end
!
subroutine ddst(n, isgn, a, ip, w)
integer n, isgn, ip(0 : *), j, nw, nc
real*8 a(0 : n - 1), w(0 : *), xr
nw = ip(0)
if (n .gt. 4 * nw) then
nw = n / 4
call makewt(nw, ip, w)
end if
nc = ip(1)
if (n .gt. nc) then
nc = n
call makect(nc, ip, w(nw))
end if
if (isgn .lt. 0) then
xr = a(n - 1)
do j = n - 2, 2, -2
a(j + 1) = -a(j) - a(j - 1)
a(j) = a(j) - a(j - 1)
end do
a(1) = a(0) + xr
a(0) = a(0) - xr
if (n .gt. 4) then
call rftbsub(n, a, nc, w(nw))
call cftbsub(n, a, ip, nw, w)
else if (n .eq. 4) then
call cftbsub(n, a, ip, nw, w)
end if
end if
call dstsub(n, a, nc, w(nw))
if (isgn .ge. 0) then
if (n .gt. 4) then
call cftfsub(n, a, ip, nw, w)
call rftfsub(n, a, nc, w(nw))
else if (n .eq. 4) then
call cftfsub(n, a, ip, nw, w)
end if
xr = a(0) - a(1)
a(0) = a(0) + a(1)
do j = 2, n - 2, 2
a(j - 1) = -a(j) - a(j + 1)
a(j) = a(j) - a(j + 1)
end do
a(n - 1) = -xr
end if
end
!
subroutine dfct(n, a, t, ip, w)
integer n, ip(0 : *), j, k, l, m, mh, nw, nc
real*8 a(0 : n), t(0 : n / 2), w(0 : *), xr, xi, yr, yi
nw = ip(0)
if (n .gt. 8 * nw) then
nw = n / 8
call makewt(nw, ip, w)
end if
nc = ip(1)
if (n .gt. 2 * nc) then
nc = n / 2
call makect(nc, ip, w(nw))
end if
m = n / 2
yi = a(m)
xi = a(0) + a(n)
a(0) = a(0) - a(n)
t(0) = xi - yi
t(m) = xi + yi
if (n .gt. 2) then
mh = m / 2
do j = 1, mh - 1
k = m - j
xr = a(j) - a(n - j)
xi = a(j) + a(n - j)
yr = a(k) - a(n - k)
yi = a(k) + a(n - k)
a(j) = xr
a(k) = yr
t(j) = xi - yi
t(k) = xi + yi
end do
t(mh) = a(mh) + a(n - mh)
a(mh) = a(mh) - a(n - mh)
call dctsub(m, a, nc, w(nw))
if (m .gt. 4) then
call cftfsub(m, a, ip, nw, w)
call rftfsub(m, a, nc, w(nw))
else if (m .eq. 4) then
call cftfsub(m, a, ip, nw, w)
end if
a(n - 1) = a(0) - a(1)
a(1) = a(0) + a(1)
do j = m - 2, 2, -2
a(2 * j + 1) = a(j) + a(j + 1)
a(2 * j - 1) = a(j) - a(j + 1)
end do
l = 2
m = mh
do while (m .ge. 2)
call dctsub(m, t, nc, w(nw))
if (m .gt. 4) then
call cftfsub(m, t, ip, nw, w)
call rftfsub(m, t, nc, w(nw))
else if (m .eq. 4) then
call cftfsub(m, t, ip, nw, w)
end if
a(n - l) = t(0) - t(1)
a(l) = t(0) + t(1)
k = 0
do j = 2, m - 2, 2
k = k + 4 * l
a(k - l) = t(j) - t(j + 1)
a(k + l) = t(j) + t(j + 1)
end do
l = 2 * l
mh = m / 2
do j = 0, mh - 1
k = m - j
t(j) = t(m + k) - t(m + j)
t(k) = t(m + k) + t(m + j)
end do
t(mh) = t(m + mh)
m = mh
end do
a(l) = t(0)
a(n) = t(2) - t(1)
a(0) = t(2) + t(1)
else
a(1) = a(0)
a(2) = t(0)
a(0) = t(1)
end if
end
!
subroutine dfst(n, a, t, ip, w)
integer n, ip(0 : *), j, k, l, m, mh, nw, nc
real*8 a(0 : n - 1), t(0 : n / 2 - 1), w(0 : *), xr, xi, yr, yi
nw = ip(0)
if (n .gt. 8 * nw) then
nw = n / 8
call makewt(nw, ip, w)
end if
nc = ip(1)
if (n .gt. 2 * nc) then
nc = n / 2
call makect(nc, ip, w(nw))
end if
if (n .gt. 2) then
m = n / 2
mh = m / 2
do j = 1, mh - 1
k = m - j
xr = a(j) + a(n - j)
xi = a(j) - a(n - j)
yr = a(k) + a(n - k)
yi = a(k) - a(n - k)
a(j) = xr
a(k) = yr
t(j) = xi + yi
t(k) = xi - yi
end do
t(0) = a(mh) - a(n - mh)
a(mh) = a(mh) + a(n - mh)
a(0) = a(m)
call dstsub(m, a, nc, w(nw))
if (m .gt. 4) then
call cftfsub(m, a, ip, nw, w)
call rftfsub(m, a, nc, w(nw))
else if (m .eq. 4) then
call cftfsub(m, a, ip, nw, w)
end if
a(n - 1) = a(1) - a(0)
a(1) = a(0) + a(1)
do j = m - 2, 2, -2
a(2 * j + 1) = a(j) - a(j + 1)
a(2 * j - 1) = -a(j) - a(j + 1)
end do
l = 2
m = mh
do while (m .ge. 2)
call dstsub(m, t, nc, w(nw))
if (m .gt. 4) then
call cftfsub(m, t, ip, nw, w)
call rftfsub(m, t, nc, w(nw))
else if (m .eq. 4) then
call cftfsub(m, t, ip, nw, w)
end if
a(n - l) = t(1) - t(0)
a(l) = t(0) + t(1)
k = 0
do j = 2, m - 2, 2
k = k + 4 * l
a(k - l) = -t(j) - t(j + 1)
a(k + l) = t(j) - t(j + 1)
end do
l = 2 * l
mh = m / 2
do j = 1, mh - 1
k = m - j
t(j) = t(m + k) + t(m + j)
t(k) = t(m + k) - t(m + j)
end do
t(0) = t(m + mh)
m = mh
end do
a(l) = t(0)
end if
a(0) = 0
end
!
! -------- initializing routines --------
!
subroutine makewt(nw, ip, w)
integer nw, ip(0 : *), j, nwh, nw0, nw1
real*8 w(0 : nw - 1), delta, wn4r, wk1r, wk1i, wk3r, wk3i
ip(0) = nw
ip(1) = 1
if (nw .gt. 2) then
nwh = nw / 2
delta = atan(1.0d0) / nwh
wn4r = cos(delta * nwh)
w(0) = 1
w(1) = wn4r
if (nwh .eq. 4) then
w(2) = cos(delta * 2)
w(3) = sin(delta * 2)
else if (nwh .gt. 4) then
call makeipt(nw, ip)
w(2) = 0.5d0 / cos(delta * 2)
w(3) = 0.5d0 / cos(delta * 6)
do j = 4, nwh - 4, 4
w(j) = cos(delta * j)
w(j + 1) = sin(delta * j)
w(j + 2) = cos(3 * delta * j)
w(j + 3) = -sin(3 * delta * j)
end do
end if
nw0 = 0
do while (nwh .gt. 2)
nw1 = nw0 + nwh
nwh = nwh / 2
w(nw1) = 1
w(nw1 + 1) = wn4r
if (nwh .eq. 4) then
wk1r = w(nw0 + 4)
wk1i = w(nw0 + 5)
w(nw1 + 2) = wk1r
w(nw1 + 3) = wk1i
else if (nwh .gt. 4) then
wk1r = w(nw0 + 4)
wk3r = w(nw0 + 6)
w(nw1 + 2) = 0.5d0 / wk1r
w(nw1 + 3) = 0.5d0 / wk3r
do j = 4, nwh - 4, 4
wk1r = w(nw0 + 2 * j)
wk1i = w(nw0 + 2 * j + 1)
wk3r = w(nw0 + 2 * j + 2)
wk3i = w(nw0 + 2 * j + 3)
w(nw1 + j) = wk1r
w(nw1 + j + 1) = wk1i
w(nw1 + j + 2) = wk3r
w(nw1 + j + 3) = wk3i
end do
end if
nw0 = nw1
end do
end if
end
!
subroutine makeipt(nw, ip)
integer nw, ip(0 : *), j, l, m, m2, p, q
ip(2) = 0
ip(3) = 16
m = 2
l = nw
do while (l .gt. 32)
m2 = 2 * m
q = 8 * m2
do j = m, m2 - 1
p = 4 * ip(j)
ip(m + j) = p
ip(m2 + j) = p + q
end do
m = m2
l = l / 4
end do
end
!
subroutine makect(nc, ip, c)
integer nc, ip(0 : *), j, nch
real*8 c(0 : nc - 1), delta
ip(1) = nc
if (nc .gt. 1) then
nch = nc / 2
delta = atan(1.0d0) / nch
c(0) = cos(delta * nch)
c(nch) = 0.5d0 * c(0)
do j = 1, nch - 1
c(j) = 0.5d0 * cos(delta * j)
c(nc - j) = 0.5d0 * sin(delta * j)
end do
end if
end
!
! -------- child routines --------
!
subroutine cftfsub(n, a, ip, nw, w)
integer n, ip(0 : *), nw
real*8 a(0 : n - 1), w(0 : nw - 1)
if (n .gt. 8) then
if (n .gt. 32) then
call cftf1st(n, a, w(nw - n / 4))
if (n .gt. 512) then
call cftrec4(n, a, nw, w)
else if (n .gt. 128) then
call cftleaf(n, 1, a, nw, w)
else
call cftfx41(n, a, nw, w)
end if
call bitrv2(n, ip, a)
else if (n .eq. 32) then
call cftf161(a, w(nw - 8))
call bitrv216(a)
else
call cftf081(a, w)
call bitrv208(a)
end if
else if (n .eq. 8) then
call cftf040(a)
else if (n .eq. 4) then
call cftx020(a)
end if
end
!
subroutine cftbsub(n, a, ip, nw, w)
integer n, ip(0 : *), nw
real*8 a(0 : n - 1), w(0 : nw - 1)
if (n .gt. 8) then
if (n .gt. 32) then
call cftb1st(n, a, w(nw - n / 4))
if (n .gt. 512) then
call cftrec4(n, a, nw, w)
else if (n .gt. 128) then
call cftleaf(n, 1, a, nw, w)
else
call cftfx41(n, a, nw, w)
end if
call bitrv2conj(n, ip, a)
else if (n .eq. 32) then
call cftf161(a, w(nw - 8))
call bitrv216neg(a)
else
call cftf081(a, w)
call bitrv208neg(a)
end if
else if (n .eq. 8) then
call cftb040(a)
else if (n .eq. 4) then
call cftx020(a)
end if
end
!
subroutine bitrv2(n, ip, a)
integer n, ip(0 : *), j, j1, k, k1, l, m, nh, nm
real*8 a(0 : n - 1), xr, xi, yr, yi
m = 1
l = n / 4
do while (l .gt. 8)
m = m * 2
l = l / 4
end do
nh = n / 2
nm = 4 * m
if (l .eq. 8) then
do k = 0, m - 1
do j = 0, k - 1
j1 = 4 * j + 2 * ip(m + k)
k1 = 4 * k + 2 * ip(m + j)
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 - nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nh
k1 = k1 + 2
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 + nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + 2
k1 = k1 + nh
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 - nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nh
k1 = k1 - 2
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 + nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
end do
k1 = 4 * k + 2 * ip(m + k)
j1 = k1 + 2
k1 = k1 + nh
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 - nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - 2
k1 = k1 - nh
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nh + 2
k1 = k1 + nh + 2
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nh + nm
k1 = k1 + 2 * nm - 2
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
end do
else
do k = 0, m - 1
do j = 0, k - 1
j1 = 4 * j + ip(m + k)
k1 = 4 * k + ip(m + j)
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nh
k1 = k1 + 2
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + 2
k1 = k1 + nh
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nh
k1 = k1 - 2
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
end do
k1 = 4 * k + ip(m + k)
j1 = k1 + 2
k1 = k1 + nh
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + nm
xr = a(j1)
xi = a(j1 + 1)
yr = a(k1)
yi = a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
end do
end if
end
!
subroutine bitrv2conj(n, ip, a)
integer n, ip(0 : *), j, j1, k, k1, l, m, nh, nm
real*8 a(0 : n - 1), xr, xi, yr, yi
m = 1
l = n / 4
do while (l .gt. 8)
m = m * 2
l = l / 4
end do
nh = n / 2
nm = 4 * m
if (l .eq. 8) then
do k = 0, m - 1
do j = 0, k - 1
j1 = 4 * j + 2 * ip(m + k)
k1 = 4 * k + 2 * ip(m + j)
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 - nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nh
k1 = k1 + 2
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 + nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + 2
k1 = k1 + nh
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 - nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nh
k1 = k1 - 2
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 + nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
end do
k1 = 4 * k + 2 * ip(m + k)
j1 = k1 + 2
k1 = k1 + nh
a(j1 - 1) = -a(j1 - 1)
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
a(k1 + 3) = -a(k1 + 3)
j1 = j1 + nm
k1 = k1 + 2 * nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 - nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - 2
k1 = k1 - nh
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nh + 2
k1 = k1 + nh + 2
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nh + nm
k1 = k1 + 2 * nm - 2
a(j1 - 1) = -a(j1 - 1)
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
a(k1 + 3) = -a(k1 + 3)
end do
else
do k = 0, m - 1
do j = 0, k - 1
j1 = 4 * j + ip(m + k)
k1 = 4 * k + ip(m + j)
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nh
k1 = k1 + 2
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + 2
k1 = k1 + nh
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 + nm
k1 = k1 + nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nh
k1 = k1 - 2
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
j1 = j1 - nm
k1 = k1 - nm
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
end do
k1 = 4 * k + ip(m + k)
j1 = k1 + 2
k1 = k1 + nh
a(j1 - 1) = -a(j1 - 1)
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
a(k1 + 3) = -a(k1 + 3)
j1 = j1 + nm
k1 = k1 + nm
a(j1 - 1) = -a(j1 - 1)
xr = a(j1)
xi = -a(j1 + 1)
yr = a(k1)
yi = -a(k1 + 1)
a(j1) = yr
a(j1 + 1) = yi
a(k1) = xr
a(k1 + 1) = xi
a(k1 + 3) = -a(k1 + 3)
end do
end if
end
!
subroutine bitrv216(a)
real*8 a(0 : 31), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i
real*8 x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i
real*8 x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i
x1r = a(2)
x1i = a(3)
x2r = a(4)
x2i = a(5)
x3r = a(6)
x3i = a(7)
x4r = a(8)
x4i = a(9)
x5r = a(10)
x5i = a(11)
x7r = a(14)
x7i = a(15)
x8r = a(16)
x8i = a(17)
x10r = a(20)
x10i = a(21)
x11r = a(22)
x11i = a(23)
x12r = a(24)
x12i = a(25)
x13r = a(26)
x13i = a(27)
x14r = a(28)
x14i = a(29)
a(2) = x8r
a(3) = x8i
a(4) = x4r
a(5) = x4i
a(6) = x12r
a(7) = x12i
a(8) = x2r
a(9) = x2i
a(10) = x10r
a(11) = x10i
a(14) = x14r
a(15) = x14i
a(16) = x1r
a(17) = x1i
a(20) = x5r
a(21) = x5i
a(22) = x13r
a(23) = x13i
a(24) = x3r
a(25) = x3i
a(26) = x11r
a(27) = x11i
a(28) = x7r
a(29) = x7i
end
!
subroutine bitrv216neg(a)
real*8 a(0 : 31), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i
real*8 x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i
real*8 x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i
real*8 x13r, x13i, x14r, x14i, x15r, x15i
x1r = a(2)
x1i = a(3)
x2r = a(4)
x2i = a(5)
x3r = a(6)
x3i = a(7)
x4r = a(8)
x4i = a(9)
x5r = a(10)
x5i = a(11)
x6r = a(12)
x6i = a(13)
x7r = a(14)
x7i = a(15)
x8r = a(16)
x8i = a(17)
x9r = a(18)
x9i = a(19)
x10r = a(20)
x10i = a(21)
x11r = a(22)
x11i = a(23)
x12r = a(24)
x12i = a(25)
x13r = a(26)
x13i = a(27)
x14r = a(28)
x14i = a(29)
x15r = a(30)
x15i = a(31)
a(2) = x15r
a(3) = x15i
a(4) = x7r
a(5) = x7i
a(6) = x11r
a(7) = x11i
a(8) = x3r
a(9) = x3i
a(10) = x13r
a(11) = x13i
a(12) = x5r
a(13) = x5i
a(14) = x9r
a(15) = x9i
a(16) = x1r
a(17) = x1i
a(18) = x14r
a(19) = x14i
a(20) = x6r
a(21) = x6i
a(22) = x10r
a(23) = x10i
a(24) = x2r
a(25) = x2i
a(26) = x12r
a(27) = x12i
a(28) = x4r
a(29) = x4i
a(30) = x8r
a(31) = x8i
end
!
subroutine bitrv208(a)
real*8 a(0 : 15), x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i
x1r = a(2)
x1i = a(3)
x3r = a(6)
x3i = a(7)
x4r = a(8)
x4i = a(9)
x6r = a(12)
x6i = a(13)
a(2) = x4r
a(3) = x4i
a(6) = x6r
a(7) = x6i
a(8) = x1r
a(9) = x1i
a(12) = x3r
a(13) = x3i
end
!
subroutine bitrv208neg(a)
real*8 a(0 : 15), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i
real*8 x5r, x5i, x6r, x6i, x7r, x7i
x1r = a(2)
x1i = a(3)
x2r = a(4)
x2i = a(5)
x3r = a(6)
x3i = a(7)
x4r = a(8)
x4i = a(9)
x5r = a(10)
x5i = a(11)
x6r = a(12)
x6i = a(13)
x7r = a(14)
x7i = a(15)
a(2) = x7r
a(3) = x7i
a(4) = x3r
a(5) = x3i
a(6) = x5r
a(7) = x5i
a(8) = x1r
a(9) = x1i
a(10) = x6r
a(11) = x6i
a(12) = x2r
a(13) = x2i
a(14) = x4r
a(15) = x4i
end
!
subroutine cftf1st(n, a, w)
integer n, j, j0, j1, j2, j3, k, m, mh
real*8 a(0 : n - 1), w(0 : *)
real*8 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i
real*8 wd1r, wd1i, wd3r, wd3i
real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
mh = n / 8
m = 2 * mh
j1 = m
j2 = j1 + m
j3 = j2 + m
x0r = a(0) + a(j2)
x0i = a(1) + a(j2 + 1)
x1r = a(0) - a(j2)
x1i = a(1) - a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(0) = x0r + x2r
a(1) = x0i + x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
a(j2) = x1r - x3i
a(j2 + 1) = x1i + x3r
a(j3) = x1r + x3i
a(j3 + 1) = x1i - x3r
wn4r = w(1)
csc1 = w(2)
csc3 = w(3)
wd1r = 1
wd1i = 0
wd3r = 1
wd3i = 0
k = 0
do j = 2, mh - 6, 4
k = k + 4
wk1r = csc1 * (wd1r + w(k))
wk1i = csc1 * (wd1i + w(k + 1))
wk3r = csc3 * (wd3r + w(k + 2))
wk3i = csc3 * (wd3i + w(k + 3))
wd1r = w(k)
wd1i = w(k + 1)
wd3r = w(k + 2)
wd3i = w(k + 3)
j1 = j + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j) + a(j2)
x0i = a(j + 1) + a(j2 + 1)
x1r = a(j) - a(j2)
x1i = a(j + 1) - a(j2 + 1)
y0r = a(j + 2) + a(j2 + 2)
y0i = a(j + 3) + a(j2 + 3)
y1r = a(j + 2) - a(j2 + 2)
y1i = a(j + 3) - a(j2 + 3)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
y2r = a(j1 + 2) + a(j3 + 2)
y2i = a(j1 + 3) + a(j3 + 3)
y3r = a(j1 + 2) - a(j3 + 2)
y3i = a(j1 + 3) - a(j3 + 3)
a(j) = x0r + x2r
a(j + 1) = x0i + x2i
a(j + 2) = y0r + y2r
a(j + 3) = y0i + y2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
a(j1 + 2) = y0r - y2r
a(j1 + 3) = y0i - y2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2) = wk1r * x0r - wk1i * x0i
a(j2 + 1) = wk1r * x0i + wk1i * x0r
x0r = y1r - y3i
x0i = y1i + y3r
a(j2 + 2) = wd1r * x0r - wd1i * x0i
a(j2 + 3) = wd1r * x0i + wd1i * x0r
x0r = x1r + x3i
x0i = x1i - x3r
a(j3) = wk3r * x0r + wk3i * x0i
a(j3 + 1) = wk3r * x0i - wk3i * x0r
x0r = y1r + y3i
x0i = y1i - y3r
a(j3 + 2) = wd3r * x0r + wd3i * x0i
a(j3 + 3) = wd3r * x0i - wd3i * x0r
j0 = m - j
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0) + a(j2)
x0i = a(j0 + 1) + a(j2 + 1)
x1r = a(j0) - a(j2)
x1i = a(j0 + 1) - a(j2 + 1)
y0r = a(j0 - 2) + a(j2 - 2)
y0i = a(j0 - 1) + a(j2 - 1)
y1r = a(j0 - 2) - a(j2 - 2)
y1i = a(j0 - 1) - a(j2 - 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
y2r = a(j1 - 2) + a(j3 - 2)
y2i = a(j1 - 1) + a(j3 - 1)
y3r = a(j1 - 2) - a(j3 - 2)
y3i = a(j1 - 1) - a(j3 - 1)
a(j0) = x0r + x2r
a(j0 + 1) = x0i + x2i
a(j0 - 2) = y0r + y2r
a(j0 - 1) = y0i + y2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
a(j1 - 2) = y0r - y2r
a(j1 - 1) = y0i - y2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2) = wk1i * x0r - wk1r * x0i
a(j2 + 1) = wk1i * x0i + wk1r * x0r
x0r = y1r - y3i
x0i = y1i + y3r
a(j2 - 2) = wd1i * x0r - wd1r * x0i
a(j2 - 1) = wd1i * x0i + wd1r * x0r
x0r = x1r + x3i
x0i = x1i - x3r
a(j3) = wk3i * x0r + wk3r * x0i
a(j3 + 1) = wk3i * x0i - wk3r * x0r
x0r = y1r + y3i
x0i = y1i - y3r
a(j3 - 2) = wd3i * x0r + wd3r * x0i
a(j3 - 1) = wd3i * x0i - wd3r * x0r
end do
wk1r = csc1 * (wd1r + wn4r)
wk1i = csc1 * (wd1i + wn4r)
wk3r = csc3 * (wd3r - wn4r)
wk3i = csc3 * (wd3i - wn4r)
j0 = mh
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0 - 2) + a(j2 - 2)
x0i = a(j0 - 1) + a(j2 - 1)
x1r = a(j0 - 2) - a(j2 - 2)
x1i = a(j0 - 1) - a(j2 - 1)
x2r = a(j1 - 2) + a(j3 - 2)
x2i = a(j1 - 1) + a(j3 - 1)
x3r = a(j1 - 2) - a(j3 - 2)
x3i = a(j1 - 1) - a(j3 - 1)
a(j0 - 2) = x0r + x2r
a(j0 - 1) = x0i + x2i
a(j1 - 2) = x0r - x2r
a(j1 - 1) = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2 - 2) = wk1r * x0r - wk1i * x0i
a(j2 - 1) = wk1r * x0i + wk1i * x0r
x0r = x1r + x3i
x0i = x1i - x3r
a(j3 - 2) = wk3r * x0r + wk3i * x0i
a(j3 - 1) = wk3r * x0i - wk3i * x0r
x0r = a(j0) + a(j2)
x0i = a(j0 + 1) + a(j2 + 1)
x1r = a(j0) - a(j2)
x1i = a(j0 + 1) - a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(j0) = x0r + x2r
a(j0 + 1) = x0i + x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2) = wn4r * (x0r - x0i)
a(j2 + 1) = wn4r * (x0i + x0r)
x0r = x1r + x3i
x0i = x1i - x3r
a(j3) = -wn4r * (x0r + x0i)
a(j3 + 1) = -wn4r * (x0i - x0r)
x0r = a(j0 + 2) + a(j2 + 2)
x0i = a(j0 + 3) + a(j2 + 3)
x1r = a(j0 + 2) - a(j2 + 2)
x1i = a(j0 + 3) - a(j2 + 3)
x2r = a(j1 + 2) + a(j3 + 2)
x2i = a(j1 + 3) + a(j3 + 3)
x3r = a(j1 + 2) - a(j3 + 2)
x3i = a(j1 + 3) - a(j3 + 3)
a(j0 + 2) = x0r + x2r
a(j0 + 3) = x0i + x2i
a(j1 + 2) = x0r - x2r
a(j1 + 3) = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2 + 2) = wk1i * x0r - wk1r * x0i
a(j2 + 3) = wk1i * x0i + wk1r * x0r
x0r = x1r + x3i
x0i = x1i - x3r
a(j3 + 2) = wk3i * x0r + wk3r * x0i
a(j3 + 3) = wk3i * x0i - wk3r * x0r
end
!
subroutine cftb1st(n, a, w)
integer n, j, j0, j1, j2, j3, k, m, mh
real*8 a(0 : n - 1), w(0 : *)
real*8 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i
real*8 wd1r, wd1i, wd3r, wd3i
real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
mh = n / 8
m = 2 * mh
j1 = m
j2 = j1 + m
j3 = j2 + m
x0r = a(0) + a(j2)
x0i = -a(1) - a(j2 + 1)
x1r = a(0) - a(j2)
x1i = -a(1) + a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(0) = x0r + x2r
a(1) = x0i - x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i + x2i
a(j2) = x1r + x3i
a(j2 + 1) = x1i + x3r
a(j3) = x1r - x3i
a(j3 + 1) = x1i - x3r
wn4r = w(1)
csc1 = w(2)
csc3 = w(3)
wd1r = 1
wd1i = 0
wd3r = 1
wd3i = 0
k = 0
do j = 2, mh - 6, 4
k = k + 4
wk1r = csc1 * (wd1r + w(k))
wk1i = csc1 * (wd1i + w(k + 1))
wk3r = csc3 * (wd3r + w(k + 2))
wk3i = csc3 * (wd3i + w(k + 3))
wd1r = w(k)
wd1i = w(k + 1)
wd3r = w(k + 2)
wd3i = w(k + 3)
j1 = j + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j) + a(j2)
x0i = -a(j + 1) - a(j2 + 1)
x1r = a(j) - a(j2)
x1i = -a(j + 1) + a(j2 + 1)
y0r = a(j + 2) + a(j2 + 2)
y0i = -a(j + 3) - a(j2 + 3)
y1r = a(j + 2) - a(j2 + 2)
y1i = -a(j + 3) + a(j2 + 3)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
y2r = a(j1 + 2) + a(j3 + 2)
y2i = a(j1 + 3) + a(j3 + 3)
y3r = a(j1 + 2) - a(j3 + 2)
y3i = a(j1 + 3) - a(j3 + 3)
a(j) = x0r + x2r
a(j + 1) = x0i - x2i
a(j + 2) = y0r + y2r
a(j + 3) = y0i - y2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i + x2i
a(j1 + 2) = y0r - y2r
a(j1 + 3) = y0i + y2i
x0r = x1r + x3i
x0i = x1i + x3r
a(j2) = wk1r * x0r - wk1i * x0i
a(j2 + 1) = wk1r * x0i + wk1i * x0r
x0r = y1r + y3i
x0i = y1i + y3r
a(j2 + 2) = wd1r * x0r - wd1i * x0i
a(j2 + 3) = wd1r * x0i + wd1i * x0r
x0r = x1r - x3i
x0i = x1i - x3r
a(j3) = wk3r * x0r + wk3i * x0i
a(j3 + 1) = wk3r * x0i - wk3i * x0r
x0r = y1r - y3i
x0i = y1i - y3r
a(j3 + 2) = wd3r * x0r + wd3i * x0i
a(j3 + 3) = wd3r * x0i - wd3i * x0r
j0 = m - j
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0) + a(j2)
x0i = -a(j0 + 1) - a(j2 + 1)
x1r = a(j0) - a(j2)
x1i = -a(j0 + 1) + a(j2 + 1)
y0r = a(j0 - 2) + a(j2 - 2)
y0i = -a(j0 - 1) - a(j2 - 1)
y1r = a(j0 - 2) - a(j2 - 2)
y1i = -a(j0 - 1) + a(j2 - 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
y2r = a(j1 - 2) + a(j3 - 2)
y2i = a(j1 - 1) + a(j3 - 1)
y3r = a(j1 - 2) - a(j3 - 2)
y3i = a(j1 - 1) - a(j3 - 1)
a(j0) = x0r + x2r
a(j0 + 1) = x0i - x2i
a(j0 - 2) = y0r + y2r
a(j0 - 1) = y0i - y2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i + x2i
a(j1 - 2) = y0r - y2r
a(j1 - 1) = y0i + y2i
x0r = x1r + x3i
x0i = x1i + x3r
a(j2) = wk1i * x0r - wk1r * x0i
a(j2 + 1) = wk1i * x0i + wk1r * x0r
x0r = y1r + y3i
x0i = y1i + y3r
a(j2 - 2) = wd1i * x0r - wd1r * x0i
a(j2 - 1) = wd1i * x0i + wd1r * x0r
x0r = x1r - x3i
x0i = x1i - x3r
a(j3) = wk3i * x0r + wk3r * x0i
a(j3 + 1) = wk3i * x0i - wk3r * x0r
x0r = y1r - y3i
x0i = y1i - y3r
a(j3 - 2) = wd3i * x0r + wd3r * x0i
a(j3 - 1) = wd3i * x0i - wd3r * x0r
end do
wk1r = csc1 * (wd1r + wn4r)
wk1i = csc1 * (wd1i + wn4r)
wk3r = csc3 * (wd3r - wn4r)
wk3i = csc3 * (wd3i - wn4r)
j0 = mh
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0 - 2) + a(j2 - 2)
x0i = -a(j0 - 1) - a(j2 - 1)
x1r = a(j0 - 2) - a(j2 - 2)
x1i = -a(j0 - 1) + a(j2 - 1)
x2r = a(j1 - 2) + a(j3 - 2)
x2i = a(j1 - 1) + a(j3 - 1)
x3r = a(j1 - 2) - a(j3 - 2)
x3i = a(j1 - 1) - a(j3 - 1)
a(j0 - 2) = x0r + x2r
a(j0 - 1) = x0i - x2i
a(j1 - 2) = x0r - x2r
a(j1 - 1) = x0i + x2i
x0r = x1r + x3i
x0i = x1i + x3r
a(j2 - 2) = wk1r * x0r - wk1i * x0i
a(j2 - 1) = wk1r * x0i + wk1i * x0r
x0r = x1r - x3i
x0i = x1i - x3r
a(j3 - 2) = wk3r * x0r + wk3i * x0i
a(j3 - 1) = wk3r * x0i - wk3i * x0r
x0r = a(j0) + a(j2)
x0i = -a(j0 + 1) - a(j2 + 1)
x1r = a(j0) - a(j2)
x1i = -a(j0 + 1) + a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(j0) = x0r + x2r
a(j0 + 1) = x0i - x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i + x2i
x0r = x1r + x3i
x0i = x1i + x3r
a(j2) = wn4r * (x0r - x0i)
a(j2 + 1) = wn4r * (x0i + x0r)
x0r = x1r - x3i
x0i = x1i - x3r
a(j3) = -wn4r * (x0r + x0i)
a(j3 + 1) = -wn4r * (x0i - x0r)
x0r = a(j0 + 2) + a(j2 + 2)
x0i = -a(j0 + 3) - a(j2 + 3)
x1r = a(j0 + 2) - a(j2 + 2)
x1i = -a(j0 + 3) + a(j2 + 3)
x2r = a(j1 + 2) + a(j3 + 2)
x2i = a(j1 + 3) + a(j3 + 3)
x3r = a(j1 + 2) - a(j3 + 2)
x3i = a(j1 + 3) - a(j3 + 3)
a(j0 + 2) = x0r + x2r
a(j0 + 3) = x0i - x2i
a(j1 + 2) = x0r - x2r
a(j1 + 3) = x0i + x2i
x0r = x1r + x3i
x0i = x1i + x3r
a(j2 + 2) = wk1i * x0r - wk1r * x0i
a(j2 + 3) = wk1i * x0i + wk1r * x0r
x0r = x1r - x3i
x0i = x1i - x3r
a(j3 + 2) = wk3i * x0r + wk3r * x0i
a(j3 + 3) = wk3i * x0i - wk3r * x0r
end
!
subroutine cftrec4(n, a, nw, w)
integer n, nw, cfttree, isplt, j, k, m
real*8 a(0 : n - 1), w(0 : nw - 1)
m = n
do while (m .gt. 512)
m = m / 4
call cftmdl1(m, a(n - m), w(nw - m / 2))
end do
call cftleaf(m, 1, a(n - m), nw, w)
k = 0
do j = n - m, m, -m
k = k + 1
isplt = cfttree(m, j, k, a, nw, w)
call cftleaf(m, isplt, a(j - m), nw, w)
end do
end
!
integer function cfttree(n, j, k, a, nw, w)
integer n, j, k, nw, i, isplt, m
real*8 a(0 : j - 1), w(0 : nw - 1)
if (mod(k, 4) .ne. 0) then
isplt = mod(k, 2)
if (isplt .ne. 0) then
call cftmdl1(n, a(j - n), w(nw - n / 2))
else
call cftmdl2(n, a(j - n), w(nw - n))
end if
else
m = n
i = k
do while (mod(i, 4) .eq. 0)
m = m * 4
i = i / 4
end do
isplt = mod(i, 2)
if (isplt .ne. 0) then
do while (m .gt. 128)
call cftmdl1(m, a(j - m), w(nw - m / 2))
m = m / 4
end do
else
do while (m .gt. 128)
call cftmdl2(m, a(j - m), w(nw - m))
m = m / 4
end do
end if
end if
cfttree = isplt
end
!
subroutine cftleaf(n, isplt, a, nw, w)
integer n, isplt, nw
real*8 a(0 : n - 1), w(0 : nw - 1)
if (n .eq. 512) then
call cftmdl1(128, a, w(nw - 64))
call cftf161(a, w(nw - 8))
call cftf162(a(32), w(nw - 32))
call cftf161(a(64), w(nw - 8))
call cftf161(a(96), w(nw - 8))
call cftmdl2(128, a(128), w(nw - 128))
call cftf161(a(128), w(nw - 8))
call cftf162(a(160), w(nw - 32))
call cftf161(a(192), w(nw - 8))
call cftf162(a(224), w(nw - 32))
call cftmdl1(128, a(256), w(nw - 64))
call cftf161(a(256), w(nw - 8))
call cftf162(a(288), w(nw - 32))
call cftf161(a(320), w(nw - 8))
call cftf161(a(352), w(nw - 8))
if (isplt .ne. 0) then
call cftmdl1(128, a(384), w(nw - 64))
call cftf161(a(480), w(nw - 8))
else
call cftmdl2(128, a(384), w(nw - 128))
call cftf162(a(480), w(nw - 32))
end if
call cftf161(a(384), w(nw - 8))
call cftf162(a(416), w(nw - 32))
call cftf161(a(448), w(nw - 8))
else
call cftmdl1(64, a, w(nw - 32))
call cftf081(a, w(nw - 8))
call cftf082(a(16), w(nw - 8))
call cftf081(a(32), w(nw - 8))
call cftf081(a(48), w(nw - 8))
call cftmdl2(64, a(64), w(nw - 64))
call cftf081(a(64), w(nw - 8))
call cftf082(a(80), w(nw - 8))
call cftf081(a(96), w(nw - 8))
call cftf082(a(112), w(nw - 8))
call cftmdl1(64, a(128), w(nw - 32))
call cftf081(a(128), w(nw - 8))
call cftf082(a(144), w(nw - 8))
call cftf081(a(160), w(nw - 8))
call cftf081(a(176), w(nw - 8))
if (isplt .ne. 0) then
call cftmdl1(64, a(192), w(nw - 32))
call cftf081(a(240), w(nw - 8))
else
call cftmdl2(64, a(192), w(nw - 64))
call cftf082(a(240), w(nw - 8))
end if
call cftf081(a(192), w(nw - 8))
call cftf082(a(208), w(nw - 8))
call cftf081(a(224), w(nw - 8))
end if
end
!
subroutine cftmdl1(n, a, w)
integer n, j, j0, j1, j2, j3, k, m, mh
real*8 a(0 : n - 1), w(0 : *)
real*8 wn4r, wk1r, wk1i, wk3r, wk3i
real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
mh = n / 8
m = 2 * mh
j1 = m
j2 = j1 + m
j3 = j2 + m
x0r = a(0) + a(j2)
x0i = a(1) + a(j2 + 1)
x1r = a(0) - a(j2)
x1i = a(1) - a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(0) = x0r + x2r
a(1) = x0i + x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
a(j2) = x1r - x3i
a(j2 + 1) = x1i + x3r
a(j3) = x1r + x3i
a(j3 + 1) = x1i - x3r
wn4r = w(1)
k = 0
do j = 2, mh - 2, 2
k = k + 4
wk1r = w(k)
wk1i = w(k + 1)
wk3r = w(k + 2)
wk3i = w(k + 3)
j1 = j + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j) + a(j2)
x0i = a(j + 1) + a(j2 + 1)
x1r = a(j) - a(j2)
x1i = a(j + 1) - a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(j) = x0r + x2r
a(j + 1) = x0i + x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2) = wk1r * x0r - wk1i * x0i
a(j2 + 1) = wk1r * x0i + wk1i * x0r
x0r = x1r + x3i
x0i = x1i - x3r
a(j3) = wk3r * x0r + wk3i * x0i
a(j3 + 1) = wk3r * x0i - wk3i * x0r
j0 = m - j
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0) + a(j2)
x0i = a(j0 + 1) + a(j2 + 1)
x1r = a(j0) - a(j2)
x1i = a(j0 + 1) - a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(j0) = x0r + x2r
a(j0 + 1) = x0i + x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2) = wk1i * x0r - wk1r * x0i
a(j2 + 1) = wk1i * x0i + wk1r * x0r
x0r = x1r + x3i
x0i = x1i - x3r
a(j3) = wk3i * x0r + wk3r * x0i
a(j3 + 1) = wk3i * x0i - wk3r * x0r
end do
j0 = mh
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0) + a(j2)
x0i = a(j0 + 1) + a(j2 + 1)
x1r = a(j0) - a(j2)
x1i = a(j0 + 1) - a(j2 + 1)
x2r = a(j1) + a(j3)
x2i = a(j1 + 1) + a(j3 + 1)
x3r = a(j1) - a(j3)
x3i = a(j1 + 1) - a(j3 + 1)
a(j0) = x0r + x2r
a(j0 + 1) = x0i + x2i
a(j1) = x0r - x2r
a(j1 + 1) = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
a(j2) = wn4r * (x0r - x0i)
a(j2 + 1) = wn4r * (x0i + x0r)
x0r = x1r + x3i
x0i = x1i - x3r
a(j3) = -wn4r * (x0r + x0i)
a(j3 + 1) = -wn4r * (x0i - x0r)
end
!
subroutine cftmdl2(n, a, w)
integer n, j, j0, j1, j2, j3, k, kr, m, mh
real*8 a(0 : n - 1), w(0 : *)
real*8 wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i
real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
real*8 y0r, y0i, y2r, y2i
mh = n / 8
m = 2 * mh
wn4r = w(1)
j1 = m
j2 = j1 + m
j3 = j2 + m
x0r = a(0) - a(j2 + 1)
x0i = a(1) + a(j2)
x1r = a(0) + a(j2 + 1)
x1i = a(1) - a(j2)
x2r = a(j1) - a(j3 + 1)
x2i = a(j1 + 1) + a(j3)
x3r = a(j1) + a(j3 + 1)
x3i = a(j1 + 1) - a(j3)
y0r = wn4r * (x2r - x2i)
y0i = wn4r * (x2i + x2r)
a(0) = x0r + y0r
a(1) = x0i + y0i
a(j1) = x0r - y0r
a(j1 + 1) = x0i - y0i
y0r = wn4r * (x3r - x3i)
y0i = wn4r * (x3i + x3r)
a(j2) = x1r - y0i
a(j2 + 1) = x1i + y0r
a(j3) = x1r + y0i
a(j3 + 1) = x1i - y0r
k = 0
kr = 2 * m
do j = 2, mh - 2, 2
k = k + 4
wk1r = w(k)
wk1i = w(k + 1)
wk3r = w(k + 2)
wk3i = w(k + 3)
kr = kr - 4
wd1i = w(kr)
wd1r = w(kr + 1)
wd3i = w(kr + 2)
wd3r = w(kr + 3)
j1 = j + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j) - a(j2 + 1)
x0i = a(j + 1) + a(j2)
x1r = a(j) + a(j2 + 1)
x1i = a(j + 1) - a(j2)
x2r = a(j1) - a(j3 + 1)
x2i = a(j1 + 1) + a(j3)
x3r = a(j1) + a(j3 + 1)
x3i = a(j1 + 1) - a(j3)
y0r = wk1r * x0r - wk1i * x0i
y0i = wk1r * x0i + wk1i * x0r
y2r = wd1r * x2r - wd1i * x2i
y2i = wd1r * x2i + wd1i * x2r
a(j) = y0r + y2r
a(j + 1) = y0i + y2i
a(j1) = y0r - y2r
a(j1 + 1) = y0i - y2i
y0r = wk3r * x1r + wk3i * x1i
y0i = wk3r * x1i - wk3i * x1r
y2r = wd3r * x3r + wd3i * x3i
y2i = wd3r * x3i - wd3i * x3r
a(j2) = y0r + y2r
a(j2 + 1) = y0i + y2i
a(j3) = y0r - y2r
a(j3 + 1) = y0i - y2i
j0 = m - j
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0) - a(j2 + 1)
x0i = a(j0 + 1) + a(j2)
x1r = a(j0) + a(j2 + 1)
x1i = a(j0 + 1) - a(j2)
x2r = a(j1) - a(j3 + 1)
x2i = a(j1 + 1) + a(j3)
x3r = a(j1) + a(j3 + 1)
x3i = a(j1 + 1) - a(j3)
y0r = wd1i * x0r - wd1r * x0i
y0i = wd1i * x0i + wd1r * x0r
y2r = wk1i * x2r - wk1r * x2i
y2i = wk1i * x2i + wk1r * x2r
a(j0) = y0r + y2r
a(j0 + 1) = y0i + y2i
a(j1) = y0r - y2r
a(j1 + 1) = y0i - y2i
y0r = wd3i * x1r + wd3r * x1i
y0i = wd3i * x1i - wd3r * x1r
y2r = wk3i * x3r + wk3r * x3i
y2i = wk3i * x3i - wk3r * x3r
a(j2) = y0r + y2r
a(j2 + 1) = y0i + y2i
a(j3) = y0r - y2r
a(j3 + 1) = y0i - y2i
end do
wk1r = w(m)
wk1i = w(m + 1)
j0 = mh
j1 = j0 + m
j2 = j1 + m
j3 = j2 + m
x0r = a(j0) - a(j2 + 1)
x0i = a(j0 + 1) + a(j2)
x1r = a(j0) + a(j2 + 1)
x1i = a(j0 + 1) - a(j2)
x2r = a(j1) - a(j3 + 1)
x2i = a(j1 + 1) + a(j3)
x3r = a(j1) + a(j3 + 1)
x3i = a(j1 + 1) - a(j3)
y0r = wk1r * x0r - wk1i * x0i
y0i = wk1r * x0i + wk1i * x0r
y2r = wk1i * x2r - wk1r * x2i
y2i = wk1i * x2i + wk1r * x2r
a(j0) = y0r + y2r
a(j0 + 1) = y0i + y2i
a(j1) = y0r - y2r
a(j1 + 1) = y0i - y2i
y0r = wk1i * x1r - wk1r * x1i
y0i = wk1i * x1i + wk1r * x1r
y2r = wk1r * x3r - wk1i * x3i
y2i = wk1r * x3i + wk1i * x3r
a(j2) = y0r - y2r
a(j2 + 1) = y0i - y2i
a(j3) = y0r + y2r
a(j3 + 1) = y0i + y2i
end
!
subroutine cftfx41(n, a, nw, w)
integer n, nw
real*8 a(0 : n - 1), w(0 : nw - 1)
if (n .eq. 128) then
call cftf161(a, w(nw - 8))
call cftf162(a(32), w(nw - 32))
call cftf161(a(64), w(nw - 8))
call cftf161(a(96), w(nw - 8))
else
call cftf081(a, w(nw - 8))
call cftf082(a(16), w(nw - 8))
call cftf081(a(32), w(nw - 8))
call cftf081(a(48), w(nw - 8))
end if
end
!
subroutine cftf161(a, w)
real*8 a(0 : 31), w(0 : *), wn4r, wk1r, wk1i
real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
real*8 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i
real*8 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i
wn4r = w(1)
wk1r = w(2)
wk1i = w(3)
x0r = a(0) + a(16)
x0i = a(1) + a(17)
x1r = a(0) - a(16)
x1i = a(1) - a(17)
x2r = a(8) + a(24)
x2i = a(9) + a(25)
x3r = a(8) - a(24)
x3i = a(9) - a(25)
y0r = x0r + x2r
y0i = x0i + x2i
y4r = x0r - x2r
y4i = x0i - x2i
y8r = x1r - x3i
y8i = x1i + x3r
y12r = x1r + x3i
y12i = x1i - x3r
x0r = a(2) + a(18)
x0i = a(3) + a(19)
x1r = a(2) - a(18)
x1i = a(3) - a(19)
x2r = a(10) + a(26)
x2i = a(11) + a(27)
x3r = a(10) - a(26)
x3i = a(11) - a(27)
y1r = x0r + x2r
y1i = x0i + x2i
y5r = x0r - x2r
y5i = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
y9r = wk1r * x0r - wk1i * x0i
y9i = wk1r * x0i + wk1i * x0r
x0r = x1r + x3i
x0i = x1i - x3r
y13r = wk1i * x0r - wk1r * x0i
y13i = wk1i * x0i + wk1r * x0r
x0r = a(4) + a(20)
x0i = a(5) + a(21)
x1r = a(4) - a(20)
x1i = a(5) - a(21)
x2r = a(12) + a(28)
x2i = a(13) + a(29)
x3r = a(12) - a(28)
x3i = a(13) - a(29)
y2r = x0r + x2r
y2i = x0i + x2i
y6r = x0r - x2r
y6i = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
y10r = wn4r * (x0r - x0i)
y10i = wn4r * (x0i + x0r)
x0r = x1r + x3i
x0i = x1i - x3r
y14r = wn4r * (x0r + x0i)
y14i = wn4r * (x0i - x0r)
x0r = a(6) + a(22)
x0i = a(7) + a(23)
x1r = a(6) - a(22)
x1i = a(7) - a(23)
x2r = a(14) + a(30)
x2i = a(15) + a(31)
x3r = a(14) - a(30)
x3i = a(15) - a(31)
y3r = x0r + x2r
y3i = x0i + x2i
y7r = x0r - x2r
y7i = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
y11r = wk1i * x0r - wk1r * x0i
y11i = wk1i * x0i + wk1r * x0r
x0r = x1r + x3i
x0i = x1i - x3r
y15r = wk1r * x0r - wk1i * x0i
y15i = wk1r * x0i + wk1i * x0r
x0r = y12r - y14r
x0i = y12i - y14i
x1r = y12r + y14r
x1i = y12i + y14i
x2r = y13r - y15r
x2i = y13i - y15i
x3r = y13r + y15r
x3i = y13i + y15i
a(24) = x0r + x2r
a(25) = x0i + x2i
a(26) = x0r - x2r
a(27) = x0i - x2i
a(28) = x1r - x3i
a(29) = x1i + x3r
a(30) = x1r + x3i
a(31) = x1i - x3r
x0r = y8r + y10r
x0i = y8i + y10i
x1r = y8r - y10r
x1i = y8i - y10i
x2r = y9r + y11r
x2i = y9i + y11i
x3r = y9r - y11r
x3i = y9i - y11i
a(16) = x0r + x2r
a(17) = x0i + x2i
a(18) = x0r - x2r
a(19) = x0i - x2i
a(20) = x1r - x3i
a(21) = x1i + x3r
a(22) = x1r + x3i
a(23) = x1i - x3r
x0r = y5r - y7i
x0i = y5i + y7r
x2r = wn4r * (x0r - x0i)
x2i = wn4r * (x0i + x0r)
x0r = y5r + y7i
x0i = y5i - y7r
x3r = wn4r * (x0r - x0i)
x3i = wn4r * (x0i + x0r)
x0r = y4r - y6i
x0i = y4i + y6r
x1r = y4r + y6i
x1i = y4i - y6r
a(8) = x0r + x2r
a(9) = x0i + x2i
a(10) = x0r - x2r
a(11) = x0i - x2i
a(12) = x1r - x3i
a(13) = x1i + x3r
a(14) = x1r + x3i
a(15) = x1i - x3r
x0r = y0r + y2r
x0i = y0i + y2i
x1r = y0r - y2r
x1i = y0i - y2i
x2r = y1r + y3r
x2i = y1i + y3i
x3r = y1r - y3r
x3i = y1i - y3i
a(0) = x0r + x2r
a(1) = x0i + x2i
a(2) = x0r - x2r
a(3) = x0i - x2i
a(4) = x1r - x3i
a(5) = x1i + x3r
a(6) = x1r + x3i
a(7) = x1i - x3r
end
!
subroutine cftf162(a, w)
real*8 a(0 : 31), w(0 : *)
real*8 wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i
real*8 x0r, x0i, x1r, x1i, x2r, x2i
real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
real*8 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i
real*8 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i
wn4r = w(1)
wk1r = w(4)
wk1i = w(5)
wk3r = w(6)
wk3i = -w(7)
wk2r = w(8)
wk2i = w(9)
x1r = a(0) - a(17)
x1i = a(1) + a(16)
x0r = a(8) - a(25)
x0i = a(9) + a(24)
x2r = wn4r * (x0r - x0i)
x2i = wn4r * (x0i + x0r)
y0r = x1r + x2r
y0i = x1i + x2i
y4r = x1r - x2r
y4i = x1i - x2i
x1r = a(0) + a(17)
x1i = a(1) - a(16)
x0r = a(8) + a(25)
x0i = a(9) - a(24)
x2r = wn4r * (x0r - x0i)
x2i = wn4r * (x0i + x0r)
y8r = x1r - x2i
y8i = x1i + x2r
y12r = x1r + x2i
y12i = x1i - x2r
x0r = a(2) - a(19)
x0i = a(3) + a(18)
x1r = wk1r * x0r - wk1i * x0i
x1i = wk1r * x0i + wk1i * x0r
x0r = a(10) - a(27)
x0i = a(11) + a(26)
x2r = wk3i * x0r - wk3r * x0i
x2i = wk3i * x0i + wk3r * x0r
y1r = x1r + x2r
y1i = x1i + x2i
y5r = x1r - x2r
y5i = x1i - x2i
x0r = a(2) + a(19)
x0i = a(3) - a(18)
x1r = wk3r * x0r - wk3i * x0i
x1i = wk3r * x0i + wk3i * x0r
x0r = a(10) + a(27)
x0i = a(11) - a(26)
x2r = wk1r * x0r + wk1i * x0i
x2i = wk1r * x0i - wk1i * x0r
y9r = x1r - x2r
y9i = x1i - x2i
y13r = x1r + x2r
y13i = x1i + x2i
x0r = a(4) - a(21)
x0i = a(5) + a(20)
x1r = wk2r * x0r - wk2i * x0i
x1i = wk2r * x0i + wk2i * x0r
x0r = a(12) - a(29)
x0i = a(13) + a(28)
x2r = wk2i * x0r - wk2r * x0i
x2i = wk2i * x0i + wk2r * x0r
y2r = x1r + x2r
y2i = x1i + x2i
y6r = x1r - x2r
y6i = x1i - x2i
x0r = a(4) + a(21)
x0i = a(5) - a(20)
x1r = wk2i * x0r - wk2r * x0i
x1i = wk2i * x0i + wk2r * x0r
x0r = a(12) + a(29)
x0i = a(13) - a(28)
x2r = wk2r * x0r - wk2i * x0i
x2i = wk2r * x0i + wk2i * x0r
y10r = x1r - x2r
y10i = x1i - x2i
y14r = x1r + x2r
y14i = x1i + x2i
x0r = a(6) - a(23)
x0i = a(7) + a(22)
x1r = wk3r * x0r - wk3i * x0i
x1i = wk3r * x0i + wk3i * x0r
x0r = a(14) - a(31)
x0i = a(15) + a(30)
x2r = wk1i * x0r - wk1r * x0i
x2i = wk1i * x0i + wk1r * x0r
y3r = x1r + x2r
y3i = x1i + x2i
y7r = x1r - x2r
y7i = x1i - x2i
x0r = a(6) + a(23)
x0i = a(7) - a(22)
x1r = wk1i * x0r + wk1r * x0i
x1i = wk1i * x0i - wk1r * x0r
x0r = a(14) + a(31)
x0i = a(15) - a(30)
x2r = wk3i * x0r - wk3r * x0i
x2i = wk3i * x0i + wk3r * x0r
y11r = x1r + x2r
y11i = x1i + x2i
y15r = x1r - x2r
y15i = x1i - x2i
x1r = y0r + y2r
x1i = y0i + y2i
x2r = y1r + y3r
x2i = y1i + y3i
a(0) = x1r + x2r
a(1) = x1i + x2i
a(2) = x1r - x2r
a(3) = x1i - x2i
x1r = y0r - y2r
x1i = y0i - y2i
x2r = y1r - y3r
x2i = y1i - y3i
a(4) = x1r - x2i
a(5) = x1i + x2r
a(6) = x1r + x2i
a(7) = x1i - x2r
x1r = y4r - y6i
x1i = y4i + y6r
x0r = y5r - y7i
x0i = y5i + y7r
x2r = wn4r * (x0r - x0i)
x2i = wn4r * (x0i + x0r)
a(8) = x1r + x2r
a(9) = x1i + x2i
a(10) = x1r - x2r
a(11) = x1i - x2i
x1r = y4r + y6i
x1i = y4i - y6r
x0r = y5r + y7i
x0i = y5i - y7r
x2r = wn4r * (x0r - x0i)
x2i = wn4r * (x0i + x0r)
a(12) = x1r - x2i
a(13) = x1i + x2r
a(14) = x1r + x2i
a(15) = x1i - x2r
x1r = y8r + y10r
x1i = y8i + y10i
x2r = y9r - y11r
x2i = y9i - y11i
a(16) = x1r + x2r
a(17) = x1i + x2i
a(18) = x1r - x2r
a(19) = x1i - x2i
x1r = y8r - y10r
x1i = y8i - y10i
x2r = y9r + y11r
x2i = y9i + y11i
a(20) = x1r - x2i
a(21) = x1i + x2r
a(22) = x1r + x2i
a(23) = x1i - x2r
x1r = y12r - y14i
x1i = y12i + y14r
x0r = y13r + y15i
x0i = y13i - y15r
x2r = wn4r * (x0r - x0i)
x2i = wn4r * (x0i + x0r)
a(24) = x1r + x2r
a(25) = x1i + x2i
a(26) = x1r - x2r
a(27) = x1i - x2i
x1r = y12r + y14i
x1i = y12i - y14r
x0r = y13r - y15i
x0i = y13i + y15r
x2r = wn4r * (x0r - x0i)
x2i = wn4r * (x0i + x0r)
a(28) = x1r - x2i
a(29) = x1i + x2r
a(30) = x1r + x2i
a(31) = x1i - x2r
end
!
subroutine cftf081(a, w)
real*8 a(0 : 15), w(0 : *)
real*8 wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
wn4r = w(1)
x0r = a(0) + a(8)
x0i = a(1) + a(9)
x1r = a(0) - a(8)
x1i = a(1) - a(9)
x2r = a(4) + a(12)
x2i = a(5) + a(13)
x3r = a(4) - a(12)
x3i = a(5) - a(13)
y0r = x0r + x2r
y0i = x0i + x2i
y2r = x0r - x2r
y2i = x0i - x2i
y1r = x1r - x3i
y1i = x1i + x3r
y3r = x1r + x3i
y3i = x1i - x3r
x0r = a(2) + a(10)
x0i = a(3) + a(11)
x1r = a(2) - a(10)
x1i = a(3) - a(11)
x2r = a(6) + a(14)
x2i = a(7) + a(15)
x3r = a(6) - a(14)
x3i = a(7) - a(15)
y4r = x0r + x2r
y4i = x0i + x2i
y6r = x0r - x2r
y6i = x0i - x2i
x0r = x1r - x3i
x0i = x1i + x3r
x2r = x1r + x3i
x2i = x1i - x3r
y5r = wn4r * (x0r - x0i)
y5i = wn4r * (x0r + x0i)
y7r = wn4r * (x2r - x2i)
y7i = wn4r * (x2r + x2i)
a(8) = y1r + y5r
a(9) = y1i + y5i
a(10) = y1r - y5r
a(11) = y1i - y5i
a(12) = y3r - y7i
a(13) = y3i + y7r
a(14) = y3r + y7i
a(15) = y3i - y7r
a(0) = y0r + y4r
a(1) = y0i + y4i
a(2) = y0r - y4r
a(3) = y0i - y4i
a(4) = y2r - y6i
a(5) = y2i + y6r
a(6) = y2r + y6i
a(7) = y2i - y6r
end
!
subroutine cftf082(a, w)
real*8 a(0 : 15), w(0 : *)
real*8 wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i
real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
wn4r = w(1)
wk1r = w(2)
wk1i = w(3)
y0r = a(0) - a(9)
y0i = a(1) + a(8)
y1r = a(0) + a(9)
y1i = a(1) - a(8)
x0r = a(4) - a(13)
x0i = a(5) + a(12)
y2r = wn4r * (x0r - x0i)
y2i = wn4r * (x0i + x0r)
x0r = a(4) + a(13)
x0i = a(5) - a(12)
y3r = wn4r * (x0r - x0i)
y3i = wn4r * (x0i + x0r)
x0r = a(2) - a(11)
x0i = a(3) + a(10)
y4r = wk1r * x0r - wk1i * x0i
y4i = wk1r * x0i + wk1i * x0r
x0r = a(2) + a(11)
x0i = a(3) - a(10)
y5r = wk1i * x0r - wk1r * x0i
y5i = wk1i * x0i + wk1r * x0r
x0r = a(6) - a(15)
x0i = a(7) + a(14)
y6r = wk1i * x0r - wk1r * x0i
y6i = wk1i * x0i + wk1r * x0r
x0r = a(6) + a(15)
x0i = a(7) - a(14)
y7r = wk1r * x0r - wk1i * x0i
y7i = wk1r * x0i + wk1i * x0r
x0r = y0r + y2r
x0i = y0i + y2i
x1r = y4r + y6r
x1i = y4i + y6i
a(0) = x0r + x1r
a(1) = x0i + x1i
a(2) = x0r - x1r
a(3) = x0i - x1i
x0r = y0r - y2r
x0i = y0i - y2i
x1r = y4r - y6r
x1i = y4i - y6i
a(4) = x0r - x1i
a(5) = x0i + x1r
a(6) = x0r + x1i
a(7) = x0i - x1r
x0r = y1r - y3i
x0i = y1i + y3r
x1r = y5r - y7r
x1i = y5i - y7i
a(8) = x0r + x1r
a(9) = x0i + x1i
a(10) = x0r - x1r
a(11) = x0i - x1i
x0r = y1r + y3i
x0i = y1i - y3r
x1r = y5r + y7r
x1i = y5i + y7i
a(12) = x0r - x1i
a(13) = x0i + x1r
a(14) = x0r + x1i
a(15) = x0i - x1r
end
!
subroutine cftf040(a)
real*8 a(0 : 7), x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
x0r = a(0) + a(4)
x0i = a(1) + a(5)
x1r = a(0) - a(4)
x1i = a(1) - a(5)
x2r = a(2) + a(6)
x2i = a(3) + a(7)
x3r = a(2) - a(6)
x3i = a(3) - a(7)
a(0) = x0r + x2r
a(1) = x0i + x2i
a(2) = x1r - x3i
a(3) = x1i + x3r
a(4) = x0r - x2r
a(5) = x0i - x2i
a(6) = x1r + x3i
a(7) = x1i - x3r
end
!
subroutine cftb040(a)
real*8 a(0 : 7), x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
x0r = a(0) + a(4)
x0i = a(1) + a(5)
x1r = a(0) - a(4)
x1i = a(1) - a(5)
x2r = a(2) + a(6)
x2i = a(3) + a(7)
x3r = a(2) - a(6)
x3i = a(3) - a(7)
a(0) = x0r + x2r
a(1) = x0i + x2i
a(2) = x1r + x3i
a(3) = x1i - x3r
a(4) = x0r - x2r
a(5) = x0i - x2i
a(6) = x1r - x3i
a(7) = x1i + x3r
end
!
subroutine cftx020(a)
real*8 a(0 : 3), x0r, x0i
x0r = a(0) - a(2)
x0i = a(1) - a(3)
a(0) = a(0) + a(2)
a(1) = a(1) + a(3)
a(2) = x0r
a(3) = x0i
end
!
subroutine rftfsub(n, a, nc, c)
integer n, nc, j, k, kk, ks, m
real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr, xi, yr, yi
m = n / 2
ks = 2 * nc / m
kk = 0
do j = 2, m - 2, 2
k = n - j
kk = kk + ks
wkr = 0.5d0 - c(nc - kk)
wki = c(kk)
xr = a(j) - a(k)
xi = a(j + 1) + a(k + 1)
yr = wkr * xr - wki * xi
yi = wkr * xi + wki * xr
a(j) = a(j) - yr
a(j + 1) = a(j + 1) - yi
a(k) = a(k) + yr
a(k + 1) = a(k + 1) - yi
end do
end
!
subroutine rftbsub(n, a, nc, c)
integer n, nc, j, k, kk, ks, m
real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr, xi, yr, yi
m = n / 2
ks = 2 * nc / m
kk = 0
do j = 2, m - 2, 2
k = n - j
kk = kk + ks
wkr = 0.5d0 - c(nc - kk)
wki = c(kk)
xr = a(j) - a(k)
xi = a(j + 1) + a(k + 1)
yr = wkr * xr + wki * xi
yi = wkr * xi - wki * xr
a(j) = a(j) - yr
a(j + 1) = a(j + 1) - yi
a(k) = a(k) + yr
a(k + 1) = a(k + 1) - yi
end do
end
!
subroutine dctsub(n, a, nc, c)
integer n, nc, j, k, kk, ks, m
real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr
m = n / 2
ks = nc / n
kk = 0
do j = 1, m - 1
k = n - j
kk = kk + ks
wkr = c(kk) - c(nc - kk)
wki = c(kk) + c(nc - kk)
xr = wki * a(j) - wkr * a(k)
a(j) = wkr * a(j) + wki * a(k)
a(k) = xr
end do
a(m) = c(0) * a(m)
end
!
subroutine dstsub(n, a, nc, c)
integer n, nc, j, k, kk, ks, m
real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr
m = n / 2
ks = nc / n
kk = 0
do j = 1, m - 1
k = n - j
kk = kk + ks
wkr = c(kk) - c(nc - kk)
wki = c(kk) + c(nc - kk)
xr = wki * a(k) - wkr * a(j)
a(k) = wkr * a(k) + wki * a(j)
a(j) = xr
end do
a(m) = c(0) * a(m)
end
!