fft.for

此程序用Fortran语言编写

	program FFT
	  integer p,n,nu
	  integer,allocatable::rev_array(:)
	  complex,allocatable::x(:)
	  character*(40)cmdfile,infile,outfile
	  write(*,*)'input cmdfilename'
	  read(*,*)cmdfile
	  call cmd(cmdfile,n,p,nu,infile,outfile)
	  allocate (rev_array(0:nu),x(0:nu))
	  call read_data(nu,x,infile)
	  call reverse(rev_array,x,nu,p)
	  call sub_fft(nu,p,x)
	  call output(nu,x,outfile)
	  deallocate(rev_array,x)
	end program
	  
	subroutine cmd(cmdfile,n,p,nu,infile,outfile)
	  character*(*)cmdfile,infile,outfile
	  integer n,p,nu
	  open(15,file=cmdfile)
	  read(15,*)n
	  read(15,*)infile
        read(15,*)outfile
	  p=log(real(n))/log(2.0)+0.01
	  nu=n-1
	  close(15)
	end subroutine

	subroutine read_data(nu,x,infile)
	  integer nu
	  character *(*) infile
	  real a(0:nu)
	  complex x(0:nu)
	  open (11,file=infile)
	  read(11,*)(x(i),i=0,nu)
        close(11)
	end subroutine

	subroutine reverse(rev,x,nu,p)
	  integer nu,p
	  integer rev(0:nu)
	  complex x(0:nu)
	  do i=0,nu
	    nt=i;rev(i)=0
	    do j=1,p
	      rev(i)=rev(i)+mod(nt,2)*2**(p-j)
	      nt=nt/2
	    end do
          if(rev(i)>i)then
	      t=x(i)
		  x(i)=x(rev(i)) 
            x(rev(i))=t
	    end if 
	  end do
	end subroutine

	subroutine sub_ fft(nu,p,x)
	  parameter pi=3.1415926
	  integer nu,p
	  complex x(0:nu)
	  complex s,w
	  do i=1,p
	    do j=1,(nu+1)/2**i
	      do k=0,2**(i-1)-1
	        w=cmplx(cos(2*pi*k/2**i),-sin(2*pi*k/2**i))
	        s=x((j-1)*2**i+k+2**(i-1))*w
	        x((j-1)*2**i+k+2**(i-1))=x((j-1)*2**i+k)-s
              x((j-1)*2**i+k)=x((j-1)*2**i+k)+s
	      end do
	    end do
	  end do
	  write(*,*) x
	end subroutine

	subroutine output(nu,x,outfile)
	  integer nu
	  complex x(0:nu)
	  character *(*) outfile
	  open(20,file=outfile)
	  close(20,status='delete')
	  open(20,file=outfile)
	  do i=0,nu
	  write(20,10)i,x(i)
 10	  format(1x,'x(',i1,')=','('f9.6,',',f9.6')')
        end do
	  close(20)
	end subroutine