Fortran Coder

查看: 10288|回复: 5

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

[复制链接]

12

帖子

4

主题

0

精华

入门

F 币
52 元
贡献
38 点
发表于 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] 纯文本查看 复制代码
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

convolution_1.rar

101.48 KB, 下载次数: 1

710

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
600 元
贡献
307 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 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] 纯文本查看 复制代码
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

12

帖子

4

主题

0

精华

入门

F 币
52 元
贡献
38 点
 楼主| 发表于 2015-8-16 15:05:02 | 显示全部楼层
恩恩,基本上是想做一维或者二维卷积,
程序还有一点小问题:错误 error #8000:  There is a conflict between local interface block and external interface block.   [CONVLV]       
两部分调试都出现了这种情况,不知道是不是我编译器的问题

710

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
600 元
贡献
307 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 2015-8-17 08:47:05 | 显示全部楼层
本帖最后由 楚香饭 于 2015-8-17 08:49 编辑

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

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

12

帖子

4

主题

0

精华

入门

F 币
52 元
贡献
38 点
 楼主| 发表于 2015-8-17 11:38:09 | 显示全部楼层
恩,使用CVF就没问题,谢谢群主,IVF我再看看我是不是哪设置出问题了

12

帖子

4

主题

0

精华

入门

F 币
52 元
贡献
38 点
 楼主| 发表于 2015-8-17 15:15:48 | 显示全部楼层
我把楚帅的代码中的外部函数都扔到模组里,然后删掉了接口块,在我的Intel版本里也能运行了,应该没啥问题,希望以后继续向论坛里的朋友学习,等姿势水平提高后,帮助更多热爱科学计算的朋友
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-4-16 08:59

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表