16bit 44.1kHz の wav ファイルを扱うことになりました。
スケール変換しようと思いハードウェアの資料を探しましたが見当たりません。レコーダーですので、ないのかもしれません。
2^15 で割って正規化し FFTをかけて、振幅を抽出します。で、この後が不明。ある基準値をもって dB値に変換するのですが、音圧と異なるので決まりがないようです。最大値を基準にしても良いし、最小値でも問題ないのでしょう。今回は見た目で真ん中あたりを選択しましたが、決まったやり方があれば知りたいところです。
出来上がったランニングスペクトルを見ると、現象によって奇麗に判別できました。これで全ての音声を聞かなくて済みそうです。時短に有効活用です。
位相の unwrapping について調べていました。
SAR ではなく、単なる信号の再現（逆フーリエ）でこの話題が出てきました。
調べてみると、ありますね。
位相アンラッピング(Phase unwrapping)
http://retrofocus28.blogspot.com/2013/12/phase-unwrapping_26.html
def mywrap(y): for i in range(1, len(y)): div = y[i] - y[i-1] if div < -pi: y[i:] += 2*pi elif div > pi: y[i:] -= 2*pi else: pass return y
numpy なら1行です。
np.unwrap()
簡単ですね。
理想的な条件なら再現できるようですが、実際のデータでは難しいようです。
Timothy Clements, Marine A. Denolle
SeisNoise.jl: Ambient Seismic Noise Cross Correlation on the CPU and GPU in Julia
Seismological Research Letters (2021) 92 (1): 517–527.
https://doi.org/10.1785/0220200192
https://docs.juliahub.com/SeisNoise
https://github.com/tclements/SeisNoise.jl
ambient seismic noise cross correlation のライブラリです。これまで見てきた passive 手法の中では最も新しく、シンプルに書けそうです。しかも CUDA 対応済み。
https://phreeqc.blogspot.com/2021/05/passive.html
https://phreeqc.blogspot.com/2021/06/blog-post_25.html
ちょうど手元の Win10 ＋ Jupyter がエラーを吐いて動かない状態でしたので、再構築がてらインストールしてみました。
まずは、既存の Julia と miniconda をアンインストール。関連フォルダをすべて削除。レジストリもクリーンにしてから再起動。
Julia 1.65 をインストール後、
pkg> add IJulia
pkg> add SeisNoise
pkg> add SeisIO; build; precompile
pkg> add Plots
miniconda 3 をインストール後、
$ conda install jupyter lab
$ jupyter lab
規定の browser だと表示されなかったので、起動 browser を変更。
$ jupyter lab --generate-config
~/.jupyter/jupyter_lab_config.py
c.ServerApp.browser = '"C:\\Program Files (x86)\\Microsoft\\Edge\\Application\\msedge.exe\" %s'
これで立ち上がりました。
Github の README.md のソースをコピペして run。
同じ絵が出ました。ＯＫです。
F-net の sac が手元になかったため、FDSN(IRIS) にある気象庁の静岡と岐阜のデータに変更。周波数を低めに設定し、run。
１日間と短いのですが、そこそこの形になりました。
DL に時間がかかりますが、仕方がないでしょう。
テストは CPU ＋1日間のデータのみでしたが、感覚的に速い方だと思います。GPU を使えるなら、大きなデータを扱えそうです（データを用意するのに時間がかかりますが）。
exsample は今後整備される予定のようです。楽しみですね。
連休中に、スペクトルの平滑化プログラムを書いておこうと思いつきました。
大崎本にサブルーチンが載っていますので、それを利用すれば完成するでしょうと高を括っていました。
が、本を持ち帰っていませんでした。
どうしましょうか、と検索してみると、大崎総合研究所から6月にソースコードが公開されていました。（正誤表もあります。）
http://www.ohsaki.co.jp/activity/download/index.html
昔の有名な図書ですので、プロは既にアレンジしたコードをお持ちです。が、私のようなアマにはありがたい。感謝です。
サブルーチンをダウンロードして、プログラムに組み込み、動かしてみました。が、正しく動きません。どこだ？と探しているうちに、連休が終わりました。
そういえば以前、倍精度に変更したサブルーチンをプロからいただいていました（今回も、連休中に頂きました）。それに組み替えると、改善。さらに、引数を変数に変えると正しく動きました。後者はともかく、前者はあるあるでしたね。感謝。
昨日、お礼がてら雑談していると、周波数によってバンド幅を変える方法もあると教えていただきました。バンド幅を決めて平滑化すると、高周波側と低周波側で見え方が変わってきます。それが改善されるのでしょう。良いですね。
といっても、それを使うのに妥当なスペクトルか判断できない知識レベルですので、今回は通常の平滑化で完成。バグ取りに数日かかりましたが、喰わせるデータや単位のチェックもあわせてできたので良しとしましょう。
文献を読んでいると、「Hilbert-Huang Transform（ヒルベルトファン変換）」が出てきました。
初めて知りました。何だろうと調べてみると、少ないながら日本語の解説がありました。
大塚ほか「Hilbert-Huang Transformによる非定常・非線形信号解析」寒地土木研究所月報, 2011
FFT では時間分解能を上げると周波数分解能が下がり、周波数分解能を上げると時間分解能が落ちます。両方上げたい場合、衝撃など短時間の波形解析の場合に利点でしょうか？
どこかの分野ではよく使われているのでしょうか？
Python にも HHT ライブラリはあるので対応可能でしょう。が、少し使いづらい？異なるデータセットに対して同じ周波数の値を得られた方が、あとで使いやすいかな。
ま、頭の片隅に残しておきましょう。
ω=2πf=2π/T・・・角振動数ωは振動数ｆの２π倍
\begin{align*}cos{\omega t},\ \sin{\omega t}・・・ω=1(=cost,\ \sin{t})\end{align*}で周期T=２πの周期運動
\begin{align*}x=x_0cos{\left(\omega t+\alpha\right)}\end{align*}調和運動：()の中身が ｔ の一次式
自由振動の一般解
\begin{align*}m\ddot{x}+kx=0\end{align*}
【表現1】解の形を\begin{align*}x=Xe^{\lambda t}, λ＝iω\end{align*}とおいて展開
\begin{align*}x=X_1e^{i\omega t}+X_2e^{-i\omega t}\end{align*}２つの任意定数を含んだ解＝一般解
【表現2】オイラーの公式を用いて書き換え
\begin{align*}x=\left(X_1+X_2\right)cos{\omega t}+i\left(X_1-X_2\right)sin{\omega t}\end{align*}
\begin{align*}=acos{\omega t}+bsin{\omega t}\end{align*}
【表現3】三角関数の合成の公式を用いて書き換え
\begin{align*}x=a_0cos{\left(\omega t+\alpha\right)}\end{align*}
周期関数のフーリエ級数
【表現2】の利用
\begin{align*}f\left(t\right)=\frac{1}{2}a_0+a_1cos{\omega t}+a_2cos{2\omega t+\ldots}\end{align*}\begin{align*}+b1sinωt+b2sin2ωt+⋯\end{align*}
オイラーの公式を利用して表現を変えると
\begin{align*}cos{\theta}=\frac{e^{i\theta}+e^{-i\theta}}{2},\ sin{\theta}=-i\frac{e^{i\theta}-e^{-i\theta}}{2}\end{align*}
\begin{align*}f\left(t\right)=\frac{1}{2}a_0+\sum_{r=1}^{\infty}\left(a_r\frac{e^{ir\omega t}+e^{-ir\omega t}}{2}{-jb}_r\frac{e^{ir\omega t}-e^{-ir\omega t}}{2}\right)\end{align*}
複素フーリエ級数
\begin{align*}f\left(t\right)=\sum_{r=‐∞}^{∞}{F_re^{ir\omega t}}\end{align*}
複素フーリエ係数
\begin{align*}F_r=\frac{1}{T}\int_{-T/2}^{T/2}{f\left(t\right)e^{-ir\omega t}dt}\end{align*}
非周期関数のフーリエ積分
\begin{align*}f\left(t\right)=\int_{-\infty}^{\infty}{F\left(\omega\right)e^{i\omega t}d\omega}\end{align*}
フーリエ変換・・・フーリエ級数のフーリエ係数に相当
\begin{align*}F\left(\omega\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{f\left(t\right)e^{-i\omega t}dt}\end{align*}
デルタ関数のフーリエ変換
t＝0の時、e^0=1より
\begin{align*}F\left(\omega\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{\delta\left(t\right)e^{-i\omega t}dt}=\frac{1}{2\pi}\end{align*}
微動計を用いた衝撃の測定に呼ばれました。
今まで他の方が頑張っていました。が、得られたデータがイマイチだったということが判明。で、地盤を対象に微動計を使っていた私も御一緒することに。
どうも、微動計に備わっていたフィルタの周波数特性をうまく活用できていなかったのが一因のようでした。
その選定から見直そうとしたのですが、私も応用が利くほど詳しくありません。調べていくうちにわからない用語が出てきました。標準モード、群遅延特性、デシメーション。
使うツールの詳細を理解していないと、正しい調査ができないのは共通。いままでプロの指導でうまくいっていただけのことであり、場を変えると途端にダメになるのは同じ。
ちょうど週末で時間があり、この機会に理解することにしました。
用語自体は信号処理分野で広く使われているようです。いくつかの図書を見ましたが、基礎中の基礎といった感じでした。
わかりやすかったのが図書よりも以下のサイト。確かに、音響ととらえると理解しやすい（ホールの設計者はこのようなことまで考えていらっしゃるのですね！）。
以下、解釈した内容を( ..)φメモメモ。
**************************************
ヤマハサウンドシステム
短期集中連載 超解説FIR！第1回～第3回
https://www.yamaha-ss.co.jp/published-articles/journals-04.html
1．用語
1－1．周波数特性・位相特性・群遅延特性
入力したインパルス信号の応答を周波数領域で見る（フーリエ変換）と各周波数のエネルギー応答がわかる。これが周波数特性。
各周波数が時間的にどのくらい遅れて到達しているかを見たものが位相特性。
位相を周波数で微分したものが群遅延特性。
周波数とその位相の変化量（遅れ）が比例していれば群遅延は一定（直線位相特性）。
ex.
1Hz→1波長（2π、1秒）
2Hz→2波長（4π、1秒）
3Hz→3波長（6π、1秒）
1－2．インパルス応答、伝達関数
システム（ブラックボックス）にインパルスを入力した時に出力された信号がインパルス応答。
インパルスは全周波数で同じエネルギー量持っているため、インパルス応答を見ただけでブラックボックスの特性=伝達関数が分かる。
インパルス応答を周波数の観点で見ると周波数特性が分かり、周波数ごとの到達時間を調べれば位相特性が分かる。
インパルス応答、伝達関数、周波数特性、位相特性等は、データの見方を変えただけ。
2．特性を作る（リバーブエフェクトなど）
2－1．FIR フィルタ
FIR の Tap にある乗算器の係数はフィルターのインパルス応答そのもの。
係数操作でインパルス応答を人工的に加工することが可能。IIR のように最小位相として位相が変化するインパルス応答などを作ることも可能。
目的とする周波数カーブを描いてそれを逆フーリエ変換すれば、カーブのインパルス応答が得られる。
そのインパルス応答を乗算器の係数に入れることで目的の周波数カーブのフィルターを作ることができる。
FIR フィルタは入力信号（サンプルデータ）とインパルス応答波形の２つを畳み込んで新たなデータを作り出す。つまり「畳み込みをするための演算器」。
2－2．FFT と窓関数
入力信号も IR も共に FFT して周波数領域のデータに変換。
畳み込み積分は周波数領域に変換すると単純な掛け算として非常に高速の計算が可能。
その結果を逆フーリエ変換（IFFT)して元の時間領域のデータに戻す。
ただし、FFT を行なうため処理するデータを切り出す際にサイドローブ（雑音）が発生。
窓関数でサイドローブの軽減処理をしているが、軽減しようとすると反対に周波数分解能が悪くなるという二律背反の関係がある。
周波数分解能： 〇 △ ×
サイドローブ軽減： × △ 〇
窓関数： ハミング ハニング ブラックマン
**************************************
インパルスにδ関数が使われているのは当然なのに新鮮。つながりますね。
次は、これらの処理順序について。
周波数の低下とともにパワースペクトラムが増加 するような信号の中で、パワースペクトラムの振幅が周波数に対して反比例する信号を「１/ｆゆらぎ」と呼ぶ。これでしたか！と目が覚めた気分。
g(x) のフーリエ変換 G(f) の f
1/距離：空間周波数（spatial frequency）
1/時間：周波数（frequency）
空間周波数 fの代わりに k=2πf がしばしば用いられる。先日の重力探査の文献では「kは波数」等と明記されていました。一方、Fourier domain の wave number と表記された文献もあったように思います。もう一度読みましょうか。
k-空間（またはフーリエ空間（Fourier domain/space)）：k=2πf
kは電磁波等の波数(k=2π/λ)と混同しないように。
