- 积分
- 754
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-1-13
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
以下是来自做个霸气的木头的分享,弱问想看年代际变化,这个平滑可以用吗?话说有没有7次平滑一说
!-----*----------------------------------------------------6---------7--
program ex
! M-term Guassian-Type Filter
! Input variables: n, x(n), m
! m: the term number used to running mean
! it must be an odd number.
! Output variables: y(n),z(n-m-1)
! y: the filtered series of x.
! Work parameters and array: c, cgm and ck(-(m-1)/2:(m-1)/2)
! c: a tunable parameter, generally, c>2.0.
! cgm: variance of Guassian distribution.
parameter(n=30,m=7) !n为数据长度,m为几点平滑,必须为奇数
dimension x(n),y(n)
dimension xw((-(m-1)/2+1):(n+(m-1)/2)),ck(-(m-1)/2:(m-1)/2) !work array
open(100,file='12.txt') !数据文件(单列数据)
open(200,file='13pot2.txt') !平滑后文件输出
do i=1,n
read(100,*) x(i)
enddo
undef=-99.99
pi=3.1415926
c=2.15
nl=(m-1)/2
cgm=float(nl)/c
x1=x(1)
xn=x(n)
if(x(1).eq.undef)x1=x(2)
if(x(n).eq.undef)xn=x(n-1)
do i=-nl+1,1
xw(i)=x1
enddo
do i=2,n-1
xw(i)=x(i)
enddo
do i=n,n+nl
xw(i)=xn
enddo
c1=1./(cgm*sqrt(2.*pi))
ck(0)=c1
do 10 i=1,nl
ck(i)=c1*exp(-i*i/(2.*cgm*cgm))
ck(-i)=ck(i)
10 continue
do 20 i=1,n
y(i)=0.
do j=-nl,nl
y(i)=y(i)+ck(j)*xw(i+j)
enddo
20 continue
if(x(1).eq.undef)y(1)=undef
if(x(n).eq.undef)y(n)=undef
do i=1,n
write(200,'(f10.2)') y(i)
enddo
end
|
|