c 03/10/8版
c ppm.f の出力として得られる画像の横縞成分yokosima.dを
c   貴家仁志『よくわかるディジタル画像処理』(CQ出版社)
c   付属の1次元FFTプログラム（fftvec）
c に入力して得られる周波数スペクトルデータを yokospe.d とし、
c ppm.fの出力として得られる画像の縦縞成分tatesima.dを
c fftvec に入力して得られる周波数スペクトルデータをtatespe.dとする。
c
c yokospe.d,tatespe.d それぞれのスペクトルは、角周波数:-π<Ω<πの
c データをΩ=0の軸で折り返してΩの正負について平均を取り、
c かつスペクトルは２乗してパワーにする。
c Ngraphで作図するためのデータとして、0<Ω<πの角周波数と
c 横縞パワーと縦縞パワーをそれぞれ常用対数を取ってloglog.dに出力。
c loglogプロットした場合の横縞パワーの回帰直線の傾きと
c 縦縞パワーの回帰直線の傾きの平均、傾きの差、標準偏差の平均を出力。
c 例えば、gazou01.ppm という画像データに対しては、
c ppm < gazou01.ppm > rgb.01
c fftvec yokosima.d > yokospe.d
c fftvec tatesima.d > tatespe.d
c spe2dg > spe.01
c みたいな手順で実行する。
c
c
       program spe2dg
       implicit real*8 (a-h, o-z)
       dimension yoko(1024),yoko2(1024),tate(1024),tate2(1024),
     &             x(512),y(512),t(512)
c
c 求まる周波数スペクトルのデータの両端（Ω=0,2π）のプロットは、
c 他のプロット群から数桁のオーダーでかけ離れて（絶対値が）大きな値に
c なることがc あり、その1点だけで回帰直線の傾きを左右することがあるので、
c 両端のデータを捨てる場合は、端から何個までのデータを捨てるかを
c nskipで指定する。捨てない場合は nskip=0
       nskip=1
c
       open(7,file='yokospe.d')
         read(7,*) n
        do i=1,n
         read(7,*) yoko(i),yoko2(i)
        end do
       close(7)
c
       open(7,file='tatespe.d')
         read(7,*) n
        do i=1,n
         read(7,*) tate(i),tate2(i)
        end do
       close(7)
c
         pi=2.d0*asin(1.d0)
         n2=n/2
       open(7,file='loglog.d')
          n1=nskip+1
         do i=n1,n2
          y2=( yoko(i)**2 + yoko(n-i+1)**2 )/2.d0
          t2=( tate(i)**2 + tate(n-i+1)**2 )/2.d0
          omega=float(i)/float(n)*pi*2.d0
           in=i-nskip
          y(in)=log10(y2)
          t(in)=log10(t2)
          x(in)=log10(omega)
          write(7,200) x(in),y(in),t(in)
         end do
       close(7)
c
          mn=n2-nskip
          call heikin(mn,x,xm,xs)
          call heikin(mn,y,ym,ys)
          call heikin(mn,t,tm,ts)
c
           xy=0.d0
           xt=0.d0
          do i=1,n2
           xy=xy+x(i)*y(i)
           xt=xt+x(i)*t(i)
          end do
           sxy=xy/n2-xm*ym
           sxt=xt/n2-xm*tm
           ayoko=sxy/xs**2
           atate=sxt/xs**2
           a=(ayoko+atate)/2.d0
           hensa=(ys+ts)/2.d0
           sa=ayoko-atate
            write(*,*) ' 両端データスキップ数：',nskip
            write(*,*) ' 横縞スペクトル回帰直線傾き：',ayoko
            write(*,*) ' 横縞スペクトル回帰直線偏差：',ys
            write(*,*) ' 縦縞スペクトル回帰直線傾き：',atate
            write(*,*) ' 縦縞スペクトル回帰直線偏差：',ts
            write(*,*) '     縦横平均傾き：',a
            write(*,*) ' 縦横平均標準偏差：',hensa
            write(*,*) '       縦横傾き差：',sa
200      format(1p3d13.5)
         end
c
          subroutine heikin(mn,x,a,s)
       implicit real*8 (a-h, o-z)
       dimension x(mn)
         a=0.d0
         s=0.d0
        do i=1,mn
         a=a+x(i)
        end do
         a=a/float(mn)
        do i=1,mn
         s=s+ ( x(i)-a )**2
        end do
         s=sqrt(s/mn)
        return
        end
c