谐波分析的位相问题

南气斑织逍遥

如图所示，我已经得到了每个站点24小时平均的（24hours,lat,lon）的数组，下面文中直接说进行谐波分析，图中给出的横轴是反映的时间，那么如何进行这个谐波分析呢？
我查阅了一些资料，找到了傅里叶分解的ncl公式，但是这个一般是用来做纬向的空间分解，看一波二波什么的，函数是“ezfftf”“ezfftb”，想求助像图中给出的这种谐波分析的位相分布是个什么算法？

江奈何

用谐波分析计算位相的话，应该是取第一个 谐波，先用ezfftf进行傅里叶分级，再用ezfftb提取第一谐波的实部和虚部系数，那位相角就可以用atan(cImag/cReal)来计算，你再把位相角转化为相应的小时就可以了

南气斑织逍遥

江奈何 发表于 2019-9-19 08:32
用谐波分析计算位相的话，应该是取第一个 谐波，先用ezfftf进行傅里叶分级，再用ezfftb提取第一谐波的实部 ...

你好！我现在部分代码如下
z           = pr_ano(lat|:,lon|:,hour|:)               ；40lats*100lons*24hours
z_wave2 = dim_avg_n_Wrap(z,2)                    ;lat,lon(最后想把时间值返还到这个数组中)
z_wave2 = 0.

  CF      = ezfftf(z)                       ; ezfftf works on right dim
  CF!1    = "lat"
  CF&lat  = pr&LAT81_160
  CF!2    = "lon"
  CF&lon  = pr&LON441_720
  printVarSummary(CF)                 ; [2] x [40] x [100]*12

  cf      = CF
  cf(:,:,:,1:)  = 0.0                        ; set waves 2 and higher to 0.0 (选取1波)
；  z_wave1 = ezfftb (cf, 0.0)       ；这一步做出来的话就重新变成了 40*100*24的数组，我这里是否需要？
；  copy_VarMeta(z22,z_wave1)
  do i = 0,39                       ；这后边是我猜着写的，最后出的图和原图看上去差不少......
    do j = 0,99
       cImag = CF(1,i,j,:)       ;虚部
       cReal = CF(0,i,j,:)        ;实部
       weixiangjiao = atan(cImag/cReal)
       dims = dimsizes(weixiangjiao)
       inds = ind_resolve(maxind(weixiangjiao),dims) ;0,0
       pppp = inds(0,0)
       z_wave2(i,j) = pppp*2
       delete(weixiangjiao)
       delete(inds)
       delete(pppp)
    end do
  end do
能否再指点指点是哪里出了问题？

南气斑织逍遥

如图................
这个位相角需要怎么才能返回到时间呢？

江奈何

南气斑织逍遥 发表于 2019-9-19 11:06
你好！我现在部分代码如下
z           = pr_ano(lat|:,lon|:,hour|:)               ；40lats*100lons* ...

首先，ezfftb这个函数是表示傅里叶合成的，也就是在你提取完第一个人谐波后得到的数组，当然还是40*100*24，你只是把时间维做了滤波而已，这个你不用管它，你的问题是求位相；其次 ，在算实部和虚部的系数时，不需要循环，直接cimag=cf(1,:,:,0);  creal=cf(0,:,:,0); 第三个问题求得weixiangjiao后，将弧度转化为24小时，就是乘以24/(2pi) 不就好了

南气斑织逍遥

江奈何 发表于 2019-9-19 22:21
首先，ezfftb这个函数是表示傅里叶合成的，也就是在你提取完第一个人谐波后得到的数组，当然还是40*100*2 ...

首先感谢您的热心回复！
我马上操作了一波，现在代码如下：
  z           = pr_ano(lat|:,lon|:,hour|:)
  CF         = ezfftf(z)                          ; ezfftf works on right dim
  CF!1      = "lat"
  CF&lat   = pr&LAT81_160
  CF!2      = "lon"
  CF&lon  = pr&LON441_720
  printVarSummary(CF)                      ; [2] x [40] x [100]*12

  z_wave2 = dim_avg_n_Wrap(z,2) ;lat,lon
  z_wave2 = 0.
  cf      = CF
  cf(:,:,:,1:)  = 0.0                             ; set waves 2 and higher to 0.0 (1波)
  z_wave1 = ezfftb (cf, 0.0)   
  copy_VarMeta(z,z_wave1)
  cImag = cf(1,:,:,0)                           ;虚部 lat*lon
  cReal = cf(0,:,:,0)                            ;实部 lat*lon
  weixiangjiao = atan(cImag/cReal)      ;位相 lat*lon
  print(weixiangjiao)                           ；这里计算出的位相角是有正有负
  pi = 3.141592652
  hour  = weixiangjiao*12./pi              ；在此处计算的“时间场”的空间分布也是有正有负，原来从没接触过这方面的知识，真是搞不明白了
  hour!0    = "lat"
  hour&lat  = pr&LAT81_160
  hour!1    = "lon"
  hour&lon  = pr&LON441_720

江奈何

南气斑织逍遥 发表于 2019-9-19 22:45
首先感谢您的热心回复！
我马上操作了一波，现在代码如下：
  z           = pr_ano(lat|:,lon|:,hour| ...

好，现在我们假设一个XY轴，Y轴就代表cImag，X轴就代表cReal，分为四个象限，那计算xiaongweijiao，就转化为计算象限上（cReal,cImag）这个点的夹角了。那你有没有发现，1,3,象限和2,4象限计算的结果是一样的？所以，你在计算xiaongweijiao时要对creal和cimag分类，举个例子，若位于第二象限的点（-2,2），那atan（2/-2）=-4/pi,这当然不是你想要的结果，它被虚假地反应到了第四象限，所以你要加Pi,就变成了0.75pi, 在转化为时间，乘以12/pi, 不就是9个小时了吗，也就此时位相为9时。 以此类推

南气斑织逍遥

江奈何 发表于 2019-9-20 09:44
好，现在我们假设一个XY轴，Y轴就代表cImag，X轴就代表cReal，分为四个象限，那计算xiaongweijiao，就转 ...

按您的讲解将脚本改为如下：
  CF      = ezfftf(z)                 ; ezfftf works on right dim
  CF!1    = "lat"
  CF&lat  = pr&LAT81_160
  CF!2    = "lon"
  CF&lon  = pr&LON441_720
  printVarSummary(CF)                 ; [2] x [33] x [100]*12


  z_wave2       = dim_avg_n_Wrap(z,2) ;lat,lon
  z_wave2       = 0.
  cf            = CF
  cf(:,:,:,1:)  = 0.0                       ; set waves 2 and higher to 0.0 (1波)
  z_wave1 = ezfftb (cf, 0.0)   
  copy_VarMeta(z,z_wave1)
  cImag         = cf(1,:,:,0)            ;虚部
  cReal         = cf(0,:,:,0)             ;实部
  weixiangjiao  = z_wave2           ;lat,lon 位相角
  pi            = 3.141592652
  do i = 0,32  
    do j = 0,99
        if (.not.ismissing(cImag(i,j)) .and. .not.ismissing(cReal(i,j)) .and. cImag(i,j) .ge. 0.0 .and. cReal(i,j) .ge. 0.0) then        ;第一象限
        weixiangjiao(i,j) = atan(cImag(i,j)/cReal(i,j))
        else if(.not.ismissing(cImag(i,j)) .and. .not.ismissing(cReal(i,j)) .and. cImag(i,j) .gt. 0.0 .and. cReal(i,j) .lt. 0.0) then     ;第二象限
        weixiangjiao(i,j) = atan(cImag(i,j)/cReal(i,j)) + pi
        else if(.not.ismissing(cImag(i,j)) .and. .not.ismissing(cReal(i,j)) .and. cImag(i,j) .le. 0.0 .and. cReal(i,j) .le. 0.0) then     ;第三象限
        weixiangjiao(i,j) = atan(cImag(i,j)/cReal(i,j)) + pi
        else if(.not.ismissing(cImag(i,j)) .and. .not.ismissing(cReal(i,j)) .and. cImag(i,j) .lt. 0.0 .and. cReal(i,j) .gt. 0.0) then     ;第四象限
        weixiangjiao(i,j) = atan(cImag(i,j)/cReal(i,j)) + pi*2.  
        end if
        end if
        end if
        end if
    end do
  end do
  hour  = weixiangjiao*12./pi
  hour!0    = "lat"
  hour&lat  = pr&LAT81_160
  hour!1    = "lon"
  hour&lon  = pr&LON441_720


出图如下：

感觉现在的图和文中的存在相似特征了，下面我再尝试尝试，非常感谢您的帮助！！！希望这个帖子可以帮助更多人

江奈何

南气斑织逍遥 发表于 2019-9-20 11:19
按您的讲解将脚本改为如下：
  CF      = ezfftf(z)                 ; ezfftf works on right dim
  C ...

我也是之前遇到过求位相角的问题，比你这个还复杂我还要求相差，只是我想不出更好的方法了，所以只能分类来讨论，如果楼主有更好的方法记得私我哦；   另外，你这个第四象限是不是应该+3/2pi会好些，你这样会少算12个小时我觉得。

南气斑织逍遥

江奈何 发表于 2019-9-20 14:19
我也是之前遇到过求位相角的问题，比你这个还复杂我还要求相差，只是我想不出更好的方法了，所以只能分类 ...

考虑到 arctan() 的计算值介于 -pi/2 至 pi/2 之间，在第四象限时，假如拿（2，-2）这个点来算，就是-pi/4，其真实值应记为 7pi/4，所以我加的是2pi，后来我尝试换了几个波，效果都不如1波的pattern，可能是文中的图还用了什么其他的方法吧