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