求助，卷积的小例子

wbaby91

各位，中午好！
最近在学习fortran，学习版本IVF11，VS2008, 附件为整个项目文件
我在书上找到了卷积的function代码，还有几个必须的模块，但是在写小例子主程序的时候遇到了障碍，
请大神不吝赐教，帮我修改一下，感激不尽，或者写一个简短的能运行的小例子启发我一下也可以
即就是data1=（0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0），respns=（0 0 1 1 1 1 0 0 0）求两者卷积

[Fortran] 纯文本查看 复制代码

PROGRAM CONV_MAIN  
parameter (N=16,M=9)
dimension data1(N),respns(M),tmpr(m+n-1)
DO i=1,N
data1(I)=0.0
IF ((I>=(N/2-N/8)).and.(I<=(N/2+N/8)))&
DATA1(I)=1.0
END DO
DO I=1,M
RESPNS(I)=0.0
IF(I>2.AND.I<7)RESPNS(I)=1.0
END DO
ISIGN=1
tmpr=convlv(DATA1,RESPNS,isign)
WRITE(*,'(/1X,T4,A,T13,A)')'I','TMPR'
END PROGRAM CONV_MAIN

楚香饭

在你的代码里。realft 这个函数用来求取傅氏变换。但是代码却是空白的。
nr.f90 里面只是接口，应该提供的是静态库。这里的代码只有接口，没有内部实现，不能直接复制过来用。

关于卷积呢。有两种：线性卷积和圆周卷积。根据你的结果是  n+m-1 来看，应该是想求取线性卷积。
而如果采用 FFT 的方式来求取，则是圆周卷积。
（你可以用补一个周期的 0 做圆周卷积，来代替线性卷积）

所以，一共有4种方法：
1.直接求取圆周卷积。（考虑周期的两层循环）
2.直接求取线性卷积。（两层循环）
3.傅氏变换求取圆周卷积。（傅氏变换的乘积，等于卷积的傅氏变换）
4.傅氏变换求取线性卷积。（补一个周期的0，再求第3方法）

在这里，给出2，4两种的代码。（我自己写的，没有用到你给的 nr 代码）

方法2：

[Fortran] 纯文本查看 复制代码

Program conv_main_fft
  Implicit None
  
  Interface 
    FUNCTION convlv(x1,x2)
      REAL, DIMENSION(:), INTENT(IN) :: x1 , x2
      REAL, DIMENSION(size(x1)+size(x2)-1) :: convlv
    End Function convlv
  End Interface
  
  Integer , Parameter :: n=16, m=9
  Real :: data1(2*n), respns(m), tmpr(n+m-1)
  integer :: i
  data1(:) = 0.0
  data1( (n/2-n/8) : (n/2+n/8) ) = 1.0
  respns(:) = 0.0
  respns( 3:6 ) = 1.0
  tmpr = 0.0
  tmpr = convlv(data1, respns)
  Write (*, '(9999(f5.1,1x))' ) tmpr
End Program conv_main_fft

  Function convlv( X1 , X2 )
    Real , Intent( IN ) :: x1(:) , x2(:)
    REAL, DIMENSION(size(x1)+size(x2)-1) :: convlv
    integer :: i , j , m , n , t
	  real:: rTemp
    m = size(x1) ; n = size(x2)
    do i = 1 , m+n-1
      convlv(i) = 0.0
	    do j = 1 , m
        if( (i-j)>=1 .and. (i-j) <= n )  then
          convlv(i)=convlv(i) + x1(j) * x2(i-j)
        end if
	    end do
    end do
  End Function convlv

方法4：

[Fortran] 纯文本查看 复制代码

Module FFT_Mod
  Implicit None
  Integer , parameter :: DP = Selected_Real_Kind( 5 )
  Integer , parameter :: FFT_Forward = 1
  Integer , parameter :: FFT_Inverse = -1
contains

  Subroutine fcFFT( x , forback )
    !//Subroutine FFT , Cooley-Tukey , radix-2
    !// [url]www.fcode.cn[/url]
    Real(Kind=DP) , parameter :: PI = 3.141592654_DP
    Complex(Kind=DP) , Intent(INOUT) :: x(:)
    Integer , Intent(IN) :: forback
    Integer :: n
    integer :: i , j , k , ncur , ntmp , itmp
    real(Kind=DP) :: e
    complex(Kind=DP) :: ctmp
    n = size(x)
    ncur = n
    Do
      ntmp = ncur
      e = 2.0_DP * PI / ncur
      ncur = ncur / 2
      if ( ncur < 1 ) exit
      Do j = 1 , ncur
        Do i = j , n , ntmp
          itmp = i + ncur
          ctmp = x(i) - x(itmp)
          x(i) = x(i) + x(itmp)
          x(itmp) = ctmp * exp( forback * cmplx( 0.0_DP , e*(j-1) ) )
        End Do   
      End Do
    End Do
    j = 1
    Do i = 1, n - 1
      If ( i < j ) then
        ctmp = x(j)
        x(j) = x(i)
        x(i) = ctmp
      End If
      k = n/2
      Do while( k < j )
        j = j - k
        k = k / 2
      End Do
      j = j + k
    End Do
    Return
  End Subroutine fcFFT

End Module FFT_Mod
  
Module typedef
  Integer , parameter :: SP   = Selected_Real_Kind( p = 5 )
  Integer , parameter :: I4b  = Selected_Int_Kind( 5 )
End Module typedef

Function convlv(x, respns, isign)
  use typedef
  use FFT_Mod
  Implicit None
  Real (Kind=sp), Dimension (:), Intent (Inout) :: x
  Real (Kind=sp), Dimension (:), Intent (In) :: respns
  Integer (Kind=i4b), Intent (In) :: isign
  Real (Kind=sp), Dimension (size(x)) :: convlv
  Integer (Kind=i4b) :: n, m
  Complex (Kind=SP), Dimension (size(x)) :: tmpd, tmpr
  tmpd = 0.0_sp
  tmpr = 0.0_sp
  n = size(x)
  m = size(respns)
  tmpd = cmplx( x , 0.0_sp )
  call fcFFT( tmpd , FFT_Forward )
  tmpr = cmplx( respns , 0.0_sp )
  call fcFFT( tmpr , FFT_Forward )
  If (isign==1) Then
    tmpr(1) = cmplx(real(tmpd(1))*real(tmpr(1))/n, aimag(tmpd(1))*aimag(tmpr(1))/n, kind=sp)
    tmpr(2:) = tmpd(2:)*tmpr(2:)/n
  Else If (isign==-1) Then
    If (any(abs(tmpr(2:))==0.0) .Or. real(tmpr(1))==0.0 .Or. aimag(tmpr(1))==0.0) write(*,*) 'deconvolving at response zero in convlv'
    tmpr(1) = cmplx(real(tmpd(1))/real(tmpr(1))/n, aimag(tmpd(1))/aimag(tmpr(1))/n, kind=sp)
    tmpr(2:) = tmpd(2:)/tmpr(2:)/n
  Else
    write(*,*) 'no meaning for isign in convlv'
  End If
  Call fcFFT( tmpr , FFT_Inverse )
  convlv = real( tmpr )  
End Function convlv

Program conv_main_fft
  use typedef
  Implicit None
  
  Interface 
    FUNCTION convlv(data,respns,isign)
      use typedef
      REAL(SP), DIMENSION(:), INTENT(INOUT) :: data
      REAL(SP), DIMENSION(:), INTENT(IN) :: respns
      INTEGER(I4B), INTENT(IN) :: isign
      REAL(SP), DIMENSION(size(data)) :: convlv
    End Function convlv
  End Interface
  
  Integer , Parameter :: n=16, m=16
  Real(Kind=SP) :: data1(2*n), respns(m), tmpr(2*n)
  integer :: i , isign
  data1(:) = 0.0_SP
  data1( (n/2-n/8) : (n/2+n/8) ) = 1.0_SP
  respns(:) = 0.0_SP
  respns( 3:6 ) = 1.0_SP
  isign = 1
  tmpr = 0.0_SP
  tmpr = convlv(data1, respns, isign)
  Write (*, '(9999(f5.1,1x))' ) tmpr
End Program conv_main_fft

wbaby91

恩恩，基本上是想做一维或者二维卷积，
程序还有一点小问题：错误 error #8000:  There is a conflict between local interface block and external interface block.   [CONVLV]       
两部分调试都出现了这种情况，不知道是不是我编译器的问题

楚香饭

我用多个编译器都编译过，没有问题。你再自己看看吧。

也许清理一下工程对你有帮助。

wbaby91

恩，使用CVF就没问题，谢谢群主，IVF我再看看我是不是哪设置出问题了

wbaby91

我把楚帅的代码中的外部函数都扔到模组里，然后删掉了接口块，在我的Intel版本里也能运行了，应该没啥问题，希望以后继续向论坛里的朋友学习，等姿势水平提高后，帮助更多热爱科学计算的朋友