13.58 spectrum1d

官方文档:spectrum1d
简介:计算一个时间序列的自功率谱,或两个时间序列的互功率谱

spectrum1d 从标准输入流或数据文件中读取一列或两列数据。这些数据被当作是采样间隔为 <dt> 的等间隔的时间序列。spectrum1d 采用Welch方法,即加窗多段平均周期图法,计算输出自功率或互功率谱密度。其输出的功率谱的标准差,是用 Bendat和Piersol提供的算法。

spectrum1d 的输出文件有三列:

f|w   p   e

其中, f代表频率,w代表波长,p代表计算的功率谱密度,e代表一个标准差的值。

spectrum1d 的输出文件的文件名是使用统一的前缀 name_stem 。如果使用了 -C 选项,那么将会有8个文件输出,否则只生成一个功率谱文件( .xpower )。这些文件默认是以ASCII码格式,除非用 -bo 选项指定为二进制格式输出。这8个文件介绍如下:

  1. name_stem.xpower: X(t)的功率谱。单位是 X*X*dt
  2. name_stem.ypower: Y(t)的功率谱。单位是 Y*Y*dt
  3. name_stem.cpower: 一致性(coherent)的功率谱。单位和 ypower 一样。
  4. name_stem.npower: 噪声的功率谱。单位和 ypower 一样。
  5. name_stem.gain: 增益谱,或传输函数的模。单位是 Y/X
  6. name_stem.phase: 相位谱,或传输函数的相位。单位是弧度。
  7. name_stem.admit: 导纳(Admittance)谱,或传输函数的实部。单位是 Y/X
  8. name_stem.coh: (平方)相干谱,或者线性相关系数(它是频率的函数)。i 无单位,取值范围为 [0,1] 。信噪比 SNR=coh/(1-coh) 。当 coh=0.5 时,SNR=1。

除非使用 -T 选项,否则以上文件会以单个文件单列的形式输出。

13.58.1 选项

-S<segment_size>
<segment_size> 是一个2的指数数值,用于控制Welch方法中分段平均时的窗口长度。它也决定了功率谱密度的最小频率分辨率和最大频率分辨率,即 1.0/(segment_size*dt)1.0/(2*dt) (即Nyquist频率)。在功率谱密度中的一个标准误差大约为 1.0/(n_data/segment_size) ,比如 segment_size=256,那么就需要25600个数据点去计算一个误差棒的10%。互功率谱误差棒的计算则需要更多数据点,而且是相干性的函数,比较复杂。
table
输入文件名。它是ASCII类型的一列数据或两列数据。如果是一列数据文件,就计算自功率谱;如果是两列,就计算互功率谱。若未指定文件名, spectrum1d 会从标准输入流中读取数据。
-C[xycnpago]
默认会输出全部8个文件。使用该选项可以指定输出8个文件中的某些文件。x=xpower、y=ypower、c=cpower、n=npower、p=phase、a=admit、g=gain、o=coh。
-D<dt>
设置读入的时间序列的时间采样间隔,默认值是1。
-L[m|h]
不去除信号中的线性趋势。默认情况下,在对信号进行变换处理前会先去掉其中的线性趋势。 m 表示去掉数据的均值, h 表示去掉数据的中值。
-N[name_stem]
输出文件名的前缀,默认为 spectrum 。若不使用此选项,则输出的8个文件会合到一个文件里。
-T
不让单个分量的结果输出到标准输出流。
-W
输出文件中第一列是波长而不是频率。默认输出时第一列是频率。

13.58.2 示例

假设 data.g 是重力数据,单位为 mGal,空间采样间隔为1.5 km。如下命令会输出数据的功率谱,单位为 mGal^2 km 表示:

gmt spectrum1d data.g -S256 -D1.5 -Ndata

假设你除了有重力数据 data.g 之外,还有在相同地点测得的地形数据 data.t ,单位为 m。计算二者之间的传输函数,即 data.t 是输入, data.g 是输出:

paste data.t data.g | gmt spectrum1d -S256 -D1.5 -Ndata -C > results.txt