scal2frq

音階から周波数

構文

frq = scal2frq(A,wname,delta)

frq = scal2frq(A,wname)

説明

frq = scal2frq(A,wname,delta) は、A で与えられた音階および wname とサンプリング周期 delta で指定されたウェーブレットに対応する疑似周波数を返します。出力 frq は実数値で、次元は A と同じになります。

frq = scal2frq(A,wname) は frq = scal2frq(A,wname,1) と等価です。

例

すべて折りたたむ

音階と疑似周波数

ライブスクリプトを開く

この例では、音階を 2 倍にしたときに疑似周波数がどのように変化するかを示します。

5 オクターブにわたってオクターブあたりの音が 10 となる音階のベクトルを作成します。

vpo = 10;
no = 5;
a0 = 2^(1/vpo);
ind = 0:vpo*no;
sc = a0.^ind;

音階の範囲が 5 オクターブをカバーしていることを確認します。

log2(max(sc)/min(sc))

ans = 5.0000

音階をプロットする場合、データカーソルを使用して、インデックス $n + 10$ の音階がインデックス $n$ の音階の 2 倍になることを確認できます。各オクターブをマークするように y 軸目盛りを設定します。

plot(ind,sc)
title('Scales')
xlabel('Index')
ylabel('Scale')
grid on
set(gca,'YTick',2.^(0:5))

Figure contains an axes object. The axes object with title Scales, xlabel Index, ylabel Scale contains an object of type line.

音階を実数値 Morlet ウェーブレットの疑似周波数に変換します。まず、サンプリング周期を 1 と仮定します。

pf = scal2frq(sc,"morl");
T = [sc(:) pf(:)];
T = array2table(T,'VariableNames',{'Scale','Pseudo-Frequency'});
disp(T)

    Scale     Pseudo-Frequency
    ______    ________________

         1          0.8125    
    1.0718         0.75809    
    1.1487         0.70732    
    1.2311         0.65996    
    1.3195         0.61576    
    1.4142         0.57452    
    1.5157         0.53605    
    1.6245         0.50015    
    1.7411         0.46666    
    1.8661         0.43541    
         2         0.40625    
    2.1435         0.37904    
    2.2974         0.35366    
    2.4623         0.32998    
     2.639         0.30788    
    2.8284         0.28726    
    3.0314         0.26803    
     3.249         0.25008    
    3.4822         0.23333    
    3.7321          0.2177    
         4         0.20313    
    4.2871         0.18952    
    4.5948         0.17683    
    4.9246         0.16499    
     5.278         0.15394    
    5.6569         0.14363    
    6.0629         0.13401    
     6.498         0.12504    
    6.9644         0.11666    
    7.4643         0.10885    
         8         0.10156    
    8.5742        0.094761    
    9.1896        0.088415    
    9.8492        0.082494    
    10.556         0.07697    
    11.314        0.071816    
    12.126        0.067006    
    12.996        0.062519    
    13.929        0.058332    
    14.929        0.054426    
        16        0.050781    
    17.148        0.047381    
    18.379        0.044208    
    19.698        0.041247    
    21.112        0.038485    
    22.627        0.035908    
    24.251        0.033503    
    25.992         0.03126    
    27.858        0.029166    
    29.857        0.027213    
        32        0.025391

100 Hz でデータがサンプリングされると仮定します。音階、対応する疑似周波数、および周期を含む table を作成します。オクターブあたり 10 音なので、表の 10 行目ごとに表示します。音階が 2 倍になるたびに疑似周波数が半分になっていることを確認します。

Fs = 100;
DT = 1/Fs;
pf = scal2frq(sc,"morl",DT);
T = [sc(:)/Fs pf(:) 1./pf(:)];
T = array2table(T,'VariableNames',{'Scale','Pseudo-Frequency','Period'});
T(1:vpo:end,:)

ans=6×3 table
    Scale    Pseudo-Frequency     Period 
    _____    ________________    ________

    0.01           81.25         0.012308
    0.02          40.625         0.024615
    0.04          20.313         0.049231
    0.08          10.156         0.098462
    0.16          5.0781          0.19692
    0.32          2.5391          0.39385

scal2frq に $Δ t = \frac{1}{F s}$ の係数があることに注意してください。これは、音階から周波数への適切な変換を実現するために必要です。 $Δ t$ は生の音階を適切に調整するために必要です。たとえば、次のようにするとします。

f = scal2frq(1,'morl',0.01);

これは、マザー Morlet ウェーブレットを 0.01 倍に拡大した場合、ウェーブレットの中心周波数はどうなるかということです。言い換えれば、 $ψ (t)$ の代わりに $ψ (t / 0.01)$ を見た場合、中心周波数にどのような影響があるかということです。 $Δ t$ によって、音階の正しい調整係数が提供されます。

調整したサイズに音階をまず変換してから scal2frq を指定せずに $Δ t$ を使用することでも、同じ結果が得られます。

scadjusted = sc.*0.01;
pf2 = scal2frq(scadjusted,'morl');
max(pf-pf2)

ans = 0

周波数を使用した等高線図での CWT のプロット

ライブスクリプトを開く

この例では、Hz 単位の近似周波数を使用して CWT の等高線図を作成する方法を示します。

加法性ノイズのサポートが互いに素である 2 つの正弦波で構成される信号を作成します。1 kHz で信号がサンプリングされると仮定します。

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = 1.5*cos(2*pi*100*t).*(t<0.25)+1.5*cos(2*pi*50*t).*(t>0.5 & t<=0.75);
x = x+0.05*randn(size(t));

入力信号の CWT を取得し、結果をプロットします。

[cfs,f] = cwt(x,Fs);
contour(t,f,abs(cfs).^2); 
axis tight;
grid on;
xlabel('Time');
ylabel('Approximate Frequency (Hz)');
title('CWT with Time vs Frequency');

Figure contains an axes object. The axes object with title CWT with Time vs Frequency, xlabel Time, ylabel Approximate Frequency (Hz) contains an object of type contour.

入力引数

すべて折りたたむ

`A` — 音階
正の実数値のベクトル

音階。正の実数値のベクトルとして指定します。

`wname` — ウェーブレット
文字ベクトル | string スカラー

ウェーブレット。文字ベクトルまたは string スカラーとして指定します。詳細については、wavefun を参照してください。

`delta` — サンプリング周期
`1` (既定値) | 正の実数値スカラー

サンプリング周期。実数値のスカラーとして指定します。

例: pf = scal2frq([1:5],"db4",0.01)

詳細

すべて折りたたむ

疑似周波数

音階と周波数の関係については、おおよその答えしかありません。

ウェーブレット解析において、音階を周波数に関連付ける方法は、ウェーブレットの中心周波数 F_c を決定し、次の関係を使用することです。

$F_{a} = \frac{F_{c}}{a}$

ここで、

a は音階。
F_c はウェーブレットの中心周波数 (Hz)。
F_a は音階 a に対応する疑似周波数 (Hz)。

この考え方は、与えられたウェーブレットに周波数 F_c の純粋に周期的な信号を関連付けるというものです。ウェーブレット係数のフーリエ変換を最大化する周波数は F_c です。関数 centfrq は、指定されたウェーブレットの中心周波数を計算します。上述の関係から、音階は疑似周波数に反比例することがわかります。たとえば、音階が上がると、ウェーブレットはさらに広がり、疑似周波数は低くなります。

中心周波数とウェーブレットの対応の例を次の Figure に示します。

実数ウェーブレットと複素数ウェーブレットの中心周波数

このように、中心周波数ベースの近似 (赤) は、ウェーブレットの主要な振動 (青) を捉えています。中心周波数は、ウェーブレットの支配的な周波数を便利かつ簡単に特徴付けるものです。

参照

[1] Abry, P. Ondelettes et turbulence. Multirésolutions, algorithmes de décomposition, invariance d'échelles et signaux de pression. Diderot, Editeurs des sciences et des arts, Paris, 1997.