[Search for users]
[Overall Top Noters]
[List of all Conferences]
[Download this site]
| Title: | Digital Fortran |
| Notice: | Read notes 1.* for important information |
| Moderator: | QUARK::LIONEL |
|
| Created: | Thu Jun 01 1995 |
| Last Modified: | Fri Jun 06 1997 |
| Last Successful Update: | Fri Jun 06 1997 |
| Number of topics: | 1333 |
| Total number of notes: | 6734 |
1246.0. "fft again ..." by RTOMS::PARETIJ () Sun Mar 30 1997 16:46
The following code does some sort of fft; there
are 2 routines fft_ce and FFT_se calling cosft1 (sinft1) and
realft calling four1. There are no comments, which strains
my imagination. I am again looking for a way to use dxml.
May be putting the code out to the forum will generate some
questions I could refer to the customer.
subroutine fft_ce(A,Re_,Im_,nx,ny,isign)
implicit none
integer nx,ny,isign
real*8 A(0:nx,0:ny-1),Re_(0:nx,0:ny/2),Im_(0:nx,0:ny/2)
integer i,j
real*8 vec_x(0:nx), vec_y(0:ny-1)
if(isign.eq.1)then
do j=0,ny-1
do i=0,nx
vec_x(i)=A(i,j)
enddo
call cosft1(vec_x,nx)
A(0,j)=vec_x(0)/nx
do i=1,nx-1
A(i,j)=vec_x(i)*2/nx
enddo
A(nx,j)=vec_x(nx)/nx
enddo
do i=0,nx
do j=0,ny-1
vec_y(j)=A(i,j)
enddo
call realft(vec_y,ny,isign)
Re_(i,0)=vec_y(0)/ny
Im_(i,0)=0
do j=1,ny/2-1
Re_(i,j)=vec_y(2*j)/ny
Im_(i,j)=-vec_y(2*j+1)/ny
enddo
Re_(i,ny/2)=0.5*vec_y(1)/ny
Im_(i,ny/2)=0
enddo
else
do i=0,nx
vec_y(0)=2*Re_(i,0)
vec_y(1)=4*Re_(i,ny/2)
do j=1,ny/2-1
vec_y(2*j)=2*Re_(i,j)
vec_y(2*j+1)=-2*Im_(i,j)
enddo
call realft(vec_y,ny,isign)
do j=0,ny-1
A(i,j)=vec_y(j)
enddo
enddo
do j=0,ny-1
vec_x(0)=2*A(0,j)
do i=1,nx-1
vec_x(i)=A(i,j)
enddo
vec_x(nx)=2*A(nx,j)
call cosft1(vec_x,nx)
do i=0,nx
A(i,j)=vec_x(i)
enddo
enddo
endif
end
subroutine fft_se(A,Re_,Im_,nx,ny,isign)
implicit none
integer nx,ny,isign
real*8 A(0:nx,0:ny-1),Re_(0:nx,0:ny/2),Im_(0:nx,0:ny/2)
integer i,j
real*8 vec_x(0:nx-1), vec_y(0:ny-1)
if(isign.eq.1)then
do j=0,ny-1
vec_x(0)=0
do i=1,nx-1
vec_x(i)=A(i,j)
enddo
call sinft(vec_x,nx)
do i=0,nx-1
A(i,j)=vec_x(i)*2/nx
enddo
enddo
do i=1,nx-1
do j=0,ny-1
vec_y(j)=A(i,j)
enddo
call realft(vec_y,ny,isign)
Re_(i,0)=vec_y(0)/ny
Im_(i,0)=0
do j=1,ny/2-1
Re_(i,j)=vec_y(2*j)/ny
Im_(i,j)=-vec_y(2*j+1)/ny
enddo
Re_(i,ny/2)=0.5*vec_y(1)/ny
Im_(i,ny/2)=0
enddo
do j=0,ny/2
Re_(0,j)=0
Im_(0,j)=0
Re_(nx,j)=0
Im_(nx,j)=0
enddo
else
do i=1,nx-1
vec_y(0)=2*Re_(i,0)
vec_y(1)=4*Re_(i,ny/2)
do j=1,ny/2-1
vec_y(2*j)=2*Re_(i,j)
vec_y(2*j+1)=-2*Im_(i,j)
enddo
call realft(vec_y,ny,isign)
do j=0,ny-1
A(i,j)=vec_y(j)
enddo
enddo
do j=0,ny-1
vec_x(0)=0
do i=1,nx-1
vec_x(i)=A(i,j)
enddo
call sinft(vec_x,nx)
do i=0,nx-1
A(i,j)=vec_x(i)
enddo
enddo
do j=0,ny-1
A(nx,j)=0
enddo
endif
end
SUBROUTINE four1(data,nn,isign)
INTEGER isign,nn
real*8 data(2*nn)
INTEGER i,istep,j,m,mmax,n
real*8 tempi,tempr
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
n=2*nn
j=1
do 11 i=1,n,2
if(j.gt.i)then
tempr=data(j)
tempi=data(j+1)
data(j)=data(i)
data(j+1)=data(i+1)
data(i)=tempr
data(i+1)=tempi
endif
m=n/2
1 if ((m.ge.2).and.(j.gt.m)) then
j=j-m
m=m/2
goto 1
endif
j=j+m
11 continue
mmax=2
2 if (n.gt.mmax) then
istep=2*mmax
theta=6.28318530717959d0/(isign*mmax)
wpr=-2.d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
wr=1.d0
wi=0.d0
do 13 m=1,mmax,2
do 12 i=m,n,istep
j=i+mmax
tempr=(wr)*data(j)-(wi)*data(j+1)
tempi=(wr)*data(j+1)+(wi)*data(j)
data(j)=data(i)-tempr
data(j+1)=data(i+1)-tempi
data(i)=data(i)+tempr
data(i+1)=data(i+1)+tempi
12 continue
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
13 continue
mmax=istep
goto 2
endif
return
END
SUBROUTINE realft(data,n,isign)
INTEGER isign,n
real*8 data(n)
CU USES four1
INTEGER i,i1,i2,i3,i4,n2p3
real*8 c1,c2,h1i,h1r,h2i,h2r,wis,wrs
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
theta=3.141592653589793d0/dble(n/2)
c1=0.5
if (isign.eq.1) then
c2=-0.5
call four1(data,n/2,+1)
else
c2=0.5
theta=-theta
endif
wpr=-2.0d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
wr=1.0d0+wpr
wi=wpi
n2p3=n+3
do 11 i=2,n/4
i1=2*i-1
i2=i1+1
i3=n2p3-i2
i4=i3+1
wrs=(wr)
wis=(wi)
h1r=c1*(data(i1)+data(i3))
h1i=c1*(data(i2)-data(i4))
h2r=-c2*(data(i2)+data(i4))
h2i=c2*(data(i1)-data(i3))
data(i1)=h1r+wrs*h2r-wis*h2i
data(i2)=h1i+wrs*h2i+wis*h2r
data(i3)=h1r-wrs*h2r+wis*h2i
data(i4)=-h1i+wrs*h2i+wis*h2r
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
11 continue
if (isign.eq.1) then
h1r=data(1)
data(1)=h1r+data(2)
data(2)=h1r-data(2)
else
h1r=data(1)
data(1)=c1*(h1r+data(2))
data(2)=c1*(h1r-data(2))
call four1(data,n/2,-1)
endif
return
END
SUBROUTINE sinft(y,n)
INTEGER n
real*8 y(n)
CU USES realft
INTEGER j
real*8 sum,y1,y2
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
theta=3.141592653589793d0/dble(n)
wr=1.0d0
wi=0.0d0
wpr=-2.0d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
y(1)=0.0
do 11 j=1,n/2
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
y1=wi*(y(j+1)+y(n-j+1))
y2=0.5*(y(j+1)-y(n-j+1))
y(j+1)=y1+y2
y(n-j+1)=y1-y2
11 continue
call realft(y,n,+1)
sum=0.0
y(1)=0.5*y(1)
y(2)=0.0
do 12 j=1,n-1,2
sum=sum+y(j)
y(j)=y(j+1)
y(j+1)=sum
12 continue
return
END
SUBROUTINE cosft1(y,n)
INTEGER n
real*8 y(n+1)
CU USES realft
INTEGER j
real*8 sum,y1,y2
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
theta=3.141592653589793d0/n
wr=1.0d0
wi=0.0d0
wpr=-2.0d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
sum=0.5*(y(1)-y(n+1))
y(1)=0.5*(y(1)+y(n+1))
do 11 j=1,n/2-1
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
y1=0.5*(y(j+1)+y(n-j+1))
y2=(y(j+1)-y(n-j+1))
y(j+1)=y1-wi*y2
y(n-j+1)=y1+wi*y2
sum=sum+wr*y2
11 continue
call realft(y,n,+1)
y(n+1)=y(2)
y(2)=sum
do 12 j=4,n,2
sum=sum+y(j)
y(j)=sum
12 continue
return
END
| T.R | Title | User | Personal
Name | Date | Lines |
|---|
| 1246.1 | | RTL::HANEK | | Mon Mar 31 1997 11:46 | 29 |
| First off, it would be best to move this discussion to the (NICCTR::)DXML
notes file.
It appears to me that the routines cosft1, sinft1, realft and four1 are
performing text book implementations of cosine, sine, real Fourier and
complex Fourier transforms. four1( data, nn, isign ) is a standard
Cooley-Tukey in-place complex Fourier transform of length nn, with isign
indicating the direction of the transform (forward or backward). The same
calling convention extends to the other routines. Each of the routines cosft1,
sinf1, realft and four1 have direct counter parts in the DXML, so you should be
able to replace them with calls into the DXML. There may be some issues with
scaling the results by a constant depending on how the DXML routines are
implemented
There seems to be lots of opportunities for performance improvements here beyond
simply using the DXML. Two in particular seem worth mentioning:
All fft codes use "tables" of sin and cos values. The user code presentd in .-1
generates these tables on the fly which can be very inefficient if the length of
the fft remains constant over several invocations of the fft routines. The DXML
provides routines that can either generate the table on the fly, or pre-compute
it. Perhaps you should look into pre-computing the table.
Second, the fft_ce routine appears to be performing a kind of 2-dimensional
cosine tranform by transforming each of the columns and then each of the
rows. It seems to me that this operation can be implemented using a
2-dimensional real fft (which the DXML supports) by either pre-conditioning the
input data or "fixing up" the result of the 2-dimensional.
|