yProcessingClub

すみません、許してください

Signal Kernelの計算方法

信号処理


66-67ページから解説．


ある時間信号z(t)の時間-周波数表現\rho_z(t,f)があるとして，\rho_z(t,f)を逆フーリエ変換したものをK_z(t,\tau)と書く．ここでK_z(t,\tau)をSignal Kernelと呼ぶことにする．K_z(t,\tau)をフーリエ変換すると\rho_z(t,f)になるわけだから，数式で書くと，
\rho_z(t,f) = \int K_z(t,\tau) \mathrm{e}^{-j 2 \pi f \tau} d\tau


ここでK_z(t,\tau)はどのような形になるか．

ある信号z(t)の例として，解析信号z(t)=\mathrm{e}^{j \phi (t)}を考える．

瞬時周波数f_i (t) = \frac{\phi ' (t)}{2 \pi}より，

\rho_z(t,f) = \delta (f - f_i(t))である．

例えば5秒時点において100Hzの部分にスペクトルがあるとすれば\rho_z(5,100) = 1となるわけである．

K_z(t,\tau)は\rho_z(t,f) = \delta (f - f_i(t))の逆フーリエ変換なので，
\begin{eqnarray*} K_z(t,\tau) &=& \int \rho_z(t,f) \mathrm{e}^{j 2 \pi f \tau} df\\ &=& \int \delta (f - f_i(t)) \mathrm{e}^{j 2 \pi f \tau} df\\ &=& \mathrm{e}^{j 2 \pi f_i(t) \tau}\\ &=& \mathrm{e}^{j \phi ' (t) \tau}\\ \end{eqnarray*}

ここの計算は \mathcal{F}^{-1} [ \delta (f - f_0) ] = \mathrm{e}^{j 2 \pi f_0 \tau}を使えば行ける．

さて，瞬時周波数\phi ' (t)だが，
\displaystyle \phi ' (t) = \lim _{\tau \to 0} \frac{ \phi(t + \tau / 2) - \phi(t - \tau / 2)}{\tau}
と書ける．これは中心有限差分などで検索すると色々出てくるが，ある点 \phi(t)における傾き（ =\phi ' (t)）は， \phi(t)の前後の2点 \phi(t + \tau / 2)及び\phi(t - \tau / 2)を通る直線の傾きに等しい（\tau \to 0において）と言える．ここで\phi ' (t) が直線ならば、言い換えると\phi ' (t) が一次式（または定数）ならば，さらに言い換えると \bf {\phi (t)}が二次式（または一次式）ならば，\tau \to 0でなくても傾きは等しくなる．つまり，

\phi ' (t) \approx \frac{ \phi(t + \tau / 2) - \phi(t - \tau / 2)}{\tau} \ \ (\mathrm{if} \ \ \phi (t) \mathrm{\ \ is \ \ quadratic \ \ or \ \ linear})

これを代入すると，
\begin{eqnarray*} K_z(t,\tau) &=& \mathrm{e}^{j \phi ' (t) \tau}\\ &=& \mathrm{e}^{j(\phi(t + \tau / 2) - \phi(t - \tau / 2))}\\ &=& \mathrm{e}^{j\phi(t + \tau / 2)}\mathrm{e}^{-j\phi(t + \tau / 2)}\\ &=& z(t+\tau / 2)z^*(t-\tau / 2) \end{eqnarray*}


よって，
\begin{eqnarray*} \rho_z(t,f) &=& \int K_z(t,\tau) \mathrm{e}^{-j 2 \pi f \tau} d\tau\\ &=& \int z(t+\tau / 2)z^*(t-\tau / 2) \mathrm{e}^{-j 2 \pi f \tau} d\tau\\ \end{eqnarray*}

この形はなんじゃらほいと考えると，まさしくウィグナービレ分布なわけである．

くどいようだが，この話が正しくなる条件として\phi ' (t) が直線である必要がある．そうでなければ極限の近似が上手くいかないからである．

s(t)=cos(2πft)のウィグナービレ分布を求めよ

信号処理

ウィグナービレ分布とは？

ウィグナービレ分布とは時間周波数解析法の一種である．「解析信号z(t)=\mathcal{A}(s(t))のウィグナー分布」をウィグナービレ分布といい，ウィグナー分布と比較してクロス項の抑制が期待できる．

問題

s(t)=\cos(2 \pi f_c t)のウィグナービレ分布を求めよ．


解法

まず解析信号z(t)を求める．
オイラーの公式より，
s(t)=\cos(2 \pi f_c t) = \frac{1}{2} \left( \mathrm{e}^{j 2 \pi f t} + \mathrm{e}^{-j 2 \pi f t} \right)

また，\mathrm{e}^{j 2 \pi ft}の解析信号は
\mathcal{A}(\mathrm{e}^{j 2 \pi ft})= \begin{cases} 0 \ \ (\mathrm{if} \ \ f < 0)\\ 2 \mathrm{e}^{j 2 \pi ft} \ \ (\mathrm{if} \ \ f > 0)\\ \end{cases}
であるので，

\begin{eqnarray*} z(t)&=&\mathcal{A}(s)=\frac{1}{2} \mathcal{A}(\mathrm{e}^{j 2 \pi f_c t})+\frac{1}{2} \mathcal{A}(\mathrm{e}^{-j 2 \pi f_c t})=\frac{1}{2} 2 \mathrm{e}^{j 2 \pi f_c t}+\frac{1}{2}0\\ &=&\mathrm{e}^{j 2 \pi f_c t} \end{eqnarray*}


また，z^*(t)=\mathrm{e}^{-j 2 \pi f_c t}．
よって，
\begin{eqnarray*} z\left(t+\frac{\tau}{2}\right)z^*\left(t-\frac{\tau}{2}\right) &=& \mathrm{e}^{j 2 \pi f_c \left( t + \frac{\tau}{2} \right)} \mathrm{e}^{-j 2 \pi f_c \left( t - \frac{\tau}{2} \right)}\\ &=& \mathrm{e}^{j 2 \pi f_c t} \mathrm{e}^{j 2 \pi f_c \frac{\tau}{2}} \mathrm{e}^{-j 2 \pi f_c t} \mathrm{e}^{j 2 \pi f_c \frac{\tau}{2}}\\ &=& \mathrm{e}^{j 2 \pi f_c \tau} \end{eqnarray*}



次にz(t)のウィグナー分布（＝ウィグナービレ分布）W_z(t, f)を計算する．
\begin{eqnarray*} W_z(t, f)&=&\int z\left(t+\frac{\tau}{2}\right)z^*\left(t-\frac{\tau}{2}\right) \mathrm{e}^{-j 2 \pi f \tau}d\tau\\ &=& \int \mathrm{e}^{j 2 \pi f_c \tau} \mathrm{e}^{-j 2 \pi f \tau}d\tau\\ &=& \delta(f-f_c) \end{eqnarray*}
となり，計算終了である．


答え

s(t)=\cos(2 \pi f_c t)のウィグナービレ分布W_z(t,f)は，
W_z(t,f) = \delta(f-f_c)


クロス項抑制

解析信号は周波数f_cの信号で負の周波数を持たないものであり，そのウィグナービレ分布はW_z(t,f) = \delta(f-f_c)と期待通りのスペクトルを持つ．ウィグナー分布ではf=\pm f_cの中間のf=0の部分にクロス項が出ていたが，ウィグナー分布は信号を解析信号に変換することで負のスペクトルを消したため，2信号の中間に出るクロス項の発生を抑制することができた．



yuri-processing-club.hatenablog.com

s(t)=cos(2πft)のウィグナー分布を求めよ

信号処理

ウィグナー分布とは？

ウィグナー分布とは時間周波数解析法の一種である．

問題

s(t)=\cos(2 \pi f_c t)のウィグナー分布を求めよ．

解法

ウィグナー分布は以下で定義される．
W_s(t,f)=\int s\left(t+\frac{\tau}{2}\right)s^*\left(t-\frac{\tau}{2}\right) \mathrm{e}^{-j 2 \pi f \tau}d\tau


s(t)は実信号なので，s^*(t)=s(t)．
よって，
s\left(t+\frac{\tau}{2}\right)s^*\left(t-\frac{\tau}{2}\right) = s\left(t+\frac{\tau}{2}\right)s\left(t-\frac{\tau}{2}\right)

s(t)=\cos(2 \pi f_c t)を代入して，
s\left(t+\frac{\tau}{2}\right)s\left(t-\frac{\tau}{2}\right) = \cos 2 \pi f_c \left(t+\frac{\tau}{2}\right) \cos 2 \pi f_c \left(t-\frac{\tau}{2}\right)


ここで\alpha = 2 \pi f_c \left( t + \frac{\tau}{2} \right)，\beta = 2 \pi f_c \left( t - \frac{\tau}{2} \right)とおくと，
\alpha - \beta = 2 \pi f_c \tau，\alpha + \beta = 2 \pi f_c 2t
加法定理： \cos \alpha \cos \beta = \frac{\cos(\alpha - \beta) + \cos(\alpha + \beta)}{2}を用いて，

\cos 2 \pi f_c \left(t+\frac{\tau}{2}\right) \cos 2 \pi f_c \left(t-\frac{\tau}{2}\right) = \frac{1}{2}cos ( 2 \pi f_c \tau ) + \frac{1}{2}cos ( 2 \pi f_c 2t )


これをW_s(t,f)=\int s\left(t+\frac{\tau}{2}\right)s^*\left(t-\frac{\tau}{2}\right)\mathrm{e}^{-j 2 \pi f \tau}d\tauに代入して，

\begin{eqnarray*} W_s(t,f)&=&\int {\left( \frac{1}{2}\cos ( 2 \pi f_c \tau ) + \frac{1}{2}\cos ( 2 \pi f_c 2t ) \right) \mathrm{e}^{-j 2 \pi f \tau}d\tau}\\ &=& \frac{1}{2} \int \cos ( 2 \pi f_c \tau )\mathrm{e}^{-j 2 \pi f \tau}d\tau + \frac{1}{2}\cos ( 2 \pi f_c 2t ) \int \mathrm{e}^{-j 2 \pi f \tau}d\tau \end{eqnarray*}

ここで\cos(2 \pi f_c t)のフーリエ変換は \frac{\delta(f-f_c) + \delta(f+f_c)}{2}
1のフーリエ変換は\delta(f)より，

\begin{eqnarray*} W_s(t,f)&=& \frac{1}{2}\frac{\delta(f-f_c) + \delta(f+f_c)}{2} + \frac{1}{2}\cos ( 2 \pi f_c 2t )\delta(f)\\ &=& \frac{1}{4}\delta(f-f_c) + \frac{1}{4}\delta(f+f_c) + \frac{1}{2}\cos ( 2 \pi f_c 2t )\delta(f) \end{eqnarray*} 　
となり，計算終了である．


答え

s(t)=\cos(2 \pi f_c t)のウィグナー分布W_s(t,f)は，
W_s(t,f) = \frac{1}{4}\delta(f-f_c) + \frac{1}{4}\delta(f+f_c) + \frac{1}{2}\cos ( 2 \pi f_c 2t )\delta(f)


クロス項

周波数f_cの信号であるので，両側スペクトルを考えると，f=\pm f_cにスペクトルがあるべきで，\frac{1}{4}\delta(f-f_c)及び\frac{1}{4}\delta(f+f_c)の項は期待通りである．しかし，\frac{1}{2}\cos ( 2 \pi f_c 2t )\delta(f)という謎の項が出ている．これをクロス項といい，アーティファクトである．クロス項の発生個所を見ると，f=0の位置，つまりf=f_c及びf=-f_cの中間の部分にスペクトルが出ていることが分かる．このように，クロス項は信号と信号の中間の位置に出る．



yuri-processing-club.hatenablog.com

防衛大 開校祭に行ってきた

防衛

防衛大の開校祭（学園祭）に行ってきた。

f:id:Yuri-Processing-Club:20181111231934p:plain
訓練展示。
敵と味方に分かれて、陣地を占領するという訓練を行った。FH-70などが支援射撃を行ったりして大迫力であった。


f:id:Yuri-Processing-Club:20181111225414j:plain
売店（PXまたはPと呼ぶらしい。免税店ではないが・・・？）にて。
「行軍に最適」とはこういうところでしか見られないキャッチコピーだろう。


PAC3

f:id:Yuri-Processing-Club:20181111224806j:plain
ペトリ。

f:id:Yuri-Processing-Club:20181111224745j:plain

f:id:Yuri-Processing-Club:20181111225916j:plain
フェーズドアレイレーダー。

f:id:Yuri-Processing-Club:20181111225731j:plain
ミサイルだけじゃなくて、電源車とか一式来ていた。

FH-70

f:id:Yuri-Processing-Club:20181111224914j:plain

f:id:Yuri-Processing-Club:20181111224854j:plain

f:id:Yuri-Processing-Club:20181111224839j:plain

f:id:Yuri-Processing-Club:20181111224829j:plain
砲脚ブレード。地面にぶっ刺さっており、発射の衝撃を受け止める。

f:id:Yuri-Processing-Club:20181111224822j:plain
写真だとなかなか伝わらないが、戦車の砲身より一回り大きくてとても迫力があった。

10式戦車

f:id:Yuri-Processing-Club:20181111225055j:plain

f:id:Yuri-Processing-Club:20181111225046j:plain

高機動車

f:id:Yuri-Processing-Club:20181111225008j:plain
初めて生で見たのだが、めちゃくちゃでかくて笑ってしまった。
ランクルみたいな大きさを想定していたが、予想以上だった。

f:id:Yuri-Processing-Club:20181111224947j:plain

f:id:Yuri-Processing-Club:20181111224939j:plain
後ろから。
運転席と助手席の間が相当空いていることが分かると思う。そのくらい横幅がでかいのである。

カワサキ

f:id:Yuri-Processing-Club:20181111224924j:plain



反省点

今回は儀仗隊及び訓練展示の動画撮影を行った。
動画撮影はほとんど初めてだったが、いくつか足りない部分があった。

・AF-Sで撮影した。
AFの切り替えを忘れてしまった。動画撮影前にAF-Cに切り替える必要がある。

・手振れ対策が何もない
レンズまたはカメラに手振れ補正がなく、三脚やスタビライザも持っていなかった。そのため動画が手振れでぶれぶれになってしまった。（そのため、ブログでの公開は断念）

・レンズが広角で寄れない
FH-70の射撃シーンを望遠レンズで撮りたかったが、今回は23mmの単焦点しか持ってなかったのでだめだった。こういうイベントの撮影だと、全体を写す広角レンズつけたカメラ（これは三脚でガチガチに固定しておいて、撮影開始したら一切動かさない）と、要所要所をアップで撮る望遠レンズをつけたカメラ（こっちは自分で構えて撮る）の最低２台は必要である。

・バッテリー切れ
動画撮影は物凄く電気を消費するのか、一瞬でバッテリーを使い果たしてしまった。
そのため、動画を撮影した後はバッテリー切れを心配して満足に写真を撮れなかった。
予備のバッテリーは必須である。

・風切り音対策をしていない
録音は意外と問題なかった。自分が使っているカメラの内蔵マイクは割と優秀であった。ただし、今回は無風のため上手く行ったのであって、強風時だと「ボボボ」という風切り音みたいなのが入りまくるはずである。ウィンドガードを装着すべきであった。

こういうやつね。


以上の反省点をフィードバックし（つまり冬のボーナスをカメラ機材に溶かす）、次回に生かす。

Cisco 1812J購入

ネットワーク

CCNAの勉強のために、中古のCisco 1812Jを買ってみた。 値段は5千円ほど。

f:id:Yuri-Processing-Club:20181110232448j:plain

f:id:Yuri-Processing-Club:20181111000019j:plain 外見はこんな感じ。


f:id:Yuri-Processing-Club:20181111000727j:plain PCからコンソールを叩くためのケーブルも買った。


f:id:Yuri-Processing-Club:20181111001356p:plain コンソールは無事表示できた。



今後これを使ってネットワークの勉強をしていく。



余談。

購入場所は秋葉原のNW工房というお店なのだが、入店する前はとてもとても怖かった。外から店内の様子を見たところ秋葉の古くからあるマニアックなお店という感じで、入るには敷居が高く思ったのである。勇気を出して入店し、そして曖昧に相談してみたところ、とても親切に対応していただいた。

僕「会社でCCNAの勉強しろって言われてルータを買ってみようと思うのですが、オススメのやつあります？」
店員さん「希望の機種とかあります？」
僕「まったくわかりません」
店員さん「じゃあこれとかどうですか？」
こんな感じ。

僕のようなズブの素人はボコボコに殴られて「ヨドバシカメラへ帰るんだな」って言われるかと思ってたので、非常にありがたい対応であった。





さらに余談だが、 yuri-processing-club.hatenablog.com 自分のこの記事を見ると、2年前にも「CCENTをとるために勉強してる」って記事を書いていた。そこから2年経ったがネットワークの知識はまるで増えていない。まぁ仕事ではあんまりネットワーク触らないからね・・・・・・。こんな感じだと10年後にも「今年こそネットワーク理解するぞ」って言ってそう。

今回は機材も買ったし、しっかりと勉強していきたい。

メモ：連続と離散の変換　積分

信号処理

数式の離散表現

　数学や信号処理をプログラミングで実装する際，注意すべきは連続と離散の取り違えである．数学などの専門書においては数式は連続型であるが，プログラミングする場合は離散型に適度に読み変える必要がある．今回は簡単に積分の離散表現を示す．

結論

y=s(t)の積分は以下．

Y=\int s(t)dt


これを離散的に表現すると，

y=s(nT_s)の積分は以下．

Y=\Sigma _{n=0} ^{N-1} T_s s(nT_s)

ここでT_sはサンプリング周期．

くわしい説明

y=s(t)はこういう形だとする．
f:id:Yuri-Processing-Club:20181108182004p:plain

これを時間間隔T_sでサンプリングしたのがs(nT_s)である．
f:id:Yuri-Processing-Club:20181108182223p:plain


Y=\int s(t)dtはグラフの面積を表している．

f:id:Yuri-Processing-Club:20181108182445p:plain

これを離散的に表現すると，T_s s(nT_s)で表される長方形の面積の合計値で近似することになる．

f:id:Yuri-Processing-Club:20181108182617p:plain

終わり

仕事でここを取り違えて大ポカをやらかしてしまったので，忘れないようにブログにてメモした．

瞬時周波数についてのふわっとした解釈

信号処理

記憶の整理のためにメモ．
記事の正しさに責任は持ちません．

追記。
↓qiitaに記事を投稿した。ブログ記事の焼き直しだが、より詳しく書いてるので良かったら見て欲しい。
qiita.com


瞬時周波数とは？

ある信号（音響信号など）の周波数の時間変化（変調と言ったりする）を調べたい欲求がある人もいるだろうが，そういうのを調べるのを時間周波数解析という．時間周波数解析において一番知りたいのは「（目的の信号の）周波数の時間変化」であり，これは位相を時間微分することで求められ，これを瞬時周波数と呼ぶことにする．（速度の時間変化は速度を微分すると求められ，これはどうやら加速度と言うらしいのだが，位相を微分したものは瞬時周波数と言う．）


簡単にまとめると，

ある一瞬の周波数を知りたい
　↓
位相を時間微分する

という素朴な考え方だ．


　　

周波数が一定の信号

周波数が時間変化する信号を考える前に，
まず周波数が一定の信号を考える．

振幅1，周波数f_0（一定）である信号s(t)
s(t)=\cos(2 \pi f_0 t)
ここで位相 2 \pi f_0 tを微分すると2 \pi f_0であり，瞬時周波数は時間変化せず、周波数が一定であることを再確認できる．

周波数が時間変化する信号

次に周波数が時間変化する信号を考える．
瞬時周波数をf'(t)で表す．

t=0[\mathrm s]でf'(0)=1000[\mathrm {Hz}]，
t=10[\mathrm s]でf'(10)=5000[\mathrm {Hz}]，
と周波数が線形変化する場合，
瞬時周波数の時間-周波数グラフは以下のようになる．

このグラフの式は
\begin{eqnarray*} f'(t)&=&\frac{5000-1000}{10-0}t + 1000\\ &=&400t+1000 \end{eqnarray*}

周波数の時間変化を見るには位相を微分すれば良かったので，逆向きに考えると，このような周波数の時間変化を持つ位相はf'(t)を積分すれば求められる．

\begin{eqnarray*} f(t)&=&\int f'(t)dt\\ &=&\int (400t+1000)dt\\ &=&200t^2+1000t \end{eqnarray*}

これを\cos()の中に入れると，
周波数が時間変化する信号
s(t)=\cos(2\pi(200t^2+1000t))
が出来上がる．
このスペクトログラムを求めると，

となり，期待通りの信号となった．

ちなみに，このように線形に周波数が変調する信号をLFM信号(Linear Frequency Modulated Signal，線形周波数変調信号)とか線形チャープ信号と言ったりする．

別の例

瞬時周波数が
g'(t) = 500\sin(2 \pi t) + 500t + 3000
となる信号を考える．
g'(t)のグラフは以下．

これもg'(t)を積分して，
g(t) = \int g'(t)dt = -\frac{500}{2 \pi} \cos(2 \pi t) + 250t^2 + 3000t
であり，これを\cos()に入れて，
信号s(t)=\cos\left( 2 \pi \left( -\frac{500}{2 \pi} \cos(2 \pi t) + 250t^2 + 3000t \right) \right)が出来上がる．
このスペクトログラムは

となって期待通りである．