home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Shareware Overload
/
ShartewareOverload.cdr
/
progm
/
grafx.zip
/
SMOTH.FOR
< prev
next >
Wrap
Text File
|
1989-03-14
|
2KB
|
76 lines
c program main
c call gset ! enter graphic mode
c call smoth
c call aset ! return to alpha mode
c end
c -----------------------------------------------------------------------
subroutine smoth
c
c ...demo of display of time series data
c
c (C) Copyright 1988, 1989 by Jim Farrell All Rights Reserved.
c
parameter(n=256)
real x(n),y(n)
c
parameter (a2pi=6.283185307,rate=40.0)
c
irk=1
do 200 i=1,n
wt=a2pi*real(i-1)/rate
ty=10.0*sin(wt)+10.0*cos(2.0*wt)
x(i)=ty
200 continue
call putstr(10,20,' DEMONSTRATION OF TPLOT WINDOWS')
call putstr(12,20,' FOR A SIGNAL PROCESSING USAGE')
call pause('WHEN READY ') ! display message and pause
call gcls
xmn=0.0
xmx=real(n)/rate
call tplot(xmn,xmx,x,n,1,2) ! plot clean signal
call pause ('REVIEW PLOT OF PURE SIGNAL')
do 300 i=1,n
if(x(i).ge.0.0)then
x(i)=x(i)+5.0*rand(irk) ! add noise
else
x(i)=x(i)-5.0*rand(irk) ! add noise
endif
300 continue
call tplot(xmn,xmx,x,n,2,2) ! plot signal & noise
call pause ('REVIEW PLOT OF SIGNAL AND NOISE')
a=0.0
b=0.0
c=0.0
call exsmoth(x,n,0.12,a,b,c,y) ! smoothing function
call tplot(xmn,xmx,y,n,3,2) ! plot smoothed signal
call pause ('REVIEW PLOT OF SMOOTHED SIGNAL')
return
end
c --------------------------------------------------------------------
subroutine exsmoth(x,nx,al,a,b,c,xs)
c
c ...smooth a time series
c
c al - smoothing coefficient (0 < al < 1)
c
real x(*),xs(*)
c
if(a.eq.0.0.and.b.eq.0.0.and.c.eq.0.0)then ! compute coefficients
c=x(1)-2.0*x(2)+x(3)
b=x(2)-x(1)-1.5*c
a=x(1)-b-0.5*c
endif
be=1.0-al ! al - smoothing coefficient
be3=be**3
al2=al**2
al3=al*al2
do 300 i=1,nx
xs(i)=a+b+0.5*c
dif=xs(i)-x(i)
a=x(i)+be3*dif
b=b+c-1.5*al2*(2.0-al)*dif
c=c-al3*dif
300 continue
return
end