home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / fft2 / noise.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  5.6 KB  |  258 lines
  1. /* *****************************************************************************
  2. *
  3. * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  4. * All Rights Reserved.
  5. *
  6. * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  7. * the contents of this file may not be disclosed to third parties, copied or
  8. * duplicated in any form, in whole or in part, without the prior written
  9. * permission of Silicon Graphics, Inc.
  10. *
  11. * RESTRICTED RIGHTS LEGEND:
  12. * Use, duplication or disclosure by the Government is subject to restrictions
  13. * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  14. * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  15. * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  16. * rights reserved under the Copyright Laws of the United States.
  17. *
  18. ***************************************************************************** */
  19.       subroutine add_circle( nlin, ncol, mat, lda)
  20.       real mat(lda,ncol)
  21.  
  22.     ior = nlin / 2
  23.     jor = ncol / 2
  24.  
  25.     r = MIN(ior,jor)/9.
  26.     r2 = r * r
  27.  
  28. c$doacross local(i,j,dx,dy,d2,z)
  29.     do j = 1, ncol
  30.         do i = 1, nlin
  31.         dy = j - jor
  32.         dx = i - ior
  33.         d2 = dx*dx + dy*dy
  34.         z = (.01 + d2) / r2
  35.         mat(i,j) = Sin(z) / z
  36. c        mat(i,j) = Sin(z) * Cos(dx/r) / z
  37.         end do
  38.     enddo
  39.  
  40.       return
  41.       end
  42.  
  43.  
  44.       subroutine inifil(a,lda,nlin,ncol)
  45.       real a(lda,1)
  46. c************************************************
  47. c                        *
  48. c   create a wave propagation filter        *
  49. c                        *
  50. c   The Fourier Transform of the 2D wave    *
  51. c   equation link the "time" frequency W, with    *
  52. c   the "space" frequencies KX and KY        *
  53. c                        *
  54. c    W^2 == c^2 * (KX^2 + KY^2), c velocity    *
  55. c                        *
  56. c   The wave propagation between time 0 and T     *
  57. c   simulated is given by the phase shift ope-    *
  58. c   rator   Exp( - i * W * T)            *
  59. c   which translate in :            *
  60. c    Exp( - i * c * Sqrt(KX^2+KY^2) * T    *
  61. c   if we consider the forward propagation,    *
  62. c   and its conjugate for the backward        *
  63. c                        *
  64. c************************************************
  65. c                        *
  66. c   Here we choose :                *
  67. c    c = 1                    *
  68. c    T = 1/ ncol                *
  69. c                        *
  70. c************************************************
  71.     i_lin_even= mod(nlin+1,2)
  72.     i_col_even= mod(ncol+1,2)
  73.     T = 1 / Float(MIN(ncol,nlin))
  74. c
  75. c   first line  ( KX == 0 )
  76. c
  77.     KX = 0.
  78.     a(1,1) = 1.0
  79.     do KY = 1, (ncol-1)/2
  80.         y = Float(KY)
  81.         W = Sqrt( y*y )
  82.         a(1,2*KY) = cos(W*T)
  83.         a(1,2*KY+1) = - sin(W*T)
  84.     end do
  85.     if( i_col_even) then
  86.         y = Float(ncol/2)
  87.         W = Sqrt( y*y )
  88.         a(1,ncol) = cos(W*T)
  89.     endif
  90. c
  91. c   General Values
  92. c
  93.     do KY = 0, (ncol-1)
  94.         y = Float(KY)
  95.         do KX = 1, (nlin-1)/2
  96.         x = Float(KX)
  97.         W = Sqrt(x*x + y*y)
  98.         a(2*KX,KY+1) = cos(W*T)
  99.         a(2*KX+1,KY+1) = - sin(W*T)
  100.         end do
  101.     end do
  102. c
  103. c   first line  ( KX == nlin/2 )
  104. c
  105.       if( i_lin_even ) then
  106.     KX = (nlin/2)
  107.     x = Float(KX)
  108.     a(ncol,1) = 1.0
  109.     do KY = 1, (ncol-1)/2
  110.         y = Float(KY)
  111.         W = Sqrt( y*y + x*x)
  112.         a(ncol,2*KY) = cos(W*T)
  113.         a(ncol,2*KY+1) = - sin(W*T)
  114.     end do
  115.     if( i_col_even) then
  116.         y = Float(ncol/2)
  117.         W = Sqrt( y*y + x*x)
  118.         a(ncol,ncol) = cos(W*T)
  119.     endif
  120.       endif
  121.  
  122.       return
  123.       end
  124.       
  125.       subroutine reverse(a,lda,nlin,ncol)
  126.       real a(lda,1)
  127.     i_lin_even= mod(nlin+1,2)
  128. c
  129. c   first line  
  130. c
  131.     do j = 3, (ncol-1), 2
  132.         a(1,j) = - a(1,j)
  133.     end do
  134. c
  135. c   General Values
  136. c
  137.     do j = 1, ncol
  138.         do i = 3, (nlin -1), 2
  139.         a(i,j) = - a(i,j)
  140.         end do
  141.     enddo
  142. c
  143. c   last line  
  144. c
  145.     if( i_lin_even) then
  146.         do j = 3, (ncol-1), 2
  147.         a(nlin,j) = - a(nlin,j)
  148.         end do
  149.     end if
  150.       return
  151.       end
  152.       
  153.       subroutine matnoise(a,lda,nlin,ncol)
  154.       real a(lda,1)
  155. c
  156.       init = 1325
  157.  
  158.       do 30 j = 1,ncol
  159.          do 20 i = 1,nlin
  160.             init = mod(3125*init,65536)
  161.             a(i,j) = a(i,j) + (init - 32768.0)/16384.0
  162.    20    continue
  163.    30    continue
  164.       return
  165.       end
  166.       
  167.       subroutine matpower( nlin,ncol,mat,power,lda)
  168.       real mat(lda,*), power(lda,*)
  169.     real re, im
  170.  
  171.         il = 1
  172.         power(il,1) = mat(il,1) * mat(il,1)
  173.         do j = 2, (ncol-1)/2
  174.         re = mat(il,2*j-1)
  175.         im = mat(il,2*j)
  176.         power(il,j) = re * re + im * im
  177.         power(il,ncol-j+1) = power(il,j)
  178.         end do
  179.         if( mod(ncol,2) .eq. 2) 
  180.      $        power(il,ncol) = mat(il,ncol) * mat(il,ncol)
  181. c$doacross local(i,j,re,im)
  182.     do j = 1, ncol
  183.         do i = 2, (nlin-1)/2
  184.         re = mat(2*i-2,j)
  185.         im = mat(2*i-1,j)
  186.         power(i,j) = re * re + im * im
  187.         power(nlin-i+1,j) = power(i,j)
  188.         end do
  189.     end do
  190.         il = nlin
  191.         power(il,1) = mat(il,1) * mat(il,1)
  192.         do j = 2, (ncol-1)/2
  193.         re = mat(il,2*j-1)
  194.         im = mat(il,2*j)
  195.         power(il,j) = re * re + im * im
  196.         power(il,ncol-j+1) = power(il,j)
  197.         end do
  198.         if( mod(ncol,2) .eq. 2) 
  199.      $        power(il,ncol) = mat(il,ncol) * mat(il,ncol)
  200.     if ( mod(nlin ,2) .eq. 0) then
  201.     end if
  202.       return
  203.       end
  204.  
  205.       subroutine matlog( nlin, ncol, power, lda)
  206.       real power(lda,*)
  207.     real t
  208. c$doacross local(i,j,t)
  209.     do j = 1, ncol
  210.         do i = 1, nlin
  211.         t = power(i,j) + 1
  212.         power(i,j) = LOG(t)
  213.         end do
  214.     end do
  215.       return
  216.       end
  217.         
  218.       subroutine matcolor( nlin,ncol,mat,color,lda)
  219.       real mat(lda,*)
  220.       integer color(nlin,*)
  221.     real t, a, coeff,xx,yy
  222.     integer red, green, blue
  223.  
  224.     xx = 0.0
  225.     yy = 0.0
  226.     do j = 1, ncol
  227.         do i = 1, nlin
  228.         a = mat(i,j)
  229.         if ( a .lt. xx) then
  230.             xx = a
  231.         else if(  a .gt. yy) then
  232.             yy   = a
  233.         end if
  234.         end do
  235.     end do
  236.     coeff = MAX(ABS(xx),ABS(yy))
  237.  
  238. c    print *, ' MIN =', xx, '   MAX =', yy
  239.  
  240.     if ( coeff .gt. 0.) coeff = .99 * 255. / coeff
  241.  
  242. c$doacross local(t, a, red,blue,green,i,j)
  243.     do j = 1, ncol
  244.         do i = 1, nlin
  245.         t= mat(i,j)
  246.         a = ABS(t)
  247.         red = coeff* (a - t)
  248.         if( red .gt. 255) red = 255
  249.         green = coeff* a
  250.         blue = coeff* (a + t)
  251.         if( blue .gt. 255) blue = 255
  252.         color(i,j) = (red) + (256*green) + (256*256*blue)
  253.         end do
  254.     end do
  255.       return
  256.       end
  257.         
  258.