关于谐波滤波的简单原理与方法，还有fortran程序。

lqouc · 发表于 2013-1-17 21:13:54

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 lqouc 于 2013-1-17 22:55 编辑

之前有一次发过一个帖子，把自己编的关于空间谐波滤波的程序拿出来提问，结果有些朋友（@大果粒）想要了解这种滤波方法，我就在这里小小的总结一下，而且发现家园和资料站都没有关于谐波滤波的程序和原理，就顺便把程序也发上来。这些都是个人自己的思考，没有经过高水平人士的验证，不保证正确，而且程序也并不是很通用，就当时给初学者的一块敲门砖吧。

关于谐波滤波的目的，首先我们的气象数据有空间两维（水平）和时间一维，根据波动理论我们知道很多的空间波动和时间波动都包含在数据中，但是在分析的时候我们并不需要所有的信息，只要选取最重要的一部分就可以了。例如空间波动在对流层中对我们影响最大的是行星尺度波（纬圈0-3波）和天气尺度波（纬圈6-12波）。这样为了得到某一个波段的波动，我们需要进行滤波。（对于气象来说时间和空间是相对应的）

进行滤波的手段有很多，比如平均，滑动平均构造的滤波器，谐波滤波，甚至EOF方法。以位势高度来讨论，a：最简单的是平均，例如对数据进行月平均，这样就去除了要素中的快变成分，得到的数据不含有天气尺度波和中尺度波。b：滑动平均稍微高级一点，可以用来构造带通滤波器，保留什么样的波段就看你取多少点进行滑动平均了。c：谐波滤波在下面会仔细讲一下。d：EOF其实也是可以做到滤波的，但是本质上讲不是用来滤波，它是把不同频率的波动按照方差贡献排列下来，通过不同模的时间系数来反映的，这样的好处就是，即便我们不清楚我们的数据中什么频率的波动影响最强，它也可以为你找到这个频率，并提取出来（而且也同时得到了时间变化哦）。

这里就是主要内容了，谐波滤波：

时域上的Fourier变换产生频率域上的谱分量。

空间域上的Fourier变换产生波数域上的谱分量。

谐波分析的基本原理，是将任意已知周期T的时间序列看作是由许多频率（周期）与基波频率（周期）成倍数的谐波叠加，基本周期等于序列的长度，即T=n，称为基波，其余各波称为谐波，他们的频率fk（k=1,2,3……，p）分别为基波频率的2,3,4，……，p倍，而它们的频率则分别为基波的1/2,1/3,1/4……,1/p（p为谐波数）。因此，最短的周期为观测间隔的2倍。如果序列长度为n=24，基本周期就是T=24，其他各谐波周期就是24/2,24/3,…….,24/p。根据上述原理，具有任意基本周期T的时间（以时间为例）序列xt可以由下式给出：

（1）

A0是序列的算术平均值， QQ截图20130117224801.png

是各谐波的振幅

是各谐波的波向角。每个谐波的周期为 QQ截图20130117224845.png

，称为基波周期，而 w为基频。对上面的式子做变换，有了

（2）

令

所以（2）可以写为

（3）

谐波分析的目的，就是要确定系数 QQ截图20130117225024.png

，以便利用k个谐波来逼近原始序列。根据最小二乘原理和三角函数的正交性，可以得到

其中称为傅氏系数

当然fortran是不会积分的，我们把上面的式子改写为有限求和形式，则有

这样我们就可以根据上面的原理和公式进行编程和滤波了。

!这是一个空间滤波的程序，在这个程序里面我们利用ncep的2.5*2.5*17的日平均位势高度数据
!目的是对位势高度的空间波动进行空间谐波滤波，得到空间行星尺度波和天气尺度波
!当然这个程序并不是很好，权当是作为一个参考吧，我也是菜鸟。
!强调一下，所谓的空间滤波只是对纬向方向的空间滤波，因为自由大气的平均状态是纬向风。
program lbx
implicit none
real,dimension(1:144,1:73,1:17,1:365)::hgt !这里是数据格式，x，y，z，t
real,dimension(0:12,1:73,12:17,1:365)::a,b !这个是傅氏系数,在这里只滤出0到12波
real,dimension(1:144,1:73,12:17,1:365)::lb0t3,lb6t12 !这两个是行星尺度和天气尺度波动，是多个频率波动的叠加
real,dimension(0:12,1:144,1:73,12:17,1:365)::lb
real,parameter::w=2.5 !这个是一个空间格距，通过这个可以计算数据出在空间中的位置
integer,parameter::it=365 !这里设置总时间长度，可以根据自己的需要更改，在这里我就算一年的
integer::i,j,k,t,m,irec !这里一般分别代表x，y，z，t，波数，换行
!第一步，读取数据
open(11,file='hgt1998.dat',status='old',form='unformatted',access='direct',recl=144*73)
irec=1
do t=1,it
do k=1,17
read(11,rec=irec)((hgt(i,j,k,t),i=1,144),j=1,73)
irec=irec+1
end do
end do
close(11)
print*,'*'
!这就开始空间滤波了
!第二步，对所有谐波的系数进行初始设置
a=0
b=0
print*,'**'
!第三步，计算每个谐波的系数（这里其实只算了上文最后有限求和公式的求和部分）
do t=1,it
do k=12,17
do m=0,12
do j=1,73
do i=1,144
a(m,j,k,t)=a(m,j,k,t)+hgt(i,j,k,t)*cosd(m*w*i)
b(m,j,k,t)=b(m,j,k,t)+hgt(i,j,k,t)*sind(m*w*i)
end do
end do
end do
end do
end do
!再除以空间序列的长度（这里是上边没完成的平均）
a=a/72
b=b/72
do t=1,it
do k=12,17
do j=1,73
a(0,j,k,t)=a(0,j,k,t)/2
end do
end do
end do
print*,'***'
!第四步，累加各个谐波，得到各个波段的空间滤波场（这里就是得到了公式3中的某一个频率的谐波和基波的叠加，和公式并不一样的是这里只有一个谐波）
do t=1,it
do k=12,17
do m=1,12
do j=1,73
do i=1,144
lb(m,i,j,k,t)=a(0,j,k,t)+a(m,j,k,t)*cosd(m*w*i)+b(m,j,k,t)*sind(m*w*i)
end do
end do
end do
end do
end do
!得到平均场即0波滤波
do t=1,it
do k=12,17
do j=1,73
do i=1,144
lb(0,i,j,k,t)=a(0,j,k,t)
end do
end do
end do
end do
!得到0-3波和6-12波的空间滤波
do t=1,it
do k=12,17
do j=1,73
do i=1,144
lb0t3(i,j,k,t)=lb(1,i,j,k,t)+lb(2,i,j,k,t)+lb(3,i,j,k,t)-2*lb(0,i,j,k,t)
lb6t12(i,j,k,t)=lb(6,i,j,k,t)+lb(7,i,j,k,t)+lb(8,i,j,k,t)+lb(9,i,j,k,t)+lb(10,i,j,k,t)+lb(11,i,j,k,t)+lb(12,i,j,k,t)-6*lb(0,i,j,k,t)
end do
end do
end do
end do
print*,'****'
!第五步，结果输出（这里就是按照GrADS格式写的，可以直接写个ctl画图了）
open(11,file='lb-all1998.dat',status='replace',form='unformatted',access='direct',recl=144*73)
irec=1
do t=1,it
do m=0,12
do k=12,17
write(11,rec=irec)((lb(m,i,j,k,t),i=1,144),j=1,73)
irec=irec+1
end do
end do
do k=12,17
write(11,rec=irec)((lb0t3(i,j,k,t),i=1,144),j=1,73)
irec=irec+1
end do
do k=12,17
write(11,rec=irec)((lb6t12(i,j,k,t),i=1,144),j=1,73)
irec=irec+1
end do
end do
close(11)
print*,'ok!'
!完成
end program

!这是一个空间滤波的程序，在这个程序里面我们利用ncep的2.5*2.5*17的日平均位势高度数据 
!目的是对位势高度的空间波动进行空间谐波滤波，得到空间行星尺度波和天气尺度波 
!当然这个程序并不是很好，权当是作为一个参考吧，我也是菜鸟。 
!强调一下，所谓的空间滤波只是对纬向方向的空间滤波，因为自由大气的平均状态是纬向风。 
program lbx 
implicit none 
real,dimension(1:144,1:73,1:17,1:365)::hgt  !这里是数据格式，x，y，z，t 
real,dimension(0:12,1:73,12:17,1:365)::a,b  !这个是傅氏系数,在这里只滤出0到12波 
real,dimension(1:144,1:73,12:17,1:365)::lb0t3,lb6t12   !这两个是行星尺度和天气尺度波动，是多个频率波动的叠加 
real,dimension(0:12,1:144,1:73,12:17,1:365)::lb 
real,parameter::w=2.5                       !这个是一个空间格距，通过这个可以计算数据出在空间中的位置 
integer,parameter::it=365                   !这里设置总时间长度，可以根据自己的需要更改，在这里我就算一年的 
integer::i,j,k,t,m,irec                     !这里一般分别代表x，y，z，t，波数，换行 
!第一步，读取数据 
open(11,file='hgt1998.dat',status='old',form='unformatted',access='direct',recl=144*73) 
 irec=1 
 do t=1,it 
  do k=1,17 
   read(11,rec=irec)((hgt(i,j,k,t),i=1,144),j=1,73) 
   irec=irec+1 
  end do 
 end do 
close(11) 
print*,'*' 
!这就开始空间滤波了 
 
!第二步，对所有谐波的系数进行初始设置 
a=0 
b=0 
print*,'**' 
!第三步，计算每个谐波的系数（这里其实只算了上文最后有限求和公式的求和部分） 
do t=1,it 
do k=12,17 
 do m=0,12 
  do j=1,73 
   do i=1,144 
    a(m,j,k,t)=a(m,j,k,t)+hgt(i,j,k,t)*cosd(m*w*i) 
    b(m,j,k,t)=b(m,j,k,t)+hgt(i,j,k,t)*sind(m*w*i) 
   end do 
  end do 
 end do 
end do 
end do 
!再除以空间序列的长度（这里是上边没完成的平均） 
a=a/72 
b=b/72 
do t=1,it 
do k=12,17 
 do j=1,73 
  a(0,j,k,t)=a(0,j,k,t)/2 
 end do 
end do 
end do 
print*,'***' 
!第四步，累加各个谐波，得到各个波段的空间滤波场（这里就是得到了公式3中的某一个频率的谐波和基波的叠加，和公式并不一样的是这里只有一个谐波） 
do t=1,it 
do k=12,17 
 do m=1,12 
  do j=1,73 
   do i=1,144 
   lb(m,i,j,k,t)=a(0,j,k,t)+a(m,j,k,t)*cosd(m*w*i)+b(m,j,k,t)*sind(m*w*i) 
   end do 
  end do 
 end do 
end do 
end do 
 
!得到平均场即0波滤波 
do t=1,it 
do k=12,17 
 do j=1,73 
  do i=1,144 
   lb(0,i,j,k,t)=a(0,j,k,t) 
  end do 
 end do 
end do 
end do 
 
!得到0-3波和6-12波的空间滤波 
do t=1,it 
do k=12,17 
 do j=1,73 
  do i=1,144 
   lb0t3(i,j,k,t)=lb(1,i,j,k,t)+lb(2,i,j,k,t)+lb(3,i,j,k,t)-2*lb(0,i,j,k,t) 
   lb6t12(i,j,k,t)=lb(6,i,j,k,t)+lb(7,i,j,k,t)+lb(8,i,j,k,t)+lb(9,i,j,k,t)+lb(10,i,j,k,t)+lb(11,i,j,k,t)+lb(12,i,j,k,t)-6*lb(0,i,j,k,t) 
  end do 
 end do 
end do 
end do 
print*,'****' 
!第五步，结果输出（这里就是按照GrADS格式写的，可以直接写个ctl画图了） 
open(11,file='lb-all1998.dat',status='replace',form='unformatted',access='direct',recl=144*73) 
irec=1 
do t=1,it 
 do m=0,12 
  do k=12,17 
  write(11,rec=irec)((lb(m,i,j,k,t),i=1,144),j=1,73) 
  irec=irec+1 
  end do 
 end do 
 do k=12,17 
  write(11,rec=irec)((lb0t3(i,j,k,t),i=1,144),j=1,73) 
  irec=irec+1 
 end do 
 do k=12,17 
  write(11,rec=irec)((lb6t12(i,j,k,t),i=1,144),j=1,73) 
  irec=irec+1 
 end do 
end do 
close(11) 
 
print*,'ok!' 
!完成 
 
end program

这样就完成了，顺便附上参考资料：胡基福,《气象统计原理与方法》,青岛海洋大学出版社,青岛,1996
还有我们学校老师上课时给的两个例子的图，给大家看看效果
无标题.png

这是总场，行星波和天气波的。
啰啰嗦嗦写了这么多，主要是之前自己做这个的时候，找遍了网上也没有现成的程序，花了两天的时间弄这个，送给需要急用的朋友吧，但是建议大家还是自己好好的学学原理。另外程序编的不好，我不很擅长这个，希望高人能给改进下，最好写个字程序出来，我实在是懒得写了。

ayzqs · 发表于 2013-1-17 21:34:20

真给力！

大果粒 · 发表于 2013-1-17 21:34:32

谢谢楼主，不过我根据你的程序算出来的形势是对的，可是值好像大好多倍，我试了差不多2-3倍，这个是什么问题呢？

lqouc · 发表于 2013-1-17 21:51:30

大果粒发表于 2013-1-17 21:34

登录/注册后可看大图

谢谢楼主，不过我根据你的程序算出来的形势是对的，可是值好像大好多倍，我试了差不多2-3倍，这个是什么问题 ...

这个我也不是很清楚啊，我做的都还可以啊，这个程序我确实也没有验证过。你要的是单独某个谐波的么？如果是的话我的程序计算出来的都是叠加平均场的谐波。或者不妨说说你的数据使用的什么，我们一起改进下。

lqouc · 发表于 2013-1-17 21:52:44

@言深深求版主指点一下啊，我辛辛苦苦写了很久的共识怎么就贴不上去呢？这个没公式没法看啊。

大果粒 · 发表于 2013-1-17 21:52:50

还有楼主的公式貌似没有上传成功呢，谢谢楼主~

大果粒 · 发表于 2013-1-17 22:09:38

lqouc 发表于 2013-1-17 21:51

登录/注册后可看大图

这个我也不是很清楚啊，我做的都还可以啊，这个程序我确实也没有验证过。你要的是单独某个谐波的么？如果 ...

我就是用gfs资料，网格1*1的，按你的这样算算下来还有数值有问题的，也是算0-3波的值，要是看到公式就能知道问题在哪里了应该，不知道有没有文献出处？

言深深 · 发表于 2013-1-17 22:25:52

lqouc 发表于 2013-1-17 21:52

登录/注册后可看大图

@言深深求版主指点一下啊，我辛辛苦苦写了很久的共识怎么就贴不上去呢？这个没公式没法看啊。

qq截图，做成图片，然后贴图片就ok了···论坛本身不支持公式呢······

lqouc · 发表于 2013-1-17 22:59:42

大果粒发表于 2013-1-17 22:09

登录/注册后可看大图

我就是用gfs资料，网格1*1的，按你的这样算算下来还有数值有问题的，也是算0-3波的值，要是看到公式就能知 ...

现在有公式了，不过这个水印打得我好悲伤啊，凑活看吧，根据上下文合理推断一下。至于文献上面写了啊。胡基福的气象统计。

爱的侍者 · 发表于 2013-1-17 23:12:23

为啥层数是12:17层啊？？
解释下

还有3——6 的波怎么办
你想提取行星尺度波和天气尺度波，也不包括这个

		自动登录	找回密码
密码			立即注册

[经验总结] 关于谐波滤波的简单原理与方法，还有fortran程序。

登录后查看更多精彩内容~

评分