求助，卷积的小例子

wbaby91 · 发表于 2015-8-16 11:05:55

各位，中午好！
最近在学习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] 纯文本查看 复制代码

01PROGRAM CONV_MAIN  
02parameter (N=16,M=9)
03dimension data1(N),respns(M),tmpr(m+n-1)
04DO i=1,N
05data1(I)=0.0
06IF ((I>=(N/2-N/8)).and.(I<=(N/2+N/8)))&
07DATA1(I)=1.0
08END DO
09DO I=1,M
10RESPNS(I)=0.0
11IF(I>2.AND.I<7)RESPNS(I)=1.0
12END DO
13ISIGN=1
14tmpr=convlv(DATA1,RESPNS,isign)
15WRITE(*,'(/1X,T4,A,T13,A)')'I','TMPR'
16END PROGRAM CONV_MAIN

wbaby91 · 发表于 2015-8-17 15:15:48

我把楚帅的代码中的外部函数都扔到模组里，然后删掉了接口块，在我的Intel版本里也能运行了，应该没啥问题

，希望以后继续向论坛里的朋友学习，等姿势水平提高后，帮助更多热爱科学计算的朋友

wbaby91 · 发表于 2015-8-17 11:38:09

恩，使用CVF就没问题，谢谢群主

，IVF我再看看我是不是哪设置出问题了

楚香饭 · 发表于 2015-8-17 08:47:05

本帖最后由楚香饭于 2015-8-17 08:49 编辑

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

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

wbaby91 · 发表于 2015-8-16 15:05:02

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

楚香饭 · 发表于 2015-8-16 14:03:16

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

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

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

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

方法2：

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

01Program conv_main_fft
02  Implicit None
03   
04  Interface
05    FUNCTION convlv(x1,x2)
06      REAL, DIMENSION(:), INTENT(IN) :: x1 , x2
07      REAL, DIMENSION(size(x1)+size(x2)-1) :: convlv
08    End Function convlv
09  End Interface
10   
11  Integer , Parameter :: n=16, m=9
12  Real :: data1(2*n), respns(m), tmpr(n+m-1)
13  integer :: i
14  data1(:) = 0.0
15  data1( (n/2-n/8) : (n/2+n/8) ) = 1.0
16  respns(:) = 0.0
17  respns( 3:6 ) = 1.0
18  tmpr = 0.0
19  tmpr = convlv(data1, respns)
20  Write (*, '(9999(f5.1,1x))' ) tmpr
21End Program conv_main_fft
22 
23  Function convlv( X1 , X2 )
24    Real , Intent( IN ) :: x1(:) , x2(:)
25    REAL, DIMENSION(size(x1)+size(x2)-1) :: convlv
26    integer :: i , j , m , n , t
27    real:: rTemp
28    m = size(x1) ; n = size(x2)
29    do i = 1 , m+n-1
30      convlv(i) = 0.0
31      do j = 1 , m
32        if( (i-j)>=1 .and. (i-j) <= n )  then
33          convlv(i)=convlv(i) + x1(j) * x2(i-j)
34        end if
35      end do
36    end do
37  End Function convlv

方法4：

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

001Module FFT_Mod
002  Implicit None
003  Integer , parameter :: DP = Selected_Real_Kind( 5 )
004  Integer , parameter :: FFT_Forward = 1
005  Integer , parameter :: FFT_Inverse = -1
006contains
007 
008  Subroutine fcFFT( x , forback )
009    !//Subroutine FFT , Cooley-Tukey , radix-2
010    !// [url]www.fcode.cn[/url]
011    Real(Kind=DP) , parameter :: PI = 3.141592654_DP
012    Complex(Kind=DP) , Intent(INOUT) :: x(:)
013    Integer , Intent(IN) :: forback
014    Integer :: n
015    integer :: i , j , k , ncur , ntmp , itmp
016    real(Kind=DP) :: e
017    complex(Kind=DP) :: ctmp
018    n = size(x)
019    ncur = n
020    Do
021      ntmp = ncur
022      e = 2.0_DP * PI / ncur
023      ncur = ncur / 2
024      if ( ncur < 1 ) exit
025      Do j = 1 , ncur
026        Do i = j , n , ntmp
027          itmp = i + ncur
028          ctmp = x(i) - x(itmp)
029          x(i) = x(i) + x(itmp)
030          x(itmp) = ctmp * exp( forback * cmplx( 0.0_DP , e*(j-1) ) )
031        End Do  
032      End Do
033    End Do
034    j = 1
035    Do i = 1, n - 1
036      If ( i < j ) then
037        ctmp = x(j)
038        x(j) = x(i)
039        x(i) = ctmp
040      End If
041      k = n/2
042      Do while( k < j )
043        j = j - k
044        k = k / 2
045      End Do
046      j = j + k
047    End Do
048    Return
049  End Subroutine fcFFT
050 
051End Module FFT_Mod
052   
053Module typedef
054  Integer , parameter :: SP   = Selected_Real_Kind( p = 5 )
055  Integer , parameter :: I4b  = Selected_Int_Kind( 5 )
056End Module typedef
057 
058Function convlv(x, respns, isign)
059  use typedef
060  use FFT_Mod
061  Implicit None
062  Real (Kind=sp), Dimension (:), Intent (Inout) :: x
063  Real (Kind=sp), Dimension (:), Intent (In) :: respns
064  Integer (Kind=i4b), Intent (In) :: isign
065  Real (Kind=sp), Dimension (size(x)) :: convlv
066  Integer (Kind=i4b) :: n, m
067  Complex (Kind=SP), Dimension (size(x)) :: tmpd, tmpr
068  tmpd = 0.0_sp
069  tmpr = 0.0_sp
070  n = size(x)
071  m = size(respns)
072  tmpd = cmplx( x , 0.0_sp )
073  call fcFFT( tmpd , FFT_Forward )
074  tmpr = cmplx( respns , 0.0_sp )
075  call fcFFT( tmpr , FFT_Forward )
076  If (isign==1) Then
077    tmpr(1) = cmplx(real(tmpd(1))*real(tmpr(1))/n, aimag(tmpd(1))*aimag(tmpr(1))/n, kind=sp)
078    tmpr(2:) = tmpd(2:)*tmpr(2:)/n
079  Else If (isign==-1) Then
080    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'
081    tmpr(1) = cmplx(real(tmpd(1))/real(tmpr(1))/n, aimag(tmpd(1))/aimag(tmpr(1))/n, kind=sp)
082    tmpr(2:) = tmpd(2:)/tmpr(2:)/n
083  Else
084    write(*,*) 'no meaning for isign in convlv'
085  End If
086  Call fcFFT( tmpr , FFT_Inverse )
087  convlv = real( tmpr ) 
088End Function convlv
089 
090Program conv_main_fft
091  use typedef
092  Implicit None
093   
094  Interface
095    FUNCTION convlv(data,respns,isign)
096      use typedef
097      REAL(SP), DIMENSION(:), INTENT(INOUT) :: data
098      REAL(SP), DIMENSION(:), INTENT(IN) :: respns
099      INTEGER(I4B), INTENT(IN) :: isign
100      REAL(SP), DIMENSION(size(data)) :: convlv
101    End Function convlv
102  End Interface
103   
104  Integer , Parameter :: n=16, m=16
105  Real(Kind=SP) :: data1(2*n), respns(m), tmpr(2*n)
106  integer :: i , isign
107  data1(:) = 0.0_SP
108  data1( (n/2-n/8) : (n/2+n/8) ) = 1.0_SP
109  respns(:) = 0.0_SP
110  respns( 3:6 ) = 1.0_SP
111  isign = 1
112  tmpr = 0.0_SP
113  tmpr = convlv(data1, respns, isign)
114  Write (*, '(9999(f5.1,1x))' ) tmpr
115End Program conv_main_fft

[通用算法] 求助，卷积的小例子

新人勋章

爱心勋章

水王勋章

元老勋章

热心勋章