元々の時系列データ。サンプル数は偶数でなければならない。
The original time series data.
The number of samples must be an even number.
new_dt
変更後のサンプリング間隔。seq.dtの整数分の1でなければならない。
The new sampling interval,
which must be an integer divisor of seq.dt.
◆戻り値(Return value)
引数seqが表す時系列データのサンプリング間隔をnew_dtに変更した時系列データ。
seqのナイキスト周波数以上の高周波成分を含まないという条件のもとで
サンプリング間隔を変更する。
A time series data obtained by changing the sampling interval
of the original time series data (argument seq) to new_dt.
The sampling interval is changed keeping zero spectral amplitude
for the frequencies higher than or equal to
the Nyquist frequency of seq.
引数seqが表す元々の時系列データを\(f(t)\)、
サンプリング間隔変更後の時系列データ(戻り値)を\(g(t)\)とする。
\(f(t)\)のデータサンプル数を\(2J\)(偶数)、先頭時刻を\(t_0\)、
時間刻みを\(\Delta t_{old}\)とする。
すなわち\(f(t)\)は時刻\(t=t_0+j\Delta t_{old}\), \(j=0,1,2,\cdots,2J-1\)
において与えられているものとする。
また、\(g(t)\)は同じ先頭時刻\(t_0\)から始まるものとし、
データサンプル数を\(2K\)(偶数)、時間刻みを\(\Delta t_{new}\)とする。
すなわち\(g(t)\)は時刻\(t=t_0+k\Delta t_{new}\), \(k=0,1,2,\cdots,2K-1\)
において与えられているものとする。
以下では\(f_j\equiv f(t_0+j\Delta t_{old})\),
\(g_k\equiv g(t_0+k\Delta t_{new})\)という表記を用いる。
\(\Delta t_{old}/\Delta t_{new}\)は整数になるものとし、
その整数を\(L\)とおく。
このとき\(K=JL\)とすれば\(f(t)\)と\(g(t)\)の長さは等しくなる。
Let \(f(t)\) be the original time series data
represented by argument seq,
and \(g(t)\) be the time series data after changing the sampling interval
(i.e., the return value).
Let the number of data samples, beginning time,
and sampling interval of \(f(t)\) be
\(2J\) (an even number), \(t_0\), and \(\Delta t_{old}\), respectively.
This means that \(f(t)\) is given at time
\(t=t_0+j\Delta t_{old}\), \(j=0,1,2,\cdots,2J-1\).
Let \(g(t)\) have the same beginning time \(t_0\),
and the number of data samples and sampling interval of \(g(t)\) be
\(2K\) (an even number) and \(\Delta t_{new}\), respectively.
This means that \(g(t)\) is given at time
\(t=t_0+k\Delta t_{new}\), \(k=0,1,2,\cdots,2K-1\).
Below, we use the notations \(f_j\equiv f(t_0+j\Delta t_{old})\)
and \(g_k\equiv g(t_0+k\Delta t_{new})\).
The ratio \(\Delta t_{old}/\Delta t_{new}\) is assumed to be an integer
and denoted by \(L\).
Then, using \(K=JL\) realizes the same data lengths for
\(f(t)\) and \(g(t)\).
\(f(t)\)と\(g(t)\)の離散フーリエ変換を考える。
\(f(t)\)と\(g(t)\)は長さが等しいのでフーリエ変換における周波数刻みは等しく
\[\begin{equation}
\Delta f=\frac{1}{2J\Delta t_{old}}=\frac{1}{2K\Delta t_{new}}
\label{eq.delta_f}
\end{equation}\]
である。この\(\Delta f\)を用いて\(f(t)\)の離散フーリエ変換は
\[\begin{eqnarray}
F_n
&\equiv& F(n\Delta f)
\nonumber \\
&=& \sum_{j=0}^{2J-1}f_j
\exp\left[2\pi in\Delta f(t_0+j\Delta t_{old})\right]\Delta t_{old}
\nonumber \\
&=& \exp\left(2\pi in\Delta f t_0\right)\Delta t_{old}
\sum_{j=0}^{2J-1}f_j\exp\left(2\pi inj\Delta f\Delta t_{old}\right)
\nonumber \\
&=& \exp\left(2\pi in\Delta f t_0\right)\Delta t_{old}
\sum_{j=0}^{2J-1}f_j\exp\left(\frac{2\pi inj}{2J}\right)
\nonumber \\
& & \left(n=0,1,2,\cdots,2J-1\right)
\label{eq.F_n}
\end{eqnarray}\]
と書ける。
このスペクトルにおいて\(n=J\)がナイキスト周波数に対応する。
ナイキスト周波数を挟んで高周波側と低周波側のスペクトルの間には
\[\begin{eqnarray}
F_{2J-n}
&=& \exp\left[2\pi i(2J-n)\Delta f t_0\right]\Delta t_{old}
\sum_{j=0}^{2J-1}f_j\exp\left[\frac{2\pi i(2J-n)j}{2J}\right]
\nonumber \\
&=& \exp\left(4\pi i J\Delta f t_0\right)
\exp\left(-2\pi in\Delta f t_0\right)
\Delta t_{old}
\nonumber \\
& & \sum_{j=0}^{2J-1}
f_j\exp\left(2\pi ij\right)\exp\left(-\frac{2\pi inj}{2J}\right)
\nonumber \\
&=& \exp\left(\frac{2\pi i t_0}{\Delta t_{old}}\right)
\exp\left(-2\pi in\Delta f t_0\right)\Delta t_{old}
\sum_{j=0}^{2J-1}\exp\left(-\frac{2\pi inj}{2J}\right)
\nonumber \\
&=& \left(F_n\right)^{*}
\exp\left(\frac{2\pi i t_0}{\Delta t_{old}}\right)
\label{eq.F_n.over_Nyquist}
\end{eqnarray}\]
の関係が成り立ち、特に\(t_0/\Delta t_{old}\)が整数ならば
\[\begin{equation}
F_{2J-n}=\left(F_n\right)^{*}
\label{eq.F_n.over_Nyquist.simple}
\end{equation}\]
というよく知られた関係が得られる。ここで\(^{*}\)は複素共役を表す。
なお、以下の取り扱いでは(\ref{eq.F_n.over_Nyquist})式を用いるので
\(t_0/\Delta t_{old}\)が整数であることを必要としない。
Now, let us consider the discrete Fourier transformations
of \(f(t)\) and \(g(t)\).
As the lengths of \(f(t)\) and \(g(t)\) are equal,
the frequency intervals of the Fourier transformations of them
are also equal, and are given by eq. (\ref{eq.delta_f}).
Using this \(\Delta f\), the discrete Fourier transformation of \(f(t)\)
can be written as (\ref{eq.F_n}).
In this spectrum, \(n=J\) corresponds to the Nyquist frequency.
Over the Nyquist frequency, higher and lower spectral contents
are related by eq. (\ref{eq.F_n.over_Nyquist}),
which reduces to the well-known relation
of eq. (\ref{eq.F_n.over_Nyquist.simple}) if \(t_0/\Delta t_{old}\),
where \(^{*}\) denotes the complex conjugate.
In the treatment below, we use eq. (\ref{eq.F_n.over_Nyquist}),
meaning that \(t_0/\Delta t_{old}\) is not needed to be an integer.
\(g(t)\)の離散フーリエ変換として
\[\begin{equation}
G_n
\equiv G(n\Delta f)
=\begin{cases}
F_n & (n=0,1,2,\cdots,J-1) \\
0 & (n=J,J+1,J+2,\cdots,K)
\end{cases}
\label{eq.G_n}
\end{equation}\]
という形のスペクトルを考える。
すなわち、\(f(t)\)のナイキスト周波数未満の周波数成分については
\(f(t)\)と\(g(t)\)のスペクトルは全く同じで、
\(f(t)\)のナイキスト周波数以上の周波数成分については
\(g(t)\)のスペクトルは0になるものとする。
なお\(n=K\)が\(G_n\)のナイキスト周波数に対応し、
これよりも高周波側の成分については(\ref{eq.F_n.over_Nyquist})式にならって
\[\begin{equation}
G_{2K-n}=\left(G_n\right)^{*}\exp\left(\frac{2\pi it_0}{\Delta t_{new}}\right)
\label{eq.G_n.over_Nyquist}
\end{equation}\]
で与えるものとする。このように定義した\(G_n\)を用いて
離散フーリエ逆変換は
\[\begin{eqnarray}
g_k
&=& \sum_{n=0}^{2K-1}
G_n \exp\left[-2\pi in\Delta f(t_0+k\Delta t_{new})\right] \Delta f
\nonumber \\
&=& \Delta f \sum_{n=0}^{2K-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-2\pi ink\Delta f\Delta t_{new}\right)
\nonumber \\
&=& \Delta f \sum_{n=0}^{2K-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\nonumber \\
&=& \Delta f \sum_{n=0}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\nonumber \\
& & +\Delta f \sum_{n=1}^{J-1}
G_{2K-n} \exp\left[-2\pi i(2K-n)\Delta f t_0\right]
\exp\left[-\frac{2\pi i(2K-n)k}{2K}\right]
\nonumber \\
&=& \Delta f \sum_{n=0}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\nonumber \\
& & +\Delta f \sum_{n=1}^{J-1}
\left(G_n\right)^{*}\exp\left(\frac{2\pi it_0}{\Delta t_{new}}\right)
\exp\left(-4\pi iK\Delta f t_0\right)
\exp\left(2\pi in\Delta f t_0\right)
\nonumber \\
& & \exp\left(-2\pi ik\right)
\exp\left(\frac{2\pi ink}{2K}\right)
\nonumber \\
&=& \Delta f \sum_{n=0}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\nonumber \\
& & +\Delta f \sum_{n=1}^{J-1}
\left(G_n\right)^{*}\exp\left(\frac{2\pi it_0}{\Delta t_{new}}\right)
\exp\left(\frac{-2\pi i t_0}{\Delta t_{new}}\right)
\exp\left(2\pi in\Delta f t_0\right)
\nonumber \\
& & \exp\left(\frac{2\pi ink}{2K}\right)
\nonumber \\
&=& \Delta f \sum_{n=0}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\nonumber \\
& & +\Delta f \sum_{n=1}^{J-1}
\left(G_n\right)^{*}
\exp\left(2\pi in\Delta f t_0\right)
\exp\left(\frac{2\pi ink}{2K}\right)
\nonumber \\
&=& \Delta f \sum_{n=0}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\nonumber \\
& & +\left[
\Delta f \sum_{n=1}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\right]^{*}
\nonumber \\
&=& \Delta f G_0
+\Delta f \sum_{n=1}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\nonumber \\
& & +\left[
\Delta f \sum_{n=1}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\right]^{*}
\nonumber \\
&=& \Delta f G_0
+2Re\left[
\Delta f \sum_{n=1}^{J-1}
G_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\right]
\nonumber \\
&=& \Delta f F_0
+2Re\left[
\Delta f \sum_{n=1}^{J-1}
F_n \exp\left(-2\pi in\Delta f t_0\right)
\exp\left(-\frac{2\pi ink}{2K}\right)
\right]
\label{eq.g_k}
\end{eqnarray}\]
と書ける。
As the discrete Fourier transformation of \(g(t)\),
we consider a spectrum of the form given by eq. (\ref{eq.G_n}).
Here, the spectra of \(f(t)\) and \(g(t)\) are perfectly equal
at frequencies less than the Nyquist frequency of \(f(t)\),
and at frequencies greater than or equal to the Nyquist frequency of \(f(t)\),
the spectrum of \(g(t)\) is zero.
The Nyquist frequency for \(G_n\) is represented by \(n=K\),
and the spectrum above this frequency is given by
eq. (\ref{eq.G_n.over_Nyquist}),
which was derived in the same manner as (\ref{eq.F_n.over_Nyquist}).
Using this \(G_n\),
the discrete Fourier inverse transformation is given by
eq. (\ref{eq.g_k}).
この関数では(\ref{eq.F_n})式を用いて
\(n=0,1,2,\cdots,J\)に対する\(F_n\)を計算し、
次にその結果を用いて(\ref{eq.g_k})式により\(g_k\)を計算する。
高速フーリエ変換のアルゴリズムは使用しない。
というのも、高速フーリエ変換を用いるにはサンプル数を2の巾乗にする必要があり、
その場合に2つの時系列データの周波数刻みを等しくする(\ref{eq.delta_f}式)ためには
\(\Delta t_{old}/\Delta t_{new}\)も2の巾乗でなければならず、
元の時系列データの1/5倍や1/10倍などの時間刻みに直した時系列データを
作成できなくなるからである。
In this function, \(F_n\) for \(n=0,1,2,\cdots,J\) are first calculated
using eq. (\ref{eq.F_n}),
and then \(g_k\) are calculated using eq. (\ref{eq.g_k}).
The algorithm of fast Fourier transformation is not used.
This is because the use of fast Fourier transformation requires
the number of samples to be a power of 2;
in this case, to realize the same frequency intervals
for the two time series data (eq. \ref{eq.delta_f}),
\(\Delta t_{old}/\Delta t_{new}\) is required to be a power of 2,
which makes it impossible to create a time series data
which has the sampling interval 1/5 or 1/10 of that of the original one.
◆検証(Validation)
この関数によって意図通りの波形を作成できるかを実際のデータを用いて検証した。
使用したのは2021/5/14 11:35:40からの80秒間の
御嶽山における名大観測点NU.FMI.Uの波形である。
もともと100 Hzでサンプリングされている。
これを読み込んでこの関数により1000 Hzサンプリングに変更した。
Whether this function creates an intended waveform
was checked by a real data.
Here, a 80-s-long waveform starting from 11:35:40 JST on 14 May 2021
at a station NU.FMI.U on Mt. Ontake operated by Nagoya University was used.
This data was originally sampled at 100 Hz,
and was changed to 1000 Hz sampling by this function.
下図は元の波形(上段)とこの関数によって1000 Hzサンプリングに変更した波形(下段)の
比較である。
両者は非常によく似た波形になっていることが分かる。
The figure below illustrates the original waveform (upper panel)
and the waveform resampled at 1000 Hz by this function (lower panel),
showing quite similar waveforms.
10秒間を引き伸ばして重ねた波形を下に示す。
赤が元の波形、緑が1000 Hzサンプリングの波形である。
よく合っていることが分かる。
The figure below shows 10-s-long waveforms overlain,
where red is the original waveform
and blue is the resampled (1000 Hz) waveform,
showing excellent fits.
時間軸を更に引き伸ばすと、
元のデータ(赤)は100 Hzサンプリングなのでそれらの点が折れ線で繋がっているが、
1000 Hzサンプリングの波形(緑)ではそれらが滑らかに繋がっている。
Further expanding the time exis show that
the original data (red) is a bended line due to 100 Hz sampling,
whereas the resampled (1000 Hz) data (green)
connects the points by smooth lines.
スペクトルの比較図を下に示す。
SACを用いてスペクトルを計算すると
末尾にダミーの0.0を埋めてサンプル数を2の巾乗にする関係で
少しだけずれが生じてしまう。
そこで、関数sft2_sequenceを用いてスペクトルを計算した。
全周波数範囲(上段・中断)で見ると、もともと50 Hzであったナイキスト周波数が
リサンプルによって500 Hzに変更されたが
この間の周波数成分は0であることが分かる。
また、ある周波数範囲を拡大して重ねてプロットしてみると(下段)
完全に一致していることが分かる。
The figure below shows a comparison of the spectra.
Note that calculating the spectra by SAC results in a slight difference
due to adding dammy zeroes at the end to make the number of samples
a power of 2.
To avoid it, the spectra were computed using function sft2_sequence.
The plot for the entire frequency band (upper and middle panels)
shows that the Nyquist frequency was changed from 50 Hz to 500 Hz,
keeping the spectral amplitude zero between these frequencies.
The lower panel shows an extension of a certain frequency band,
showing complete fit.
