はじめに

awkで標準出力から～(統計編)の続編になります。なので以下引用。

お仕事で、データに対してサクッと~~統計値求めたかったり、規格化したかったり、~~周波数解析したかったりする場合があります。
通常であればpythonとかRとかでサクッとやれちゃうんでしょうが、データ取得環境がオフラインでパッケージも入れられないような環境だったり、取得したデータが不等間隔のデータだったりしたので、awkでスクリプト作って標準出力のデータをパイプでつなげてサクッと求められるようにしました。

不等間隔データの周波数解析

Lomb-Scargle Periodogram

間隔が不均一なデータの場合、一般的なFFTの方法が採用できません。そこで、今回はLomb-Scargle法を採用しました。簡単に言うと任意周波数のsin関数でフィッティングすることで周波数成分求めようってやつですね。詳しくは下記とか見てください。

結構誰しもが考えそうな手法ですが、自分は$asin(wt-b)+c$の形のまま解こうとして挫折しました。頭悪かったですね。積和変換して$Asin(wt)+Bcos(wt)+C$で最小二乗すればよかったんですね。
(こんなこと言ってますけど実は真面目に式は導出してません…。)

ソース

とりあえず実装しました。

fft

awk '{x[NR]=$1;y[NR]=$2;ysum+=$2;}\
   END{yave=ysum/NR;
   for(i=1;i<=NR;i++)YY1+=(y[i]-yave)^2;
   pi=atan2(0,-0);st=1;pt=1;en=100;\
   if(arg_st>0&&arg_pt>0&&arg_en>arg_st){
       st=arg_st;pt=arg_pt;en=arg_en;\
   }
   for(f=st;f<=en;f+=pt){
       C=S=YS1=YC1=CC2=SS2=CS2=0;\
       w=2*pi*f;\
       for(i=1;i<=NR;i++){
          wt=w*x[i];coswt=cos(wt);sinwt=sin(wt);\
          C+=coswt;S+=sinwt;\
          YC1+=(y[i]-yave)*coswt;\
          YS1+=(y[i]-yave)*sinwt;\
          CC2+=(coswt)^2;\
          SS2+=(sinwt)^2;\
          CS2+=coswt*sinwt;\
       }
       CC1=CC2-C*C;SS1=SS2-S*S;\
       CS1=CS2-C*S;D=CC1*SS1-CS1^2;\
       p_u=SS1*YC1^2+CC1*YS1^2-2*CS1*YC1*YS1;\
       p_d=YY1*D;if(p_d==0)p_d=1;\
       p=p_u/p_d;\
       print f,p;\
   }
   }' arg_st=$1 arg_pt=$2 arg_en=$3

引数1で開始周波数、引数2で刻み幅、引数3で終了周波数を指定できます。
ここにも簡単に記載されていますが、最小二乗で無理くり求められるため、平均ナイキスト周波数にも制限されません。細かくしすぎると何が原因か分かりませんがめちゃくちゃ大きい値になったりするので、上手く探索するアルゴリズムが必要です。
awkで$π$を表現する方法を初めて知ったのですがpi=atan2(0,-0)なんて書くんですね。

実際に使ってみる

データ１
$sin(2πt)+sin(22πt)+sin(32πt)+sin(42πt)+sin(52πt)$
-10<t<10, n=10000で試してみました。

cat sin1-5.dat | fft 0.1 0.1 6

とりあえず等間隔データに対してはばっちりですね。ちなみにスペクトルは相対密度なので、5本ピークがあってそれぞれ0.2なので正しいです。

この時点でちょっとめんどくさくなったので後で不等間隔データの例についは追記します…
FFTとの差分みたいなものも(1)の論文に書かれたような気がするので見てみてください。
気が向けばもう少し例を追加します。

最後に

事前にどんなデータか踏まえた上で、手動でナイキスト周波数と分解能を決めてあげる必要がありますが、手軽にサクッと確認する分には十分使えそうです。よかったら使ってみてください。
また、改善箇所などあったら指摘していただけると嬉しいです。