由AR模型参数得到功率谱
subroutine ar1psd(a,ip,mfre,ep,ts)c-----------------------------------------------------------------------
c Routine AR1PSD: To compute the power spectum by AR-model parameters.
c Input parameters:
c IP : AR model order (integer)
c EP : White noise variance of model input (real)
c Ts : Sampling interval in seconds (real)
c A : Complex array of ARparameters, A(0) to A(IP)
c Output parameters:
c PSDR : Real array of power spectral density values
c PSDI : Real work array
c in chapter 12
c-----------------------------------------------------------------------
complex a(0:ip)
dimension psdr(0:4095),psdi(0:4095)
do 10 k=0,ip
psdr(k)=real(a(k))
psdi(k)=aimag(a(k))
10 continue
do 20 k=ip+1,mfre-1
psdr(k)=0.
psdi(k)=0.
20 continue
call relfft(psdr,psdi,mfre,-1)
do 30 k=0,mfre-1
p=psdr(k)**2+psdi(k)**2
psdr(k)=ep*ts/p
30 continue
call psplot(psdr,psdi,mfre,ts)
return
end PS:
relfft(psdr,psdi,mfre,-1)
psplot(psdr,psdi,mfre,ts)
是两个子程序吧，请问我可以在哪里找到？ subroutine relfft(xr,xi,n,isign)
c----------------------------------------------------------------------
croutinue RELFFT: To performthe split-radix DIF fft algorithm.
cinput parameters:
c xr,xi:real and image part of complex data for DFT/IDFT,n=0,...,N-1
c n : the dimension of xr and xi ;
cISIGN=-1: For Forward Transform;
c =+1: For Inverse Transform.
coutput parameters:
c xr,xi:real and image partof DFT/IDFT's result,n=0,...,N-1
cNote: Nmust be a power of 2 .
cFrom Ref. of chapter 5, Modified by Lao Changan
c in chapter 5
c----------------------------------------------------------------------
dimension xr(0:n-1),xi(0:n-1)
do 5 m=1,16
n2=2**m
if(n.eq.n2) go to 8
5 continue
write(*,*)'N is not a power of 2'
return
8 n2=n*2
es=-isign*atan(1.0)*8.0
do 10 k=1,m-1
n2=n2/2
n4=n2/4
e=es/n2
a=0
do 20 j=0,n4-1
a3=3*a
cc1=cos(a)
ss1=sin(a)
cc3=cos(a3)
ss3=sin(a3)
a=(j+1)*e
is=j
id=2*n2
40 do 30 i0=is,n-1,id
i1=i0+n4
i2=i1+n4
i3=i2+n4
r1=xr(i0)-xr(i2)
s1=xi(i0)-xi(i2)
r2=xr(i1)-xr(i3)
s2=xi(i1)-xi(i3)
xr(i0)=xr(i0)+xr(i2)
xi(i0)=xi(i0)+xi(i2)
xr(i1)=xr(i1)+xr(i3)
xi(i1)=xi(i1)+xi(i3)
if(isign.eq.1)goto 41
s3=r1-s2
r1=r1+s2
s2=r2-s1
r2=r2+s1
goto 42
41 s3=r1+s2
r1=r1-s2
s2=-r2-s1
r2=-r2+s1
42 xr(i2)=r1*cc1-s2*ss1
xi(i2)=-s2*cc1-r1*ss1
xr(i3)=s3*cc3+r2*ss3
xi(i3)=r2*cc3-s3*ss3
30 continue
is=2*id-n2+j
id=4*id
if(is.lt.n-1) goto 40
20 continue
10 continue
c
c ------------ special last stage -------------------------
is=0
id=4
50 do 60 i0=is,n-1,id
i1=i0+1
xtr=xr(i0)
xti=xi(i0)
xr(i0)=xtr+xr(i1)
xi(i0)=xti+xi(i1)
xr(i1)=xtr-xr(i1)
xi(i1)=xti-xi(i1)
60 continue
is=2*id-2
id=4*id
if(is.lt.n-1) goto 50
100 j=1
n1=n-1
do 104 i=1,n1
if(i.ge.j) goto 101
xtr=xr(j-1)
xti=xi(j-1)
xr(j-1)=xr(i-1)
xi(j-1)=xi(i-1)
xr(i-1)=xtr
xi(i-1)=xti
101 k=n/2
102 if(k.ge.j) goto 103
j=j-k
k=k/2
goto 102
103 j=j+k
104 continue
if(isign.eq.-1)return
do 105 i=0,n-1
xr(i)=xr(i)/n
xi(i)=xi(i)/n
105 continue
return
end subroutine psplot(psdr,psdi,mfre,ts)
c-----------------------------------------------------------------------
c Routine PSPLOT: To plot the normalized power spectum curve on the
c normalized frequency axis from -.5 to+.5 .
c mfre : Points in frequency axis and must be the power of 2.
c ts : Sampling interval in seconds (real).
c psdr : Real array of power spectral density values.
c psdi : Real work array.
c in Chapter 11
c-----------------------------------------------------------------------
dimension psdr(0:mfre-1),psdi(0:mfre-1)
m2=mfre/2
do 40 k=0,m2-1
psdi(k)=psdr(k)
psdr(k)=psdr(k+m2)
psdr(k+m2)=psdi(k)
40 continue
pmax=psdr(0)
do 50 k=1,mfre-1
if(psdr(k).gt.pmax) pmax=psdr(k)
50 continue
do 60 k=0,mfre-1
psdr(k)=psdr(k)/pmax
if(psdr(k).le.0.0)psdr(k)=.000001
60 continue
fs=1./ts
fs=fs/float(mfre)
open(3,file='psd.dat',status='new')
do 100 k=0,mfre-1
faxis=fs*(k-m2)
write(3,*)faxis,10.*alog10(psdr(k))
100 continue
close(3)
return
end 原帖由 dxjuan115 于 2006-11-13 16:00 发表
PS:
relfft(psdr,psdi,mfre,-1)
psplot(psdr,psdi,mfre,ts)
是两个子程序吧，请问我可以在哪里找到？
你看一下是否是上面给的两个程序 收藏了，谢谢 好东西，收藏，谢谢
页:
[1]