利用DFT的卷积性质求两个复序列的线性卷积

VibInfo

利用DFT的卷积性质求两个复序列的线性卷积

主程序：

1. C----------------------------------------------------------------------
   
2. C main program HCONVO2:To test subroutine CONVO2
   
3. c Please link subroutine CONVO2,SPLFFT
   
4. C----------------------------------------------------------------------
   
5.         complex h(0:7),x(0:7),y(0:7)
   
6.         data n1/4/,n2/2/,n/8/
   
7.         do 10 i=0,3
   
8.            h(i)=1.
   
9.            x(i)=i+1.
   
10. 10      continue
    
11.         h(0)=1.
    
12.         h(1)=1.
    
13.         x(0)=1.
    
14.         x(1)=2.
    
15.         x(2)=3.
    
16.         x(3)=4.
    
17.         call convo2(x,h,y,n1,n2,n)
    
18.         do 20 k=0,n1+n2-2
    
19. 20        write(*,*)k,real(y(k))
    
20.         stop
    
21.         end

复制代码

子程序：

1.         subroutine convo2(x,h,y,n1,n2,n)
   
2. c----------------------------------------------------------------------
   
3. c  Routine CONVO2: To Compute the Convolution of x(i) and h(i) by DFT.
   
4. c  x(i),i=0,...,n1-1; h(n),i=0,...,n2-1,But the dimension n of x,h,y
   
5. c       must be >=(n1+n2-1) and be the power of 2 ;
   
6. c input parameters:
   
7. c x( ):n dimensioned real array,signal data is stored in x(0) to x(n1-1).
   
8. c h( ):n dimensioned real array,impuls response is stored in h(0) to h(n2-1).
   
9. c output parameters:
   
10. c y( ):n dimensioned real array, y(i)=x(i)*h(i),i=0,...n-1.
    
11. c Notes:
    
12. c     n must be a power of 2.
    
13. c                                      in chapter 3
    
14. c----------------------------------------------------------------------
    
15.         complex x(0:n-1),y(0:n-1),h(0:n-1)
    
16.         do 10 i=1,16
    
17.            nn=2**i
    
18.            if(n.eq.nn) go to 20
    
19. 10       continue
    
20.         write(*,*)'  N is not a power of 2'
    
21.         return
    
22. 20      do 30 i=n1,n-1
    
23. 30          x(i)=0.
    
24.         do 40 i=n2,n-1
    
25. 40         h(i)=0.
    
26.         call splfft(x,n,-1)
    
27.         call splfft(h,n,-1)
    
28.         do 50 k=0,n-1
    
29. 50         y(k)=x(k)*h(k)
    
30.         call splfft(y,n,1)
    
31.         return
    
32.         end

复制代码

[ 本帖最后由 VibInfo 于 2006-8-11 07:17 编辑 ]