このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

poctave

オクターブ スペクトルの生成

構文

p = poctave(x,fs)

p = poctave(xt)

p = poctave(pxx,fs,f)

p = poctave(___,type)

p = poctave(___,Name,Value)

[p,cf] = poctave(___)

[p,cf,t] = poctave(___)

poctave(___)

説明

p = poctave(x,fs) は、fs のレートでサンプリングされた信号 x のオクターブ スペクトルを返します。オクターブ スペクトルは、ANSI S1.11 規格[2]で定義されているオクターブ帯域上の平均パワーです。x が行列の場合、この関数は各列のオクターブ スペクトルを個別に推定し、結果を対応する p の列に返します。

p = poctave(xt) は、MATLAB^® timetable xt に保存された信号のオクターブ スペクトルを返します。

p = poctave(pxx,fs,f) では、パワー スペクトル密度、pxx を 1/b オクターブ パワー スペクトルに変換することによってオクターブ平滑化を実行します (b はオクターブ帯域内のサブバンドの数)。f の周波数は、pxx の PSD 推定値に相当します。

p = poctave(___,type) は、関数で実行されるスペクトル解析の種類を指定します。type は、'power' または 'spectrogram' に指定します。

p = poctave(___,Name,Value) は、名前と値の引数を使用して、上記の任意の構文に追加オプションを指定します。

[p,cf] = poctave(___) は、オクターブ スペクトルが計算されるオクターブ帯域の中心周波数も返します。

[p,cf,t] = poctave(___) は、type が 'spectrogram' の場合、パワー スペクトルの推定の計算に使用されるセグメントの中央時間に対応する時点のベクトル t を追加で返します。

出力引数なしで poctave(___) を使用すると、現在の Figure にオクターブ スペクトルまたはスペクトログラムがプロットされます。type を 'spectrogram' に指定した場合、この関数は単一チャネルの入力でのみサポートされます。

例

ホワイトおよびピンク ノイズのオクターブ スペクトル

ホワイト ガウス ノイズを ${10}^5$ サンプル生成します。零点と極がすべて正の x 軸上にあるフィルターでホワイト ノイズをフィルター処理することによって、擬似ピンク ノイズの信号を作成します。零点と極を可視化します。

N = 1e5;
wn = randn(N,1);

z = [0.982231570015379 0.832656605953720 0.107980893771348]';
p = [0.995168968915815 0.943841773712820 0.555945259371364]';

[b,a] = zp2tf(z,p,1);
pn = filter(b,a,wn);

zplane(z,p)

ホワイトおよびピンク ノイズで構成される 2 チャネル信号を作成します。オクターブ スペクトルを計算します。サンプル レートは 44.1 kHz だと仮定します。周波数帯を 30 Hz からナイキスト周波数の範囲に設定します。

sg = [wn pn];

fs = 44100;

poctave(sg,fs,'FrequencyLimits',[30 fs/2])
legend('White noise','Pink noise','Location','SouthEast')

Warning: Negative data ignored

ホワイト ノイズには、周波数と共に増加するオクターブ スペクトルがあります。ピンク ノイズのオクターブ スペクトルは、周波数範囲全体をとおしてほぼ一定です。信号のオクターブ スペクトルは、人間の耳が信号をどのように認識するかを示します。

ホワイトおよびピンク ノイズのオクターブ平滑化

44.1 kHz でサンプリングされたホワイト ガウス ノイズを ${10}^5$ サンプル生成します。零点と極がすべて正の x 軸上にあるフィルターでホワイト ノイズをフィルター処理することによって、ピンク ノイズの信号を作成します。

N = 1e5;
fs = 44.1e3;
wn = randn(N,1);

z = [0.982231570015379 0.832656605953720 0.107980893771348]';
p = [0.995168968915815 0.943841773712820 0.555945259371364]';
[b,a] = zp2tf(z,p,1);

pn = filter(b,a,wn);

両方の信号についてパワー スペクトル密度のウェルチ推定を計算します。信号を 2048 サンプルのセグメントに分割し、隣接するセグメント間に 50% のオーバーラップを指定し、各セグメントをハミング ウィンドウでウィンドウ処理し、4096 DFT ポイントを使用します。

[pxx,f] = pwelch([wn pn],hamming(2048),1024,4096,fs);

200 Hz からナイキスト周波数までの範囲の周波数帯上のスペクトル密度を表示します。周波数軸に対数スケールを使用します。

pwelch([wn pn],hamming(2048),1024,4096,fs)
ax = gca;
ax.XScale = 'log';
xlim([200 fs/2]/1000)
legend('White','Pink')

信号のオクターブ スペクトルを計算し、表示します。前のプロットと同じ周波数範囲を使用します。1 オクターブあたり 6 帯域を指定し、8 次フィルターを使用してスペクトルを計算します。

poctave(pxx,fs,f,'BandsPerOctave',6,'FilterOrder',8,'FrequencyLimits',[200 fs/2],'psd')
legend('White','Pink')

音声信号のオクターブ スペクトログラム

電子歯ブラシの音声録音を MATLAB® に読み取ります。歯ブラシは約 1.75 秒で電源が入り、約 2 秒間持続します。

[y,fs] = audioread('toothbrush.m4a');

音声信号のオクターブ スペクトログラムを計算します。1 オクターブあたり 48 帯域、82% のオーバーラップを指定します。全体の周波数範囲を 100 Hz ～ fs/2 Hz に制限し、C 重み付けを使用します。

poctave(y,fs,'spectrogram','BandsPerOctave',48,'OverlapPercent',82,'FrequencyLimits',[100 fs/2],'Weighting','C')

オクターブ スペクトルの重み付け

44.1 kHz でサンプリングされたホワイト ガウス ノイズを ${10}^5$ サンプル生成します。零点と極がすべて正の x 軸上にあるフィルターでホワイト ノイズをフィルター処理することによって、ピンク ノイズの信号を作成します。

N = 1e5;
fs = 44.1e3;
wn = randn(N,1);

z = [0.982231570015379 0.832656605953720 0.107980893771348]';
p = [0.995168968915815 0.943841773712820 0.555945259371364]';
[b,a] = zp2tf(z,p,1);

pn = filter(b,a,wn);

信号のオクターブ スペクトルを計算します。1 オクターブあたり 3 つの帯域を指定して、全体の周波数範囲を 200 Hz ～ 20 kHz に制限します。後で利用するために cell 配列に名前と値のペアを保存します。スペクトルを表示します。

flims = [200 20e3];
bpo = 3;
opts = {'FrequencyLimits',flims,'BandsPerOctave',bpo};

poctave(pn,fs,opts{:});

同じ設定で、ただし C 重み付けを使用して信号のオクターブ スペクトルを計算します。C 重み付けのスペクトルは、6 kHz を超える周波数では低下します。

hold on
poctave(pn,fs,opts{:},'Weighting','C')

今度は A 重み付けを使用して、再度オクターブ スペクトルを計算します。A 重み付けのスペクトルは、約 3 kHz でピークになり、6 kHz を超える場合と周波数帯の下端では低下します。

poctave(pn,fs,opts{:},'Weighting','A')
hold off
legend('Pink noise','C-weighted','A-weighted','Location','SouthWest')

入力引数

x — 入力信号

ベクトル | 行列

ベクトルまたは行列として指定される入力信号。x がベクトルの場合、poctave では単一チャネルとして取り扱います。x が行列の場合、poctave は各列のオクターブ スペクトルまたはスペクトログラムを個別に計算し、結果を対応する p の列に返します。type を 'spectrogram' に設定した場合、この関数は p の 3 番目の次元に沿ってスペクトログラムを連結します。

例: sin(2*pi*(0:127)/16)+randn(1,128)/100 は、ノイズを含む正弦波を指定します。

例: [2 1].*sin(2*pi*(0:127)'./[16 64]) は、2 チャネルの正弦波を指定します。

データ型: single | double

fs — サンプル レート

正のスカラー

サンプル レート。正のスカラーとしてヘルツ単位で指定します。サンプル レートは 7 Hz より低くできません。

xt — 入力 timetable

timetable

入力 timetable。xt に増加する有限の等間隔の行時間を含めなければなりません。xt がマルチチャネル信号を表す場合は、行列を含む単一変数か、ベクトルで構成される複数の変数のどちらかをもたなければなりません。

timetable が欠損している場合や時間点が重複している場合、欠損または重複する時間および非等間隔の時間をもつ timetable の整理のヒントを使用して修正できます。

例: timetable(seconds(0:4)',randn(5,1)) は 1 Hz で 4 秒間サンプリングされたランダム過程を指定します。

pxx — パワー スペクトル密度

ベクトル | 行列

パワー スペクトル密度 (PSD)。実数の非負の要素を含むベクトルまたは行列として指定します。パワー スペクトル密度は、デシベル単位ではなく、線形単位で表さなければなりません。デシベル値をパワー値に変換するには、db2pow を使用します。type が 'spectrogram' の場合、pxx の各列は特定の時間枠またはサンプルの PSD と見なされます。

例: [pxx,f] = periodogram(cos(pi./[4;2]*(0:159))'+randn(160,2)) は、2π Hz でサンプリングされたノイズの多い 2 チャネル正弦波のピリオドグラム PSD 推定値とそれについて計算される周波数を指定します。

f — PSD 周波数

ベクトル

PSD 周波数。ベクトルとして指定します。f は有限で、厳密に増加し、線形スケールにおいて等間隔でなければなりません。

例: [pxx,f] = periodogram(cos(pi./[4;2]*(0:159))'+randn(160,2)) は、2π Hz でサンプリングされたノイズの多い 2 チャネル正弦波のピリオドグラム PSD 推定値とそれについて計算される周波数を指定します。

type — 計算するスペクトルのタイプ

'power' (既定値) | 'spectrogram'

計算するスペクトルのタイプ。'power' または 'spectrogram' として指定します。

• 'power' — 入力のオクターブ パワー スペクトルを計算します。

• 'spectrogram' — 入力のオクターブ スペクトログラムを計算します。この関数は入力をセグメントに分割し、各セグメントの短時間オクターブ パワー スペクトルを返します。

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、コンマを使用して名前と値をそれぞれ区切り、Name を引用符で囲みます。

例: 'Weighting','A','FilterOrder',8 は、A 重み付けと 8 次フィルターを使用してオクターブ スペクトルを計算します。

BandsPerOctave — オクターブ帯域におけるサブバンドの数

1 (既定値) | 3/2 | 2 | 3 | 6 | 12 | 24 | 48 | 96

オクターブ帯域におけるサブバンドの数。1、3/2、2、3、6、12、24、48、または 96 として指定します。このパラメーターは、1/N オクターブ帯域の幅を規定します。このような周波数帯において、上端の周波数は下端周波数の 2^1/b 倍です。ここで、b はサブバンドの数です。

データ型: single | double

FilterOrder — バンドパス フィルターの次数

6 (既定値) | 正の偶数の整数

バンドパス フィルターの次数。正の偶数の整数として指定します。

データ型: single | double

FrequencyLimits — 周波数帯

[max(3,3*fs/48e3) fs/2] (既定値) | 2 要素ベクトル

周波数帯。ヘルツ単位で表される 2 要素の増加ベクトルとして指定します。低い方のベクトル値は、3 Hz 以上でなければなりません。高い方のベクトル値は、ナイキスト周波数以下でなければなりません。ベクトルにオクターブ中心が含まれていない場合、poctave は指定された制限の範囲外の中心周波数を返すことがあります。安定したフィルター設計にするには、サンプル レートが 48 kHz を超えた場合に達成可能な実際の最小周波数範囲を 3*fs/48e3 に増やします。この引数が指定されていない場合、poctave は [max(3,3*fs/48e3) fs/2] の間隔を使用します。

データ型: single | double

Weighting — 周波数の重み付け

'none' (既定値) | 'A' | 'C' | ベクトル | 行列 | 1 行 2 列の cell 配列 | digitalFilter オブジェクト

周波数の重み付け。次のいずれかとして指定します。

• 'none' — poctave は入力に対してどのような周波数の重み付けも行いません。

• 'A' — poctave は入力に対して A 重み付けを実行します。ANSI S1.42 規格は、A 重み付け曲線を定義しています。IEC 61672-1 規格は、A 重み付けフィルターの最小および最大減衰範囲を定義しています。ANSI S1.42.2001 規格は、アナログの極および零点を指定することによって重み付け曲線を定義しています。

• 'C' — poctave は入力に対して C 重み付けを実行します。ANSI S1.42 規格は、C 重み付け曲線を定義しています。IEC 61672-1 規格は、C 重み付けフィルターの最小および最大減衰範囲を定義しています。ANSI S1.42.2001 規格は、アナログの極および零点を指定することによって重み付け曲線を定義しています。

• ベクトル — poctave では、入力を有限インパルス応答 (FIR) フィルターを指定する係数のベクトルとして処理します。

• 行列 — poctave では、入力を無限インパルス応答 (IIR) フィルターを指定する 2 次セクション係数の行列として処理します。行列には、少なくとも 2 つの行と厳密に 6 つの列がなくてはなりません。

• 1 行 2 列 cell 配列 — poctave では、入力を IIR フィルターの伝達関数を指定するその次数の分子と分母の係数として処理します。

• digitalFilter オブジェクト — poctave は、入力を、designfilt を使用して設計されたフィルターとして処理します。

この引数は、入力が信号の場合にのみサポートされます。オクターブ平滑化は、周波数の重み付けをサポートしていません。

例: 'Weighting',fir1(30,0.5) は、0.5π ラジアン/サンプルの正規化されたカットオフ周波数をもつ 30 次の FIR フィルターを指定します。

例: 'Weighting',[2 4 2 6 0 2;3 3 0 6 0 0] は、正規化された 3 dB の周波数 0.5π ラジアン/サンプルをもつ 3 次のバタワース フィルターを指定します。

例: 'Weighting',{[1 3 3 1]/6 [3 0 1]/3} は、正規化された 3 dB の周波数 0.5π ラジアン/サンプルをもつ 3 次のバタワース フィルターを指定します。

例: 'Weighting',designfilt('lowpassiir','FilterOrder',3,'HalfPowerFrequency',0.5) は、正規化された 3 dB の周波数 0.5π ラジアン/サンプルをもつ 3 次のバタワース フィルターを指定します。

データ型: single | double | char | string | cell

MinThreshold — 非ゼロ値の下限

-Inf (既定値) | 実数スカラー

非ゼロ値の下限。実数スカラーとして指定します。この関数は、10 log₁₀(p) ≤ 'MinThreshold' がゼロとなるように p の要素を設定します。'MinThreshold' をデシベル単位で指定します。

データ型: single | double

WindowLength — データ セグメントの長さ

非負の整数

データ セグメントの長さ。非負の整数として指定します。'WindowLength' は入力信号の長さ以下にしなければなりません。指定しない場合、データ セグメントの長さは入力信号のサイズに基づいて計算されます。この入力は、type が 'spectrogram' の場合にのみ有効です。

データ型: single | double

OverlapPercent — 隣り合ったセグメント間のオーバーラップ率

[0, 100) の範囲にある実数スカラー

隣り合ったセグメント間のオーバーラップ率。[0, 100) の範囲にある実数スカラーとして指定します。指定しない場合、'OverlapPercent' はゼロとなります。この入力は、type が 'spectrogram' の場合にのみ有効です。

データ型: single | double

出力引数

p — オクターブ スペクトルまたはスペクトログラム

ベクトル | 行列 | 3 次元配列

オクターブ スペクトルまたはスペクトログラム。ベクトル、行列、または 3 次元配列として返されます。3 番目の次元が存在する場合は、入力チャネルに対応します。

cf — 中心周波数

ベクトル

中心周波数。ベクトルとして返されます。cf には、poctave がオクターブ スペクトルを推定したオクターブ帯域の中心周波数のリストが含まれます。cf の単位は Hz です。

t — 中央時間

ベクトル

中央時間。ベクトルとして返されます。入力が PSD の場合、t は pxx の列に対応するサンプル インデックスを表します。この引数は、type が 'spectrogram' の場合にのみ適用されます。

アルゴリズム

"オクターブ解析" は、人間の耳が音を知覚するのに似たプロセスで、広い周波数範囲の音や振動レベルを特定するために使用されます。信号スペクトルはオクターブまたは 1/N オクターブ帯域に分割されます。各帯域の周波数制限は低い周波数制限の 2 倍となっており、帯域幅は周波数が高くなるほど広くなります。

オクターブ フィルターの使用

オクターブ解析を実行するために、poctave 関数は並列バンドパス フィルターのフィルター バンクを作成します。各デジタル バンドパス フィルターは、対応するバタワース ローパス アナログ フィルター [3] にマッピングされています。アナログ フィルターは、バンドパス型の bilinear 変換を使用してデジタル バンドパス フィルターにマッピングされ、その結果は 4 次セクションのカスケードとして返されます。

ANSI S1.11-2004 規格[2]では、オクターブ帯域の中心周波数を次のように定義しています。

$$f_c=\{\begin{array}{ll}{f}_r\times G^{(k-30)/b},\hfill & b\text{ is odd}\hfill \\ f_r\times G^{(2k-59)/2b},\hfill & b\text{ is even}\hfill \end{array}$$

ここで、

• f_r は、基準周波数 1000 Hz です。

• G は、オクターブ比 10^0.3 です。

• b は、BandsPerOctave で指定されるオクターブあたりのバンド数です。

• k は任意の整数であり、バンド番号を表します。

中心周波数の定義は、b が奇数であるか (1 オクターブまたは 1/3 オクターブの帯域幅の場合)、または偶数であるか (1/2 オクターブまたは 1/6 オクターブの帯域幅の場合) によって変化します。これにより、オクターブ帯域全体の帯域エッジが必ずすべての比帯域の帯域エッジとしてとどまります。

各オクターブ帯域の低域側と広域側のエッジ周波数は、次のように求められます。

$$f_l=f_c\cdot (G^{-1/2b})$$

$$f_u=f_c\cdot (G^{1/2b})$$

オクターブ フィルターの設計と実装の詳細については、Digital Filter Design (Audio Toolbox)を参照してください。

オクターブ平滑化の使用

関数 poctave は、帯域内の信号のパワー スペクトル密度 (PSD) を長方形近似法で積分することで、各オクターブ帯域の平均パワーを計算します。オクターブ帯域の平均パワーは、帯域の中心周波数における信号レベルを表しています。

• 帯域エッジがビン内にある場合、この関数は、帯域が占める周波数ビンの割合に対応するパワーの割合のみを帯域に割り当てます。たとえば、次の図はオクターブ帯域のエッジが、オレンジと青の破線の長方形で表された 2 つの異なる周波数ビン内にあることを示しています。影付きの領域内のパワーは、指定されたオクターブ帯域について計算されます。

  

• 帯域エッジが 0 またはナイキスト周波数 f_Nyquist にある場合、この関数は、帯域が占める周波数ビンの割合に対応するパワーの割合の 2 倍を帯域に割り当てます。この重複は、[–w/2, 0] および [f_Nyquist, f_Nyquist + w/2] の範囲に存在するビンのパワーの半分を占め、ここで w はビンの幅です。たとえば、次の図は右端がナイキスト周波数にあるオクターブ帯域を示しています。影付きの領域内のパワーは、指定されたオクターブ帯域について計算されます。

  

参照

[1] Smith, Julius Orion, III. "Example: Synthesis of 1/F Noise (Pink Noise)." In Spectral Audio Signal Processing. https://ccrma.stanford.edu/~jos/sasp/.

[2] Specification for Octave-Band and Fractional-Octave-Band Analog and Digital Filters. ANSI Standard S1.11-2004. Melville, NY: Acoustical Society of America, 2004.

[3] Orfanidis, Sophocles J. Introduction to Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 2010.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• コード生成時に入力の重み付けを可変サイズの行列で指定した場合、実行時にベクトルに変換してはいけません。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

使用上の注意および制限:

• 関数 poctave により、不安定で gpuArray 入力に対応していない低周波数オクターブ フィルターが生成される場合があります。オプションの名前と値の引数 FrequencyLimits を使用して、低域帯域周波数を増やすか、入力信号をダウンサンプリングしてください。

詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2018a で導入

R2022b: gpuArray オブジェクトの使用

関数 poctave は gpuArray オブジェクトをサポートしています。この機能を使用するには、Parallel Computing Toolbox™ が必要です。詳細については、GPU 配列を参照してください。

参考

pspectrum