一次元離散

一次元離散

分布の計算方法

鈴 木 宏

A   N e w A l g o r i t h m   t o   C o m p u t e   t h e   D i s c r e t e   W i g n e r   D i s t r i b u t i o n Hiroshi SUZUKI

WignerDistribution(WD)hasrecentlyattractedattentionasatoolfortime‑ frequencysignalanalysis.However,itsprocessingtimeneedsseveraltimesasmuchas Fouriertransfom (FT)that.Thispaper丘rstconsidersthatWDcanbecalculatedasthe convlutionofFourierspectra(FS).Thenwepresentafastcomputationmethodofthe discreteWD.Finally,wecomparethistechniquewithseveralmethods.

1 は じ め に

最近,時変信号の有力 な解析法の一つ として,Wigner分布 (以下WD)が注 目されてい る(1)‑(2). これは,特殊な前処理を行 った信号の フー リエスペク トル (以下FS)で, この処 理により周波数分解能が, フー リエ変換のそれに比べ優れている.反面 この前処理に よ り, 複数の周波数成分を含む信号で,その成分間に干渉項が生 じスペク トルを汚す. この うち負 周波数成分 との干渉項 お よびエイ リアシンダの影響を取除 くために,一般 にWDは解析信 号を用いて計算す る.解析信号を計算す るには2回の 7‑ リエ変換 (以下FT)が必要 とな

り,最終的にWD を求めるには3回のFTが必要 となる.

WDは,FTに比べ前処理などがあるために,計算時間がかかる.そ こで,い くつかの高 速浜算法 ･計算法(4)‑(5)が考 えられてい る.それ らは,前処理 した信号の対称性 を利 用 した

cos変換(6),信号が実部のみ よ り1/2点FTの利用(7),‑ ‑ ドウェアによる手法(5)な どがある が,解析信号のWDに対する高速計算手法の報告は少ない.

本稿では,一般的なWDの計算手法を述べた後 に,解析信号 のWDがFSの畳込み によ り計算 される⁽³⁾ことを従来のWD と比較 して確認 し,それを利用 した計算手法を提案す る.

この手法は, 0を掛 ける操作を省 くことにより,解析点数 (ウィソ ドゥ点数)が少 なけれ ば, WDの計算時間が短縮で きる.最後に,FT,一般的 なWD と畳込みによるWD との計算 時間を比較 し,本手法が従来の方法に比べ優れていることを示す.

2

分布の一般的計算方法

68 _{鈴 木 宏}

〟/2 〟/2

W(nh･ka'o)=m=‑∑M/2I(to(n

+

u o

( 二 義

で定 義 され て い る(1).す なわ ち,WDは,あ る時 刻nもに お い て,(n+m)も時 の信 号 と (n‑m)も時 の信号 の複素共役 とを掛合わ せたそのFSであ る.一般 に観測信号 は実関数で あるため に,図1に示す よ うに,単 に信号 のA部 とB部 とを掛合わせ て, このmをm‑‑

M/2‑M/2まで変化 させ,図2のよ うなnも時 で対称 となる信号R(ち,n, m)を求 め, この信号 のFSがWDとなる.

(n‑m)lo nlo (n+m)10 '図1 観測信号

f (

nto

図2 信 号 R(ち,n,m)

FTす る信号R(^ち,n,m)が対称 (偶 関数) である ことよ り,WD の虚部 はすべ て 0と なる. このために,奇関数成分 のFSを求 め る必要 がな くcos変換(6)(FSの実部 の成分 しか 求め ない変換)が利用で きる. また,FTす る信号尺(ち, n, m) が虚部 を持 た ないため に, 奇数項 を虚部 に入れて半分 の点数で計算す る高速演算法(7)が可能 とな る. しか し後者 は,負 周波数成分 との干渉項除去 の立場 か ら実際 には利用 され ていない.

WDは,負周波数成分 との干渉項 お よびェイ リアシンダの影響があるので, これ らを無 く すために,観測信号 その ものを解析 す るのではな く解析信号,

fa(nち)‑I(nち)+

l f h (

を用 いて計算す る. これ は,実 部 に観 測 信号 U(nち)) を虚 部 に観 測信 号 をHilbert変換 (位相 を900シフ トさせ る効果 を もつ変換) した信号 Uh(nち)) をそれぞれ持 つ ものであ る.

この fh(nち) 紘,

〟/2

fh(nゐ)=m

=

Sin2(^翠 ) Miq2

‑ ∑ f((n‑ m)ち)h(mゐ)

m=‑N/2 pziLo

2

ここで,m≠0,h(mち) はHilbert変換器 のインパルス応答

で定 義 されてい る(2).これ は,観測信号f(nち) と変換器 の イ ソパ ル ス応 答h(mち) の畳込 み となってい るために,FTす ることに よ り掛 算 で計算 で きる.す なわ ちfh(nち)のFS(Fh

(kaJo)) は,I(nち) のFSをF(kaJo) とす ると,

Fh(ka)o)ニーjsgnkF(kaJo) (4)

ここで,

s g n k ‑舟

これ よ り,解析信 号fa(nち)のFS(Fa(kalo)) は, Fa(kaJo)‑F(kaJo)+jFh(ka)o)

‑F(ka)o)+i(‑isgnkF(kaJo))

‑F(kwo)+sankF(kaJo)

‑(1+sgnk)F(kwo) (5) となる.すなわち,観測信号f(nh)をFTした後 に,正周波数成分 (k>0のFS)を2倍,負周 波 数 成 分 (k<0のFS)を0,直 流 成 分 (k‑ 0の FS)とナ イキス ト周波数成 分 (k‑M/2のFS) を 1倍 にそ れ ぞ れ し,そ の道 FTが解 析 信 号fa (nん) となる.

以上 よ り,解析信号 のWDを求め るためには, 図 3に示す よ うに,観測信号をFTL式(5)に よ り Fa(kaJo)を求 め逆FTを して,式(1)よ りR(ち,n,

m)を計算 し,そのFSがWDとなる.

したが って解析信号 のWDを求め るには, 3回 のFTが必要 となる.

図3 解析信号のWigner分布 の一般的計算法

3 Wigner

分布の畳込み表現 とその意味

nh時を恥 とす るウィン ドウ長 M の観潮信号f((n･m)ち)(ここで,‑孝 <‑ ≦# ) は,直流成分を持たない として,

〟/ 2

f (

i .

‑∑[ L l ⇒ i l

.

+

TVa(nt.,kw^{.) ‑ "i 2 "i 2[(A}pAq+ BpB.)8[(k‑(♪+ q))a,.]] (7)

Ll=lq=1

ここで,SlkaJ.]‑ 1(k=0),0(k≠0)

で計算す ることができる(3).すなわち,解析信号のWDはFSの実部 と虚部 の各々の畳込み の和で求めることができる.

式(7)紘,図3のC部の計算を表 している.これは,逆FTL変則的は掛算を した後FTす る計算 と,式(7)による畳込み計算が等 しいことを意味す る. これ らの等価性については付録

70 _{鈴 木 宏}

4

分布の畳込みによる計算手法

従来WDは図3の手法 (以下一般的WD:GWD)で計算 されていた. ここでは,式(7)に よるWD(以下畳込 みWD:CWD) を考 える.一般 的に,畳込み計算 は非常 に時間がかか るため,信号点数が2のn乗 の場合は,信号をFTして周波数上で掛算 して逆FTす る. こ れは,式(7)と図3の関係 と同 じである.

CWDである式(7)を行列表現す ると,(〟‑〟/2とす る)

A

A

･･･ ･･

‑

4'二..<･.'Al 0‑‑OAA.+++.I 01

0

0

0

..

q :

.

.㌧o

o

o

‑ & ･･.

cq

.

& o ･･･ o

‑AIAJ{+

B

とな る.行列AL,BLの右上 と左下が oとなっている. これ は,図3において負周波数成分 を 0とおいていることに対応 している.GWD は, この部分 まで計算 しているために,点数 が少 なければCWDの方が高速 に計算できる.

WDは信号の局所周波数を解析す る有効 な手法であ り, ウイン ドウ点数が短いほどその特 徴が得 られ る.ただ し,あま り短 いと解析 で きな くある程度 の点数 はいる.体表面心電位 の 周波数解析にWDを利用 した報告(8)では, ウィ′ソ ドゥ点数16であったが,一般 に32,64点が 最低最適 であ ると考 えられ るため, この点数での計算量低下 はWDにおいて特 に有効 であ

る.

5

計 算 例

FT,GWD,CWDの各 計 算 時 間 と各 比 を表 1に示 す.各 値 は,80386SX(16MHz)+

80387を使用 し,各点数 を1万 回実行 した ときの時間 と,各点数の計算以外 の処理を2万回

表1 各処理の計算時間およびその比

点数 各種処理法 (単位ms) 計 算 時 間 比 FT GWD CWD GWD:FT CWD:FT CWD:GWD 16 4.7 15.1 14.6 3.21 3.ll 0.97 32 10.0 33.8 26.6 3.38 2.66 0.79 64 23.4 76.8 65.3 3.28 2.79 0.85

実行 した ときの時間を測定 し,それぞれ1回分 の時間 にしたその差 を1回の計算時間 とした.

これ よ り,GWD :FTは,GWDが3回のFTを必要 とす るこ とですべ て3倍 強 にな っ ている.CWD:GWDは32,64点で80%は どの実行時間で計算が可能であるが,128点以上 ではGWDの方が早 く畳込みに よる時間がかか ることがわかる.時変信号の解析 に必要 とさ れ る最適点数 (32,64点)で は,CWDがGWD よ り優 って い るため,本手 法 がWDの解 析手法に有効 だ といえる.

6 お わ り に

本稿では,WD の一般的計算手法 を説明 した後 に,WDがFSの畳込みによ り計算 で きる こと確認す るとともに,それを利用 した畳込み による計算が,時変信号の解析 に必要 とされ る最適点数 (32,64点)では,80%ほ ど従来の手法 よ り早 く計算で きることを実測結果 よ り 示 した.本手法で あるWDの高速計算で も,FTの2.7倍 の時間が かか るため,本手法 の‑

‑ ドウエアによる高速化が今後の課題である.

7

参 考 文 献

1. ClaasenT.A.C.M.andMechlenbmkerW.F.G.;"WignerDistribution‑ToolforTime･ FrequencySignalAnalysis･PART 1.2.3",PhilipsJ.Recリ35,pp.217‑250,pp.276‑300,pp.

372‑389(1980).

2. BoashashB.:"NoteontheUseoftheWignerDistributionforTime‑FrequencySignal Analysis",IEEETrams.Acoust.,Speech&SignalProcess.,ASSP‑36,9,pp.1518‑1521(1988). 3.鈴木宏,小林史典 :"畳込みによるWigner分布の計算",第10回計測自動制御学会九州支部学

術講演会,140,pp.93‑94(1991).

4. P.A.Ramamoorthy:"AutoregressiveModelingoftheWignerSpectrum",ICASSP87,35. 9(1987).

5. H.Garudadri:"OnComputingtheSmoothedWignerDistribution",ICASSP87,35.12(1987). 6. Lee:"ANewAlgorithm toCompute也eDiscreteCosineTransform",IEEETrams.Acoust.,

Speech&SignalProcess.,ASSP‑32,6,pp.124311245(1984).

7.高橋秀俊 :"高速7‑1)‑交換(FFT)について",情報処理,Vol.14,No.8pp.616‑622(1973). 8.白井支朗 ら :̀̀ゥィグナ‑分布による体表面心電位図の局所空間周波数解析",第2回デ ィジタ

ル信号処理シンポジウム,B‑4‑1,pp.2531258(1987).

付 録

スペク トル数列Xl(A)をXl(k)‑Ah+jBAと置 いた とき,いまこのXl(k)とこの複素共役 ('で表す)Xl'(k)の畳込み W (k)は,

rLn‑il

W(k)‑pS^.Xl(♪)Xt(k‑♪) [Al] で定義 されている. ここで,複素数

X

と虚郡 はそれぞれ

72 鈴 木 宏

A(‑I

gtlw(k)]‑pS.(mlXl(b)]gtlXl'(k‑♪)]

+

b

]

〟 ‑1

=^{昌 (ApA(A‑P},+^{BpB糾 ,)}

琵[W(k)]‑0

とな り,式(8)と等 しくなる.

つ ぎに,時間軸上の掛算 として式[Al】を考 える.Xl(k)の逆FTxl(n)は, xl(n)

= 去

である.xl(n)とxl'(n)との掛算xl(n)xr(n)のFT(X(k))は,

111I

X(k)‑pS.Xl(b)Xl'(♪‑ k)

lA5]

[ A6 ]

x 2(k)で式

【 A6

〟 ‑I

X'(k)‑p^{S.Xl(♪)x2(k‑b)‑ W(A) lA7]}

とな り,式【Al]に等 しくなる. このx2(k)の逆FTした苑(n)は,xl(n)の複素共役で しか も中点を対称 として左右の値を交換 したものである.

式(8)の畳込みを求めるには,FTした ものを単 に掛 ければよいのではな く,信号

X

した もの,すなわち図3のC部の計算 と等 しくなる.