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

10 views (last 30 days)
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

Accepted Answer

Naoya
Naoya on 14 Jul 2017
まず、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 程度と解釈できます。
  3 Comments
Ryosuke Takahashi
Ryosuke Takahashi on 25 Jul 2017
ご回答ありがとうございました。 今後ともよろしくお願いいたします。

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!