| 1 | ! Transformee de Fourier complexe-reel multiple |
|---|
| 2 | |
|---|
| 3 | program tjmrfftmlt1 |
|---|
| 4 | |
|---|
| 5 | implicit none |
|---|
| 6 | |
|---|
| 7 | integer, parameter :: m = 7 |
|---|
| 8 | |
|---|
| 9 | integer, parameter :: n = 36 |
|---|
| 10 | complex, dimension(0:n/2,0:m-1) :: x, xx |
|---|
| 11 | real, dimension(0:n/2,0:m-1,2) :: rx |
|---|
| 12 | equivalence ( x, rx ) |
|---|
| 13 | real, dimension(0:2*(n/2+1)-1,0:m-1) :: y |
|---|
| 14 | equivalence(x(0,0),y(0,0)) |
|---|
| 15 | |
|---|
| 16 | integer, parameter :: ntrigs = 2*n |
|---|
| 17 | real, dimension(0:ntrigs-1) :: trigs |
|---|
| 18 | |
|---|
| 19 | integer, parameter :: nifax = 19 |
|---|
| 20 | real, dimension(0:nifax-1) :: ifax |
|---|
| 21 | |
|---|
| 22 | integer, parameter :: nwork = 2*m*n |
|---|
| 23 | real, dimension(0:nwork-1) :: work |
|---|
| 24 | |
|---|
| 25 | integer :: isign |
|---|
| 26 | integer :: i, j, k |
|---|
| 27 | real :: twopi |
|---|
| 28 | complex :: s |
|---|
| 29 | integer :: inc, jump |
|---|
| 30 | |
|---|
| 31 | twopi = 2 * acos(real(-1)) |
|---|
| 32 | |
|---|
| 33 | ! On prepare le tableau d'entree en forcant a 0 les termes necessaires |
|---|
| 34 | call random_number( rx ) |
|---|
| 35 | x(0,:) = cmplx( real( x(0,:) ), 0 ) |
|---|
| 36 | if ( mod( n, 2 ) == 0 ) x(n/2,:) = cmplx( real( x(n/2,:) ), 0 ) |
|---|
| 37 | xx = x |
|---|
| 38 | |
|---|
| 39 | call fftfax(n,ifax,trigs) |
|---|
| 40 | isign = +1 |
|---|
| 41 | print *,'jmrfftmlt ',n,m,isign |
|---|
| 42 | inc = 1 |
|---|
| 43 | jump = 2*(n/2+1) |
|---|
| 44 | call rfftmlt(x,work,trigs,ifax,inc,jump,n,m,isign) |
|---|
| 45 | |
|---|
| 46 | ! On imprime le tableau de sortie |
|---|
| 47 | open(10,file='temp1',status='unknown',form='formatted') |
|---|
| 48 | write(10,'(e25.12)') y(0:n-1,0:m-1) |
|---|
| 49 | |
|---|
| 50 | ! Ce qu'il faut trouver |
|---|
| 51 | open(11,file='temp2',status='unknown',form='formatted') |
|---|
| 52 | |
|---|
| 53 | ! On reprepare le tableau d'entree |
|---|
| 54 | x = xx |
|---|
| 55 | ! Et on calcule |
|---|
| 56 | do j = 0,m-1 |
|---|
| 57 | do i = 0,n-1 |
|---|
| 58 | s = 0 |
|---|
| 59 | do k = 0,n/2 |
|---|
| 60 | s = s+cmplx(cos(twopi*i*k/real(n)),sin(twopi*i*k/real(n)))*x(k,j) |
|---|
| 61 | end do |
|---|
| 62 | do k = n/2+1,n-1 |
|---|
| 63 | s = s+cmplx(cos(twopi*i*k/real(n)),sin(twopi*i*k/real(n)))* & |
|---|
| 64 | & conjg(x(n-k,j)) |
|---|
| 65 | end do |
|---|
| 66 | write(11,'(e25.12)') real(s) |
|---|
| 67 | end do |
|---|
| 68 | end do |
|---|
| 69 | close(11) |
|---|
| 70 | |
|---|
| 71 | end program tjmrfftmlt1 |
|---|
