コヒーレント状態とウェーブレット変換(3) Cauchy waveletとaffine coherent state -

少し勘違いしてたことなどを補足

coherent state物語
http://d.hatena.ne.jp/m-a-o/20130303#p2

コヒーレント状態の有限群での類似物
http://d.hatena.ne.jp/m-a-o/20130917#p1

が(1)と(2)のつもり

Short-time Fourier変換(STFT)の場合、Gauss窓を使うのは、Weyl-Heisenberg群のコヒーレント状態(boson coherent stateと呼ぶことにする)を与え、この場合は、Gabor変換と呼ばれて、"数学的には"、最も自然なSTFTであると言える。連続ウェーブレット変換(CWT)の場合、背後にある群は、ax+b群であって、この場合は、コヒーレント状態に相当するものがないので、Gabor変換のような特別によいCWTは存在しないのだと思っていた。が、実は、Cauchy waveletというものが存在していて、これは、ax+b群のコヒーレント状態(affine coherent stateと呼ぶ)と考えられているらしい。つまり、

CWT : STFT = ax+b群 : Weyl-Heisenberg群 = Cauchy wavelet : Gabor変換

という対応関係がある。

コヒーレント状態の理論は、コンパクト半単純連結Lie群に対しては、もはや最小不確定性のような特徴づけは出来ないけれど、非常に綺麗に出来ている(コヒーレント状態は、余随伴軌道/旗多様体の小平埋め込みの具体的な実現を与え、それは最高ウェイトベクトルに群を作用させて得られる軌道と一致する)のに対して、一般の非コンパクト群については、少し怪しいものとなる。

非コンパクト群の場合、コヒーレント状態の集合は射影化した表現空間に埋め込むようなことはできないし、Weyl-Heisenberg群の場合、単純に真空ベクトルに群を作用させると、余計な位相因子が付いてきてコヒーレント状態をはみだす。アフィン・コヒーレント状態の場合は、コヒーレント状態の集合というのは一意には決まらない。それでも、あるベクトルの集合がコヒーレント状態と見なされるための強い条件として、幾何学的量子化によって対応する表現が得られる余随伴軌道と同型で、余随伴軌道には自然なケーラー構造が入ってるということがある。Weyl-Heisenberg群の場合、対応するケーラー多様体は複素平面で、ax+b群の場合は、上半平面(z=b+iaによって、座標を入れることができる)となる

#一般的に正しいのかどうか知らないけど、多分正しいこととして、以下のようなことがある。Lie群の表現を微分すればLie環の作用も得られることに注意すると、(コンパクトとも半単純とも限らない)Lig群Gの一般化コヒーレント状態|z>は、GのLie環の元Jに対して、を対応させることで、Lie環の双対空間の元を定める。コヒーレント状態から余随伴軌道への直接的な対応は多分このようにして得られる(つまり、これはコヒーレント状態の集合からある余随伴軌道との1:1対応を与えるのだと思う)。そして、逆にJは、各コヒーレント状態|z>から実数値線形写像を定める。つまり、Jは、コヒーレント状態の空間上の関数とみなすことができ、更にLie環の生成元をJ_1, ... ,J_nとすると、これは、コヒーレント状態の集合(=余随伴軌道)上の関数環の生成元を定める(これは、小平埋め込みとは全く無関係な実アフィン空間への埋め込みを与えると言える)。そして、Lie環の元A,Bに対して、によって、関数とのPoisson括弧が定まる。生成元J_1, ... , J_n間のPoisson括弧が定まれば、一般の関数のPoisson括弧も定まる(勿論、普通≠であるから、関数の積は演算子レベルでは決定できない)。例として、PerelomovによるSU(1,1)の一般化コヒーレント状態|z,k>に対して(今の場合、コヒーレント状態系の点zは、|z|<1で単位円板上の点に対応し、kは正の整数or半整数で、正則離散系列表現を指定するindex)、K_0,K_1,K_2をLie環su(1,1)の(標準的な)生成元とすると
$y_0 = \langle z,k|K_0|z,k \rangle = k \frac{1+|z|^2}{1-|z|^2}$
$y_1 = \langle z,k|K_1|z,k \rangle = k \frac{2 Im(z)}{1-|z|^2}$
$y_2 = \langle z,k|K_2|z,k \rangle = k \frac{2 Re(z)}{1-|z|^2}$
となって、
$y_0^2-y_1^2-y_2^2 = k^2$
なので、座標(y_0,y_1,y_2)は二葉双曲面の連結成分と対応し(y_0>0に注意)、単位円板と同型である。更に、この余随伴軌道の幾何学的量子化によって、元の既約表現を復元できる(Kostant-Kirillovの軌道法)ので、コヒーレント状態は"最も古典的な状態"とか言われるけど、実際にある意味で量子化の逆を与える(まぁ、物理としては、状態のHilbert空間が、適当なLie群の(二乗可積分な)既約表現空間になっているというのは結構特殊な状況である気もするけど)。関数環の"量子化"を得たければ、普遍展開環の既約表現のannihilator(原始イデアルと呼ぶ)による商環を考えればよい？記号とか計算の詳細は、Perelomovの原論文
Coherent states for arbitrary Lie group
http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.cmp/1103858078
や
Some basics of su(1,1)
http://arxiv.org/abs/quant-ph/0407092
を参照。別のよく知られた一般化コヒーレント状態として、su(2)コヒーレント状態(スピン・コヒーレント状態)があるけども、この場合は、余随伴軌道は球面と同型で、全ての有限次元既約表現は、余随伴軌道の幾何学的量子化で得られる
geometric quantization of the 2-sphere
http://ncatlab.org/nlab/show/geometric+quantization+of+the+2-sphere

Cauchy waveletの具体的な形は
Time-frequency localisation operators - a geometric phase space approach: II. The use of dilations
http://www.math.duke.edu/~ingrid/publications/Inverse_Prob_4_p661.pdf
のAppendixに書いてある(ここに書いてあるのは、より一般的な何かで、特にγ=1の時が、Cauchy waveletで、βごとに異なるコヒーレント状態系が得られる)。Mother waveletの形は、Poisson分布の確率密度関数を思い出すけども、Poisson分布の時と、変数と定数が入れ替わってる

歴史的には、Klauderがaffine coherent stateの考案者らしい。
Path Integrals for Affine Variables
http://link.springer.com/chapter/10.1007%2F978-1-4615-7035-6_7

Unitary Representations of the Affine Group
http://jmp.aip.org/resource/1/jmapaq/v9/i2/p206_s1?isAuthorized=no

これは、Morletによるウェーブレット解析の発表より早い。Cauchy waveletという命名は、Holschneiderという人の
Wavelets: An Analysis Tool
http://www.amazon.co.jp/dp/0198505213
という本によるらしい(1995年)。何でCauchyなのかよく分からないのだけど、この名前は一般的になっているっぽい。

#別のアプローチとしてboson coherent stateの最小不確定性状態という定義を直接+b群に持ち込もうという試みが
The Affine Uncertainty Principle in One and Two Dimensions
http://www.sciencedirect.com/science/article/pii/0898122195001085
でなされているけど、多分あまりよい答えではないと、個人的には思う

ax+b群は、SU(1,1)に埋め込むことができ、ax+b群の既約ユニタリ表現は、SU(1,1)の離散系列表現をax+b群に制限することで得られる。Cauchy waveletとPerelomovによるSU(1,1)コヒーレント状態は実は、適当な関数空間で見ると、同じ形をしている
Characterization of SU(1,1) coherent states in terms of affine group wavelets
http://arxiv.org/abs/math-ph/0209006

それで数学的なことは、ともかく、実際にCauchy waveletはMorlet waveletとかよりよかったりするのかというと、謎。あんまり使ってる人いないし。

おまけ:というか、STFTでも、Gabor変換と他の窓関数使った場合で、どういう差があるのだろう
窓関数
http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/basic/chap9/index.htm
とか見ると、ガウス窓は、メインローブの幅がでかいけど、サイドローブの振幅が小さくて、"主成分の周波数分解能は劣りますが、小さい電力のスペクトルを検出するのに優れた特性を示します"とか書いてある。

最初何を言ってるか分からなかったけれど、メインローブの幅というのは、要するに単純な正弦波をSTFTした時の分散の大きさのことっぽい。これは矩形窓でやるのがベストに決まっていて、十分に大きな矩形窓を取れば、無限区間でFourier変換するのとほぼ同じ結果(デルタ関数的な何か)に近づいていく。完全に周期的な波では、サイドローブとか出てこないのだけど、現実の音とか信号は、減衰していくのが普通なので、これがサイドローブとか生み出す要因になるっぽい(直接的には、矩形関数のFourier変換がsinc関数なので、それがサイドローブの元であるらしい)。サイドローブがあると、複数個所にピークもどきが出てきてしまうので、あまり嬉しくない

Gauss窓だと、そんなにメインローブの幅が大きくなるのだろうかと思って、適当に実験して見た

#!/usr/bin/env python
import numpy as np
import scipy as sp
from scipy import signal
from math import sqrt,pi,sin

"""
import wave

def wavread(filename):
   wf = wave.open(filename , "r")
   fs = wf.getframerate()
   x = wf.readframes(wf.getnframes())
   if wf.getsampwidth()==2:
      x = np.fromstring(x,dtype=np.int16)/32768.0
   else:
      x = np.fromstring(x,dtype=np.int8)/128.0
   wf.close()
   return x,fs
"""

fs = 44100
sinewave = np.array([sin(2*pi*n*400/fs) for n in xrange(fs*2)])

window_func = {"Hamming" : (lambda sz:sp.hamming(sz)) ,
               "Hanning" : (lambda sz:sp.hanning(sz)) ,
               "Rectangular" : (lambda sz:sp.kaiser(sz,0)) ,
               "Kaiser(alpha=3)" : (lambda sz:sp.kaiser(sz,3)) ,
               "Gaussian(window size=2 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/2.0)/(sqrt(2*pi)*sz/2.0)) ,
               "Gaussian(window size=4 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/4.0)/(sqrt(2*pi)*sz/4.0)) ,
               "Gaussian(window size=5 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/5.0)/(sqrt(2*pi)*sz/5.0)) ,
               "Gaussian(window size=7 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/7.0)/(sqrt(2*pi)*sz/7.0)) ,
               "Gaussian(window size=10 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/10.0)/(sqrt(2*pi)*sz/10.0)) ,
               "Gaussian(window size=20 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/20.0)/(sqrt(2*pi)*sz/20.0))}

if __name__=="__main__":
   for (name , wfun) in window_func.iteritems():
      start = int(0.5*fs)
      framesize = int(0.04*fs)  #-- 40msec
      w = wfun(framesize)
      xs = sp.fft(w*sinewave[start:start+framesize])
      total = sum([abs(v)**2 for v in xs])
      avg = sum([min(n,framesize-n)*abs(v)**2 for (n,v) in enumerate(xs)])/total
      delta = sum([(min(n,framesize-n)-avg)**2*abs(v)**2 for (n,v) in enumerate(xs)])/total
      print "avg=%f[Hz] delta=%f name=%s" % (avg*fs/framesize,delta,name)

やってることは、400Hzの正弦波に色んな窓関数でSTFTしてるだけ(窓幅は40ms固定)。出力は

avg=399.999999[Hz] delta=0.333712 name=Hanning
avg=400.000949[Hz] delta=0.193046 name=Gaussian(window size=4 sigma)
avg=400.000000[Hz] delta=1.266515 name=Gaussian(window size=10 sigma)
avg=400.002197[Hz] delta=0.137814 name=Kaiser(alpha=3)
avg=399.999999[Hz] delta=0.266531 name=Hamming
avg=400.000000[Hz] delta=5.066059 name=Gaussian(window size=20 sigma)
avg=400.000005[Hz] delta=0.620575 name=Gaussian(window size=7 sigma)
avg=400.000658[Hz] delta=0.024861 name=Gaussian(window size=2 sigma)
avg=400.000318[Hz] delta=0.314517 name=Gaussian(window size=5 sigma)
avg=400.000000[Hz] delta=0.000000 name=Rectangular

とかなって、deltaが分散の大きさを表す(平方根とって適当にスケールすれば周波数の単位に変換できる)けど、正弦波相手だと、圧倒的に矩形窓最強っぽい。次によいのが、Gaussian(window size=2 sigma)というやつだけど、例えばsignal.gaussian(10,std=5)はarray([ 0.66697681, 0.78270454, 0.8824969 , 0.95599748, 0.99501248,0.99501248, 0.95599748, 0.8824969 , 0.78270454, 0.66697681])とかで、窓の端の値がピーク値の0.6倍くらいもあるので、ガウス窓といっても、積分区間が明らかに不十分。

signal.gaussian(10,std=2)だとarray([ 0.07955951, 0.21626517, 0.45783336, 0.7548396 , 0.96923323, 0.96923323, 0.7548396 , 0.45783336, 0.21626517, 0.07955951])で、scipy.hamming(10)がarray([ 0.08 , 0.18761956, 0.46012184, 0.77 , 0.97225861,0.97225861, 0.77 , 0.46012184, 0.18761956, 0.08 ])なので、これくらいだと、Hamming窓に近い感じがある。signal.gaussian(10,std=10.0/7)はarray([ 0.00700417, 0.04972487, 0.21626517, 0.57622907, 0.94058806, 0.94058806, 0.57622907, 0.21626517, 0.04972487, 0.00700417])でよさげではあるけど、このくらいになると、分散が他の窓関数より大きい

なんか別にGabor変換が殊更よいというわけでもないっぽい。。まぁ、そもそも何故Weyl-Heisenberg群やax+b群を選ぶのかという問いにも明確な答えがあるわけでもないと思うし