点击下方公众号，回复资料，收获惊喜
滤波的原理
滤波的含义和应用就不再赘述，先理清楚几个概念：频谱就是频率的分布曲线，复杂振荡分解为振幅不同和频率不同的谐振荡，这些谐振荡的幅值按频率排列的图形叫做频谱。设一个能量信号为f(t)，则它的频谱密度F(ω)可以由傅里叶变换求得。
即频谱（频谱密度函数/振幅密度函数）的图横坐标应该是频率，纵坐标是振幅。下述中谱即代表频谱（频谱密度函数），从数学角度上看，滤波过程实际上是原始序列经过一定的变换转化为另一序列的过程。时间序列使用的是数字滤波器，是一个线性运算系统，从输入的时间序列（时间函数）x(t)，后到输出新的时间函数y(t)，所经过运算是
其中，c(k)是脉冲函数，也称为脉冲响应，亦称权函数。若称输入时间函数的谱为X(f)，输出时间函数的谱为Y(f)，则将y(t)的运算代入，恒等变形后得到脉冲函数的谱为频率响应函数Ｈ(f)
推导过程：
一般而言，H(f)是复数，实部和虚部分别是
由此可见，H(f)的模|H(f)|是频率为f 的成分在输出序列中的振幅较之输入序列中的振幅增长的倍数，称为振幅响应函数或增益函数,|H(f)|＝1的频率成分滤波前后振幅将不变。
H(f)频率响应函数（振幅响应函数/增益函数/响应函数）的横坐标应是频率frequency，纵坐标是响应response，其峰值应该小于1，显而易见响应值越大的频率区间（截断频率区间），即是滤波后被保留下来的所需的波的区间。例如以下两张图片中的函数
定义：H (f)的辐角argH (f) 是频率为f的成分在输出序列中的初位相相对于输人序列初位相的位相移动，称为位相响应函数或相移函数。在一般气象应用中，不希望滤波后产生位相移动，这要求arg=0,这往是不能完全做到的，实际应用中可先使虚部为0，并取偶函数c(k)=c(-k)以达到接近效果。
实际应用中时间t时离散的等间距的，输入序列的长度有限，则输出序列可表达为有限项求和的形式，Ck为权重函数（权系数/权重系数）：
所以确定了m和Ck即确定了一个滤波器，即m和Ck是表征滤波器特征的两个特征量。
同时由于观测到的大气或海洋变量序列包含的频率成分很丰富，但是并非都很重要，某些成分甚至不一定是实际信号。而滤波器只能提取某些频率成分，不能识别周期性强弱和滤出的成分是否有意义。滤波后的结果总是较光滑和规则的，画成图呈现很好的正弦或余弦波型，但是，这不能说明原序列中包含着这种成分。因此，为了得到有意义的结果必须在滤波前或滤波后做统计检验。这部分涉及功率谱就先暂时略过，以后有空再记录。
Lanczos滤波模型
Lanczos滤波器模型描述如下：
式中，x(t)是t时刻的时间序列输入；yL(t)是t时刻的滤波输出；w(t)是Lanczos滤波器的权重函数（权系数，相当于Ck），描述为
由前面公式可知，它由(2n 1)个数构成。故原序列经上式的Lanczos滤波,所得输出序列{yL(t), t =n 1,n 2,…,m -n}，较原序列缩短2×n×dt。
NCL中的filwgts_lanczos函数
对于滤波，首先需要构造滤波器，根据需要分为三种：高通滤波器，低通滤波器以及带通滤波器。这里面的“通”字指的是通过滤波器的信号，也就是我们要的信号，比如对于高通滤波，即指高频的信号部分通过了滤波器，而其它则被滤掉。NCL提供了filwgts_lanczos函数来构造Lanczos滤波器，其参数总共有nWgt, ihp, fca, fcb, sigma这五个，参数含义取值参见下方或官网函数说明。构造完滤波器，使用filwgts_lanczos函数可以被用于产生一组具有由用户指定的特征权重。运用wgt_runave或 wgt_runave_Wrap功能可用于应用的权重（即平滑滤波的权重）。需要注意的是，无法避免的在使用wgt_runave 或wgt_runave_Wrap在过滤序列的开始和结束时丢失数据。
filwgts_lanczos (nwt, ihp, fca, fcb, nsigma)
nwt：
表示权重系数的总数，为标量（必须为奇数；nwt > = 3）。权重越大，过滤器越好。必须使用不同数量的权重进行迭代才能获得所需的响应。通常需要注意以下事项：过滤器越窄，所需的权重就越大，因此，丢失的每一端数据也就越多。
关于为何是奇数，是由于之前的有限项求和公式中，共有2m 1项，两端各丢失m×dt，若是最常用的等权滑动平均，例如m=2，则为2m 1=5项的等权滑动平均滤波。
ihp
指示低通滤波器的标量：ihp = 0; 高通 ihp = 1; 带通ihp = 2。
fca
标量，指示理想的高通或低通滤波器的截止频率：（0.0 < fca <0.5）。
fcb
仅在需要带通滤波器时使用的标量。它是第二截止频率，否则设为-999（fca<fcb <0.5）。
nsigma
表示sigma系数的幂的标量（nsigma > = 0）。注意：通常设为nsigma = 1（这个参数还没有搞懂，感觉应该是和教材所说取权系数Ck为正态分布概率密度数值时的正太分布的均方差，sigma取多大应该与m有关）
返回值
返回数组长度为nwt的一维数组。如果fca为double，则类型为double；否则为float。
检查滤波器：
当设置好这一系列参数之后，对于返回响应函数（频率和幅度）作为返回权重（例如wt）的属性。具体来说：属性wt @ freq和wt @ resp 是长度为（2 * nwt 3）的一维数组，并且与wt的类型相同。通过这两个属性可以来绘制频率响应函数（横坐标wt@freq,纵坐标wt@resp），通过频率响应函数来检查滤波器是否满足需求
容易出错的地方
过滤器功能按时间序列中的离散时间步长运行。他们的时间指标是时间步长。
为了获得正确的滤波器参数，必须手动通过“时间步长”计算表示截止频率（fca,fcb）。
例如，如果所需的过滤器是10到50天，并且时间序列是按3天的时间步长，则：
dt =每个时间步长3天（若为每小时6小时一次资料序列，显然dt=1/6,小于1）
t1 = 50天（低频截止时间，以周期表示）
t2 = 10天（高频截止时间，以周期表示）
fca = dt / t1 = 3/50= 0.06个时间步长（低频截止频率）
fcb = dt / t2 = 3/10= 0.30时间步长（高频截止频率）
所以计算截止频率时应注意dt,t1,t2的时间单位的统一性。
