wavelet波形の周波数分解能について

wavelet波形の周波数分解のについて質問があります。

scal2freq(スケール，マザーウェブレット，1/周波数)でマザーウェブレットの帯域幅が算出されると思います。 (間違っていたら教えてください)

例えば，サンプリング周波数 2000Hz で、マザーウェーブレット cmor1-1 とした場合、

スケールファクタが 1000 の時、下記のように 2Hz となります。 これは，中心周波数10Hzの場合，9〜11Hzまでは同様の結果になると考えております。

ここからが質問ですが，各周波数における帯域幅を算出するにはどのようにしたらいいのでしょうか？ 下記のようなsin波のモデル波形を作成した際，20Hzと100Hzにおける周波数分解能を知りたいです。

申し訳ありませんが，ご教示のほど宜しくお願い致します。

%%

fs=1000; t = 0:1/fs:2-1/fs;

S1 = sin(2*pi*20*t)+sin(2*pi*100*t);

time = (1:fs*2)/fs;

figure(1);

subplot(2,1,1);

plot(time,S1);

subplot(2,1,2);

wname = 'morl';

fc = centfrq(wname); % 中心周波数

fa = 1:200; % 擬似周波数(Hz)

scal2frq(200,wname,1/fs)%周波数分解能

sf = fc./(fa.*1/fs); % スケールファクタ

[CWTcoeffs,frq] = cwt(S1,sf,wname,1/fs);

abs_CWT = abs(CWTcoeffs);

imagesc(time,fa,abs_CWT);

colormap(jet)

axis xy

axis([0,inf, -inf, 200])

title('Scalogram')

ylabel('Hz')

colorbar

Answer 1

[Naoya]

まず、scal2frq は、帯域幅を求めるものではなく、疑似周波数に対する、中心周波数を求める機能となります。 残念ながら、疑似周波数に対しての帯域幅を確認する直接的な方法は提供されていません。

お手数をお掛けすることになりますが、下記のようなスクリプト例で帯域幅を確認することができます (原理的には、インパルス信号に対して、疑似周波数に相当するスケールファクタを設定した CWT を適用します。)

下記例は、疑似周波数が 20Hz の場合となります。

fs=2000;

t = 0:1/fs:2-1/fs;

wname = 'cmor1-1';

% すべて0の値を持つ信号の中心にインパルスを付加

sig = zeros(1,length(t));

sig(length(t)/2) = 1;

sf = 100; % scaleファクタ = 100 => Fc = 20Hz

filterd = cwt(sig, sf, wname);

Y = fft(filterd);

F = (0:3999)/4000*fs;

plot(F, abs(Y))

% 中心周波数付近を拡大

xlim([fs-2*scal2frq(sf,'cmor1-1',1/fs) fs])

中心周波数は 1980Hz と現れますが、これは20Hz のミラー成分に相当します. どのあたりまでの振幅の大きさを帯域幅と捉えるかに寄りますが、 例えば、中心周波数に対して、振幅が半分までを帯域幅と捉える場合は、+/-5Hz 程度と解釈できます。

[Ryosuke Takahashi]

返信が遅れてしまい申し訳ありませんでした。

とても助かりました。 今回は20Hzの場合でしたが，10Hzで確認したい場合はscal2frq を使用して擬似周波数を算出し，sfに算出した値を入力する方法であっていますでしょうか？

[Naoya]

はい、その理解でよろしいかとおもいます。 下記FAQ

http://www.mathworks.com/matlabcentral/answers/103712-1

Fa = Fc ./ (scalef * delta)

として 中心周波数(Fc), 疑似周波数(Fa), スケールファクタ(scalef), サンプリング時間(delta)の関係式がありますが 疑似周波数に対してのスケールファクタは

Scalef = Fc / (Fa * delta)

のようにきめられますので、

Fc = centfrq('cmor1-1'); Fa = 10; delta = 1/2000; Scalef = Fc / (Fa * delta)

Scalef =

   200

のようにスケールファクタを求めることができます。

[Ryosuke Takahashi]

ご回答ありがとうございました。 今後ともよろしくお願いいたします。