1C/MEMBR ADD NAME=CHEBY,SSI=0 2 subroutine cheby(ordr,demi,ieo,dp,x0,tam,win) 3c! 4c sous-programme: cheby 5c fenetre de dolph chebyshev 6c en double precision 7c acheve le 05/12/85 8c ecrit par philippe touron 9c 10c 11c parametres entrants 12c ------------------- 13c ordr, l'ordre du filtre (entier) 14c demi, l'ordre de la demi fenetre (entier) 15c ieo, indicateur de parite vaut 0 si ordr pair et 1 sinon (entier) 16c dp, attenuation en absolu (reelle) 17c x0, constante de la fenetre de chebyshev fonction de 18c la largeur du lobe principale (reelle) 19c tam, tableau de travail (reels) qui doit etre dans un programme 20c appelant dimensionne a 3 fois ordr) 21c 22c parametres sortants 23c ------------------- 24c win, tableau de valeurs de la fenetre (reelles) qui doit dans 25c un programme appelant etre dimensionne a ordr 26c 27c variables internes 28c ------------------ 29c c0,c1,c2, des buffers de calcul (reels) 30c compt,indi,xcompt,xindi, compteurs de boucles (entiers et reels) 31c teta,gama, coefficients de la fenetres de chebyshev (reels) 32c 33c sous programmes appele: coshin 34c 35c! 36 double precision win(*),tam(*) 37 double precision x0,dp,xordr,teta,gama,freq 38 double precision pi,twopi,twn,c0,c1,c2,somm,xcompt,xindi 39 double precision coshin 40 integer ordr,demi,compt,indi,icompt,jcompt 41c 42 xordr=dble(ordr) 43 teta=(x0+1.)/2.0d+0 44 gama=(x0-1.)/2.0d+0 45 pi=acos(-1.0d+0) 46 twopi=2.0d+0*pi 47 c0=(xordr-1.0d+0)/2.0d+0 48 do 40 compt=1,ordr 49 icompt=ordr+compt 50 jcompt=icompt+ordr 51 xcompt=dble(compt)-1.0d+0 52 c1=xcompt/xordr 53 freq=teta*cos(twopi*c1)+gama 54 CRES=abs(freq)-1.0d+0 55 if (CRES .le. 0) then 56 goto 10 57 else 58 goto 20 59 endif 6010 tam(icompt)=dp*cos(c0*acos(freq)) 61 goto 30 6220 tam(icompt)=dp*cosh(c0*coshin(freq)) 6330 tam(jcompt)=0.0d+0 64c 65c modification si filtre d'ordr pair 66c 67 if(ieo.eq.1)goto 40 68 tam(jcompt)=-tam(icompt)*sin(pi*c1) 69 tam(icompt)= tam(icompt)*cos(pi*c1) 70 if(compt.gt.(ordr/2+1))then 71 tam(icompt)=-tam(icompt) 72 tam(jcompt)=-tam(jcompt) 73 endif 7440 continue 75c 76c transformee de fourrier inverse 77c 78 twn=twopi/xordr 79 do 60 compt=1,demi 80 xcompt=dble(compt)-1.0d+0 81 somm=0.0d+0 82 do 50 indi=1,ordr 83 icompt=indi+ordr 84 jcompt=icompt+ordr 85 xindi=dble(indi)-1.0d+0 86 somm=somm+tam(icompt)*cos(twn*xindi*xcompt) 87 & +tam(jcompt)*sin(twn*xindi*xcompt) 8850 continue 89 win(compt)=somm 9060 continue 91 c2=win(1) 92 do 70 compt=1,demi 93 win(compt)=win(compt)/c2 9470 continue 95 return 96 end 97