用FFT实现相关函数快速估计

VibInfo

主程序：

1. c----------------------------------------------------------------------
   
2. c   Main program HCORRE2 : to test subroutine CORRE2
   
3. c   To compute the auto-correlation function of a simple data by FFT
   
4. C   Please link subroutine CORRE2,SPLFFT
   
5. c----------------------------------------------------------------------
   
6.         complex x(0:15)
   
7.         data m/8/,n/16/,icorre/0/
   
8.         do 10 i=0,m-1
   
9.            x(i)=i+1.
   
10. 10      continue
    
11.         call corre2(x,x,m,n,icorre,ierror)
    
12.         write(*,*)'   ierror=',ierror
    
13.         if(ierror.ne.0) stop
    
14.         write(*,*)(x(i),i=0,n-1)
    
15.         stop
    
16.         end

复制代码

子程序：

1.         subroutine corre2(x,y,m,n,icorre,ierror)
   
2. c----------------------------------------------------------------------
   
3. c  Routinue CORRE2: To Compute the Correlation  Function of x(i) and
   
4. c  y(i) by DFT. x(i),y(i),i=0,...,m-1;But the dimension n of x,y must
   
5. c  be >=2*m and be the power of 2 ;
   
6. c  If: icorre=0:For auto-correlation
   
7. c      icorre=1:For crros-correlation
   
8. c                                       in chapter 11
   
9. c----------------------------------------------------------------------
   
10.         complex x(0:n-1),y(0:n-1)
    
11.         ierror=1
    
12.         s=n
    
13. 1       s=s/2.
    
14.         if(s-2.)70,5,1
    
15. 5       do 10 i=m,n-1
    
16. 10          x(i)=0.
    
17.         call splfft(x,n,-1)
    
18.         if(icorre.eq.1)goto 30
    
19.         do 20 k=0,n-1
    
20.            x(k)=conjg(x(k))*x(k)/float(m)
    
21. 20      continue
    
22.         goto 60
    
23. 30      do 40 i=m,n-1
    
24. 40         y(i)=0.
    
25.         call splfft(y,n,-1)
    
26.         do 50 k=0,n-1
    
27.            x(k)=conjg(x(k))*y(k)/float(m)
    
28. 50      continue
    
29. 60      call splfft(x,n,1)
    
30.         ierror=0
    
31. 70      return
    
32.         end

复制代码