filt1文档

的filt1应用零相位巴特沃斯滤波器。需要Matlab的信号处理工具箱。

回到气候数据工具内容．

内容

语法

yf = filt1 (filtertype y“俱乐部”,fc) yf = filt1 (filtertype, y, Tc, Tc) yf = filt1 (filtertype y lambdac, lambdac) yf = filt1(…,fs, fs) yf = filt1(…,“x”,x) yf = filt1(…,Ts, Ts) yf = filt1…,“秩序”,FilterOrder) yf = filt1(…,“暗”,暗)[yf, filtb filta] = filt1(…)

描述

yf = filt1(filtertype,y，'fc'， fc)滤波1D信号y使用指定的filtertype和截止频率足球俱乐部．用于高通或低通滤波器足球俱乐部必须是标量。用于带通和带阻滤波器足球俱乐部必须是双元素数组。的filtertype可以

yf = filt1(filtertype,y，'Tc'，Tc)指定截止周期，而不是截止频率。该语法假设T = 1/f，与'lambdac'选项完全相同，但对于处理时间序列可能更直观一些。

Yf = filt1(filtertype,y，'lambdac'，lambdac)指定截止波长(秒)而不是截止频率。这个语法假设lambda = 1/f。

yf = filt1(…，'fs'， fs)指定采样频率Fs．如果既不“fs”，“x”,也不“t”是指定的,Fs = 1假定。

Yf = filt1(…，'x'，x)指定单调递增、等间隔采样次数或的向量x对应的位置y，用于确定采样频率。如果既不“fs”，“x”,也不“t”是指定的,Fs = 1假定。

yf = filt1(…，'Ts'，Ts)指定采样周期或采样距离Fs = 1/Ts．如果既不“fs”，“x”,也不“t”是指定的,Fs = 1假定。

Yf = filt1(…，'dim'，dim)指定要进行操作的维度。默认情况下，filt1对1D或2D数组沿第一个非单维操作，但对3D网格化数据集沿第3维操作。

yf = filt1(…，'order'，FilterOrder)指定顺序(有时称为滚边)的巴特沃斯滤波器。如果未指定的,FilterOrder = 1假定。

[yf,filtb,filta] = filt1(…)还返回过滤器分子filta和分母filtb．

例1:火车汽笛

对于这个例子，我们将使用内置的火车哨子示例文件，我们将添加一些高斯随机噪声，使事情变得有趣。

负载火车Y = Y +0.1*randn(size(Y));

原始信号可以像这样用黑色表示。首先，我们必须从采样频率建立一个时间向量Fs：

t =(0:长度(y)-1)/Fs;情节(t y“k -”，“线宽”, 1)箱从；轴紧包含的时间(秒)持有在

如果你有扬声器，你可以这样听火车汽笛声:

soundsc (y, Fs)

对火车汽笛进行高通滤波，只保留750赫兹以上的频率，并将其绘制在原始黑色信号的顶部:

Yhp = filt1(“惠普”, y,“fs”Fs,“俱乐部”, 750);持有在情节(t, yhp“r”）

或者，您可能想对原始信号进行低通滤波，以消除1100 Hz以下的频率。注:上面我们通过设置一个采样频率“fs”价值。你也可以定义一个向量作为自变量“x”．在这种情况下，自变量是时间，但对于空间滤波，它可能是沿着某些路径的累积距离。

火车汽笛的三个主频率在频率空间中间隔得很近，所以我们上面使用的默认一阶巴特沃斯滤波器不能消除750 Hz以下的所有能量。您可能希望通过指定使用更陡的滚落“秩序”,5．我们将用蓝色绘制低通过滤火车汽笛。

Ylp = filt1(“资讯”, y,“x”t“俱乐部”, 1100,“秩序”5);情节(t, ylp“b”）

使用plotpsd比较:

图plotpsd (y, Fs,“k”，“线宽”, 2)在plotpsd (yhp Fs,“r”) plotpsd (ylp Fs,“b”)包含的频率(赫兹)轴([600 1300 0 0.02])图例(原始信号的，“高通800赫兹”，.．.“低通1100赫兹”，“位置”，“西北”传说)boxoff

清晰的变量关闭所有

例2:地形概况

假设您有一个地形剖面，每10米测量40公里的高程。假设剖面有三个主要波长——761米、4公里和9.4米。概要文件可能是这样的。在例1中，我使用plotpsd绘制周期图

SpatialRes = 10;每10米采样一次x = 0:SpatialRes:40e3;%域0到40公里Lambda1 = 761;% 761米Lambda2 = 4000;% 4公里Lambda3 = 3000*pi;% ~9.4公里%生成配置文件:Y = rand(size(x)) + 5*sin(2* *x/lambda1) +.．.11*sin(2*pi*x/lambda2) + 15*sin(2*pi*x/lambda3);%地块概况:图(“位置”，[100 100 560 506])子图(211)图(x/1000,y，“k”，“线宽”, 2)在包含沿着某条路径的距离(km)ylabel的高度(米)盒子从轴紧%绘图功率谱:次要情节(212)plotpsd (y, x,“k”，“线宽”2,“数据库”，“日志”，“λ”)举行在包含的波长(公里)ylabel功率谱(dB)轴紧ylim(70年[-10])

上面你可以看到三个主要波长在周期图中的三个峰。

也许你想要消除我们添加到地形中的高频随机噪声兰德．要做到这一点，你可以低通过滤掉所有小于300米的波长:

Ylo = filt1(“资讯”, y,“x”, x,“lambdac”, 300);次要情节(211)情节(x / 1000, ylo,“r”) subplot(212) plotpsd(ylo,x，“r”，“数据库”，“日志”，“λ”）

上面，当我们低通滤波地形时，我们指定了一个数组x为路径距离，对应于y．或者，我们可以指定空间分辨率，也就是我们的采样距离“t”达到同样的效果。下面我们对原始地形进行高通滤波，去除波长超过6公里的图像。通过指定一个五阶巴特沃斯滤波器来使用一个紧密的滚转。

Yhi = filt1(“惠普”, y,“t”SpatialRes,“lambdac”, 6000,“秩序”5);次要情节(211)情节(x / 1000, yhi,“b”) subplot(212) plotpsd(yhi,x，“b”，“数据库”，“日志”，“λ”）

也许你想去除高频噪声和低频噪声。你可以通过对信号进行两次滤波，或者使用带通滤波器来实现。我们可以通过带通过滤掉所有短于3000米或长于3000米的波长，只保留功率谱中的中间峰值

Ybp = filt1(“英国石油公司”, y,“x”, x,“lambdac”(3000 5000),“秩序”3);次要情节(211)情节(x / 1000, ybp,“米”) subplot(212)“米”，“数据库”，“日志”，“λ”）

也许您只想删除一个频率范围。你可以通过从原始信号中减去一个带通信号来做到这一点:

Ybs = y - ybp;

或者你可以直接创建一个带阻滤波器，使用与我们使用带通滤波器相同的语法:

Ybs = filt1(“废话”, y,“x”, x,“lambdac”(3000 5000),“秩序”3);次要情节(211)情节(x / 1000, yb,“颜色”,(。0.45 . 82]) subplot(212) plotpsd(ybs,x，“颜色”,(。98 .45 .02]，“数据库”，“日志”，“λ”)传说(“原始”，低通滤波器的，“高反差保留”，.．.“带通”，“bandstop”，“位置”，“东北”传说)boxoff

例3:海冰范围

加载海冰范围的示例时间序列，只使用1989年至今的数据，因为旧数据的采样分辨率低于日分辨率。

负载seaice_extent.mat%获取1989年后的指数:Ind = t>datetime(1989,1,1);将数据集修剪为1989-2018:T = T (ind);extent_N = extent_N(ind);图(t,extent_N)轴紧盒子从ylabel北半球海冰范围(10^6 km^2)

我们能不能把季节周期过滤掉?当然deseason函数是一种方法，但如果我们用这个巴特沃斯过滤器呢?季节性周期的带阻滤波器要求我们定义拐角频率，那么6个月和2年之间的所有频率呢?

extent_N_filt = filt1(“废话”extent_N,“fs”, 365.25,“Tc”(0.5 - 2));持有在情节(t, extent_N_filt“线宽”, 2)

例4:网格化3D时间序列

假设您对墨西哥湾(GoM)的sst感兴趣:

负载pacific_sst图imagescn(lon,lat,mean(sst,3)) cmocean热% colormap持有在墨西哥湾的大致轮廓:Gomlon = [-91.4 -103.8 -98.8 -88.6 -82.8 -82.3];格姆拉特= [33.0 30.6 20.4 16.6 21.8 34.1];情节(gomlon gomlat,“罗”)文本(-92.4,23.9,“墨西哥湾”，“颜色”，“红色”，.．.“水平的”，“中心”，“绿色”，“机器人”）%设置文本对齐方式

我们可以用geomask而且当地的得到墨西哥湾sst的时间序列。首先，制作蒙版，并将其绘制为蓝色轮廓线，以确保蒙版在正确的位置:

[Lon,Lat] =网格(Lon,Lat);mask = geomask(Lat,Lon,gomlat,gomlon);轮廓(经度、纬度、双(面具),(0.5 - 0.5),“b”）

与面具定义时，获得墨西哥湾平均海温的一维时间序列很容易:

Sst_gom = local(sst,mask，“omitnan”）;图plot(t,sst_gom)轴紧datetick (“x”，“keeplimits”) ylabel“墨西哥湾海面温度\circC”

这个数据集是每月的时间分辨率，所以如果我们想对它进行低通滤波，只保留周期超过18个月的频率，我们可以在1d上操作sst1这样的数据集:

% lowpass filtered 1D sst1:Sst_gom_lp = filt1(“资讯”sst_gom,“Tc”, 18);持有在情节(t, sst_gom_lp)传说(“原始”，“低通滤波器过滤”）

你会注意到在1年的频率上仍然有大量的能量。那是因为巴特沃斯滤波器并不理想——在频率空间中，它的肩膀很宽。为了提高效率，在18个月以内更有效地削减能源，增加过滤器的顺序:

% lowpass filtered 1D sst1:Sst_gom_lp5 = filt1(“资讯”sst_gom,“Tc”, 18岁,“秩序”5);情节(t, sst_gom_lp5)传说(“原始”，.．.“低通滤波器过滤”，.．.“低通五阶”）

现在你可以看到一些年际变化不太受季节周期的影响。在任何应用程序中过滤掉大部分可变性可能不是一个好主意，但这里我们只是在玩弄函数的机制，所以我们只是说它很好，然后继续前进。

但是，如果不是只过滤一维数组，而是要过滤每个网格单元的时间序列呢?计算效率低的方法是嵌套几个循环，遍历网格的每一行和每一列，所有65x55网格单元格。但这意味着即使对于这个粗糙的网格，同样的操作也要做3000多次!幸运的是,filt1可以做得更有效率。如果自风场为三维矩阵，filt1已经知道如何在三维空间中操作了，尽管你总是可以指定“暗”,3如果你想更加确定的话。

下面是如何过滤整个三维网格SST时间序列使用filt1：

Sst_lp = filt1(“资讯”风场,“Tc”, 18);

现在我们可以看看墨西哥湾的时间序列过滤后的3D sst数据，应该完美匹配:

Sst_lp_gom = local(sst_lp,mask，“omitnan”）;情节(t, sst_lp_gom)传说(“原始”，.．.“低通滤波器过滤”，.．.“低通五阶”，.．.局部一阶后滤波器）

实际上，对过滤后的海表温度网格进行局部平均时间序列处理，得到的结果与对膳食局部时间序列进行滤波相同。这里有一个直接的比较，因为上面的图有点忙:

图绘制(sst_lp_gom sst_gom_lp,“。”)包含“先过滤，再局部平均”ylabel“先局部平均，再过滤”

两者之间的差异只是数值噪声。

作者信息

这个函数是气候数据工具箱Matlab．该功能和辅助文档由德克萨斯大学奥斯汀分校的Chad A. Greene编写。

发布与MATLAB®R2018b