写在前面
写材料写累了，看文献也看累了，那就来写代码。
受困于之前使用NCL函数的烦恼，以及为了加深对于smthClmDayTLL函数的理解，这里打算使用python来尝试复现一下该函数的作用
数据
使用的数据是OLR-daily日平均数据进行测试：
https://psl.noaa.gov/data/gridded/data.olrcdr.interp.html
复现函数:smthClmDayTLL
以下是该函数的相关信息：
Prototype
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; This library is automatically loaded
; from NCL V6.2.0 onward.
; No need for user to explicitly load.
function smthClmDayTLL (
clmDay [*][*][*] : float or double,
nHarm [1] : integer
)
return_val [366][*][*] : typeof(clmDay)
传入两个参数：
1、日气候态平均数据，一般由ClmDayTLL 输出的三维数组(ntim、 lat、 lon)。通常，ntim 的大小为360(360 _ day)、365(365 _ day)或366(gregorian) ，这取决于与变量关联的日历属性。 2、用于构造平滑平均年周期的谐波数。通常，nHarm 选择1。当选择2时，意味着只使用年度和半年度谐波。
整个函数的流程是：
使用 clmDayTLL 函数的输出结果计算平滑的年循环数据。主要先通过 FFT 来完成的，将所有大于 nHarm 的谐波设置为零，并执行逆变换。
啥意思呢，就是说对于计算出来的日气候态平均结果平滑了一下。
NCL源码解析
那么，具体整个函数是如何来实现上面的过程呢。这需要打开他的源码：
undef("smthClmDayTLL")
function smthClmDayTLL (clmDay[*][*][*]:numeric, nHarm:integer)
;
local nFill, dn, z, cf, clmDaySmth
begin
dn = getvardims(clmDay) ; get dimension names
if (dn(0).ne."year_day") then
print("smthClmDayTLL: Warning: Usually expect year_day to be the dimension name")
end if
z = clmDay($dn(1)$|:,$dn(2)$|:,$dn(0)$|:); reorder make time fastest varying dimension
cf = ezfftf( z ) ; [2] x [nlat] x [mlon] x [nday/2+1]
; remember NCL is 0-based
; cf(:,0:nHarm-1) are retained unaltered
cf(:,:,:,nHarm ) = 0.5*cf(:,:,:,nHarm) ; mini-taper
cf(:,:,:,nHarm+1:) = 0.0 ; set all higher coef to 0.0
z = ezfftb( cf, cf@xbar) ; reconstructed series
clmDaySmth = z($dn(0)$|:,$dn(1)$|:,$dn(2)$|:)
clmDaySmth@information = "Smoothed daily climatological averages"
clmDaySmth@smoothing = "FFT: "+nHarm+" harmonics were retained."
return(clmDaySmth)
end
下面来简单解释一下每一步的作用：
1、使用getvardims函数获取输入数据的维度，一般是year_day，lat，lon，然后进行一个简单的逻辑判断，如果第一个维度是"year_day"则进行下面的操作，否则打印警告信息
(0) year_day
(1) lat
(2) lon
2、转换输入的数据的维度，从
year_dayxlatxlon变为latxlonxyear_dayxlatxlon,这是为了方便后续进行傅里叶变换3、使用
ezfftf()函数对于转换后的数组进行傅里叶变换，返回一个拓展一维的数组，从三维变四维，即，变为：[2] x [nlat] x [mlon] x [nday/2+1],最左边新增的维度分布表示实系数和虚系数
其中，ezfftf()调用的是fortran函数，如下所示,可以看出写的版本还是较老的
SUBROUTINE EZFFTF (N,R,AZERO,A,B,WSAVE)
C
C VERSION 3 JUNE 1979
C
DIMENSION R(*),A(1),B(1),WSAVE(1)
IF (N-2) 101,102,103
101 AZERO = R(1)
RETURN
102 AZERO = .5*(R(1)+R(2))
A(1) = .5*(R(1)-R(2))
RETURN
103 DO 104 I=1,N
WSAVE(I) = R(I)
104 CONTINUE
CALL RFFTF (N,WSAVE,WSAVE(N+1))
CF = 2./FLOAT(N)
CFM = -CF
AZERO = .5*CF*WSAVE(1)
NS2 = (N+1)/2
NS2M = NS2-1
DO 105 I=1,NS2M
A(I) = CF*WSAVE(2*I)
B(I) = CFM*WSAVE(2*I+1)
105 CONTINUE
IF (MOD(N,2) .EQ. 1) RETURN
A(NS2) = .5*CF*WSAVE(N)
B(NS2) = 0.
RETURN
END
以下是在fortran调用的RFFTF的函数说明，也是fortran写的。这里就先不去解释了，因为在python里面我们可以更简单实现。
rfftf.f
SUBROUTINE RFFTF (N, R, WSAVE)
C***BEGIN PROLOGUE RFFTF
....
....
Warning: this routine is not intended to be user-callable.
....
....
C***SUBSIDIARY
C***PURPOSE Compute the forward transform of a real, periodic sequence.
C***LIBRARY SLATEC (FFTPACK)
C***CATEGORY J1A1
C***TYPE SINGLE PRECISION (RFFTF-S, CFFTF-C)
C***KEYWORDS FFTPACK, FOURIER TRANSFORM
C***AUTHOR Swarztrauber, P. N., (NCAR)
C***DESCRIPTION
C
C ********************************************************************
C * NOTICE NOTICE NOTICE NOTICE NOTICE NOTICE NOTICE *
C ********************************************************************
C * *
C * This routine uses non-standard Fortran 77 constructs and will *
C * be removed from the library at a future date. You are *
C * requested to use RFFTF1. *
C * *
C ********************************************************************
C
C Subroutine RFFTF computes the Fourier coefficients of a real
C periodic sequence (Fourier analysis). The transform is defined
C below at output parameter R.
C
C Input Arguments
C
C N the length of the array R to be transformed. The method
C is most efficient when N is a product of small primes.
C N may change so long as different work arrays are provided.
C
C R a real array of length N which contains the sequence
C to be transformed.
C
C WSAVE a work array which must be dimensioned at least 2*N+15
C in the program that calls RFFTF. The WSAVE array must be
C initialized by calling subroutine RFFTI, and a different
C WSAVE array must be used for each different value of N.
C This initialization does not have to be repeated so long as
C remains unchanged. Thus subsequent transforms can be
C obtained faster than the first. Moreover, the same WSAVE
C array can be used by RFFTF and RFFTB as long as N remains
C unchanged.
C
C Output Argument
C
C R R(1) = the sum from I=1 to I=N of R(I)
C
C If N is even set L = N/2; if N is odd set L = (N+1)/2
C
C then for K = 2,...,L
C
C R(2*K-2) = the sum from I = 1 to I = N of
C
C R(I)*COS((K-1)*(I-1)*2*PI/N)
C
C R(2*K-1) = the sum from I = 1 to I = N of
C
C -R(I)*SIN((K-1)*(I-1)*2*PI/N)
C
C If N is even
C
C R(N) = the sum from I = 1 to I = N of
C
C (-1)**(I-1)*R(I)
C
C Note: This transform is unnormalized since a call of RFFTF
C followed by a call of RFFTB will multiply the input
C sequence by N.
C
C WSAVE contains results which must not be destroyed between
C calls of RFFTF or RFFTB.
C
C***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
C Computations (G. Rodrigue, ed.), Academic Press,
C 1982, pp. 51-83.
C***ROUTINES CALLED RFFTF1
C***REVISION HISTORY (YYMMDD)
C 790601 DATE WRITTEN
C 830401 Modified to use SLATEC library source file format.
C 860115 Modified by Ron Boisvert to adhere to Fortran 77 by
C changing dummy array size declarations (1) to (*).
C 861211 REVISION DATE from Version 3.2
C 881128 Modified by Dick Valent to meet prologue standards.
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900131 Routine changed from user-callable to subsidiary
C because of non-standard Fortran 77 arguments in the
C call to CFFTB1. (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE RFFTF
3、然后这里需要用到输入的第二个参数，即谐波数，当nHarm=1时，需要进行两步处理：对于一波的所有数值乘上一个缩放系数0.5，对于高于1波的其他波设置为0 4、将上述的结果进行傅里叶逆变换
5、将逆变换结果的维度重新转换回去： [year_day,lat,lon]
以上就完成了整个平滑过程，返回的结果即:保留了年循环的一波的平滑数据,或者说，保留了1个谐波。
在后面继续计算异常时，anomaly = calcDayAnomTLL(olr, yyyyddd, climatology),会减去这个平滑的气候态日平均数据，即：去除了第一谐波
或许有人要问，如果不平滑呢，结果会有什么差异吗。在后续分析以及查询官网资料过程中，对于后续异常的计算有些许影响，但是在后续进行滤波过程中可能影响不大。以下是官网的一个对比结果：
这里输出数据，准备与python的结果进行对比：
outfile=addfile("./smth2.nc","c")
outfile->smth=climatology
python复现
上面经过解析ncl的代码，对于整个计算流程其实已经很清晰了，下面在python中尝试复现一下。主要包含以下几个过程：
1、计算日气候态年平均结果 2、进行傅里叶变换 3、保留第一谐波 4、傅里叶逆变换
这里从他的频率上也可以看出来，第一个谐波表示一年的周期，第二个谐波表示半年的周期，第三个表示1/3 年的周期，第四个表示1/4 年的周期，依次类推...
下面是python和ncl的对比结果
绘图的标签不一致是因为一个基于xarray.plot()绘制，自动带有了经纬度范围，另一个是直接使用plt.contourf()绘制，只有数组的大小。实际上两个都是144x21的网格点数。
不得不说，在数值上还是有点差异，但是应该还是能接受的。而且发现貌似去除前3个谐波和2、1个都是类似的结果，挺有意思。
总结
尝试在python里面复现了ncl的平滑函数smthClmDayTLL,对于谐波的概念进一步有了理解，虽然复现结果仍然存在差异。对于后续的滤波需要进一步的检验。相关ncl的测试数据放到了GitHub上，olr的数据过大，这里不在GitHub上分享了，有兴趣的可以可以自己去下载一下或者使用其他的数据进行检验。
https://github.com/Blissful-Jasper/jianpu_record
★https://www.ncl.ucar.edu/Document/Functions/Contributed/smthClmDayTLL.shtml
https://github.com/yyr/ncl/blob/master/ni/src/examples/gsun/contributed.ncl
https://github.com/mansourmoufid/udsp/blob/master/fftpack/ezfftf.f
http://www.cs.yorku.ca/~roumani/fortran/slatecAPI/rfftf.f.html
https://www.ncl.ucar.edu/Document/Functions/Contributed/calcDayAnomTLL.shtml
给我点在看的人
越来越好看
