1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
|
C/MEMBR ADD NAME=PROJ2,SSI=0
c Copyright INRIA
subroutine proj2(f,nn,am,n,np1,np2,pf,w)
c ce sous programme calcule les produits scalaires:
c <f,(z**(n-v)/am>=pf(np+1-v) v=np1...np2 ;np=np2-np1+1
c les pf(v) sont les dernieres valeurs de la filtree de f(nn+1-i) par
c (z**n-1)/am,am est un polynome de degre n range par puissance
c croissante
c
c w:tableau de travail de taille n
c
implicit double precision (a-h,o-z)
dimension f(nn),am(*),pf(*),w(n)
c
c w contient l'etat du filtre
c
np=np2
if(np1.gt.1) np=np2-np1+1
c
nf=0
nb=nn-np
an=am(n+1)
if(n.eq.1) goto 50
call dset(n,0.0d+0,w,1)
do 20 i=1,nb
wn=w(n)/an
nf=nf+1
iw=n+1
do 10 j=2,n
iw=iw-1
w(iw)=w(iw-1)-am(iw)*wn
10 continue
w(1)=-am(1)*wn
w(n)=w(n)+f(nf)
20 continue
c les n valeurs suivantes de la sortie du filtre donnent les produits
c scalaires
do 31 i=1,np
wn=w(n)/an
iw=n+1
nf=nf+1
do 30 j=2,n
iw=iw-1
w(iw)=w(iw-1)-am(iw)*wn
30 continue
w(1)=-am(1)*wn
w(n)=w(n)+f(nf)
pf(i)=w(n)/an
31 continue
if(np1.ge.1) return
do 41 i=np+1,np2-np1+1
wn=w(n)/an
iw=n+1
do 40 j=2,n
iw=iw-1
w(iw)=w(iw-1)-am(iw)*wn
40 continue
w(1)=-am(1)*wn
pf(i)=w(n)/an
41 continue
return
c cas particulier n=1
50 wn=0.0d+0
do 60 i=1,nb
nf=nf+1
wn=-am(1)/an*wn+f(nf)
60 continue
do 70 i=1,np
nf=nf+1
wn=-am(1)/an*wn+f(nf)
pf(i)=wn/an
70 continue
if(np1.ge.1) return
do 71 i=np+1,np2-np1+1
nf=nf+1
wn=-am(1)/an*wn
pf(i)=wn/an
71 continue
return
end
|