MATLABで学ぶ信号処理・音響工学 音場制御の評価指標について1

はじめに

これまでに音場制御のプログラミング方法について説明しました。下記記事内容を理解していることを前提に本記事を書き進めていますので、まだ読んでいない方は、まず最初に下記の記事を読んで頂けると幸いです。

MATLABで学ぶ信号処理・音響工学 音場制御のプログラミングについて1
本記事ではフィードフォワードの音場制御方法の理論について、数式とMATLABのコードを交えながら説明いたします。理論としてはすごくシンプルなので、高校生でも理解できると思います。
MATLABで学ぶ信号処理・音響工学 音場制御のプログラミングについて2 (最小二乗法を用いた最適化)
最小二乗法を用いた音場最適化制御方法について、理論式やMATLABプログラムを示しながら説明します。ここまで簡単に詳しく説明してくれている文献はほかにないので、ぜひ読んでみてください。
MATLABで学ぶ信号処理・音響工学 音場制御のプログラミングについて3 (効果的な目標音場の設定)
最小二乗法を用いた音場最適化制御方法について、理論式やMATLABプログラムを示しながら説明します。この音場最適化手法では、目標音場の設定が重要となります。この目標音場の設定方法について詳しく説明いたします。ここまで簡単に詳しく説明してくれている文献はほかにないので、ぜひ読んでみてください。

これまでの記事で音場制御方法について説明しましたが、”音を伝搬させるブライトゾーンと、音を伝搬させないダークゾーンの2つの音場をうまく制御できているかの指標がないと、音場制御方法を評価することが難しい”ということがわかりました。

そこで、ここでは音場制御を評価するための指標について説明します。

スポンサーリンク

音場制御を評価するための指標

音場制御を評価するための指標は大きく下記の3つです。

    1. Array Effort
    2. Acoustic Contrast
    3. ブライトゾーン内の音圧のばらつき
1. Array Effort

音場制御をするには、スピーカーから音を出します。このとき、各スピーカーに入力する入力電圧の総和(入力パワー)は小さい方が、より効率的な制御であると言えます。この効率性を示す指標がArray Effortであり、下式で定義されます。

$$
Array Effort = \frac{q^H q}{q_r^H q_r}
$$

ここで、qは入力ベクトルである。
また、入力ベクトルqのときのブライトゾーンの平均音圧を求め、この平均音圧と同じ値になるように、全入力点から同位相で入力したときに入力ベクトルがq_rである。

2. Acoustic Contrast

Acoustic Contrastはブライトゾーンとダークゾーンの平均音圧の比です。式ではΣオートパワーの比になります。このとき、ブライトゾーンの制御点(応答点)をL_Bとし、ダークゾーンの制御点をL_Dとすると、この制御点の数が多いほどΣオートパワーが大きくなってしまうので、制御点数で平均化します。

$$
Acoustic Contrast = \frac{\frac{p_B^H p_B}{L_B}}{\frac{p_D^H p_D}{L_D}}=\frac{L_D q^H G_B^H G_B q}{L_B q^H G_D^H G_D q}
$$

Acoustic Contrastの値が高いほど、ブライトゾーンとダークゾーンの音圧差が大きいので、音場制御がうまくいっていると言えます。

3. ブライトゾーン内の音圧のばらつき

Acoustic Contrastはブライトゾーンとダークゾーンの平均音圧の比であると述べました。そのため、極端な例ですが、ブライトゾーンのある一か所だけの音圧が高く、それ以外の音圧は0であっても、Acoustic Contrastが高くなる可能性があります。この場合、音場制御はうまくいっているとは言えません。そこで、ブライトゾーン内の音圧のばらつきσが指標として扱われることがあります。

$$
σ = (\frac{1}{L_B-1}\sum_{n=1}^{L_B} (10\log{10} (\frac{(p_B)_n^* (p_B)_n}{p_{ave}^* p_{ave}})^2))^{1/2}
$$

スポンサーリンク

MATLABでの検証

では、プログラムを組んで理論を検証しましょう。
今回はArray EffortとAcoustic Contrastの2つ指標を求めました。なお、ここではただの逆行列計算と、最小二乗法を用いた最適化手法を比較しました。結果が下図になります。

細かくは分析しませんが、”最小二乗法を用いた最適化手法”ではArray Effortが低下できています。一方Acoustic Contrastは”ただの逆行列計算”よりも小さい周波数帯域が広いですね。

個人的にAcoustic Contrastが高くできない音場制御方法は価値がないと思っています。それは「別々の音場を提供する」という目的が達成できないからです。

では、Acostic Contrastを最大化する音場制御方法はないのでしょうか?
それがあるのです。
その名も”Acoustic Contrast Optimization”です。そのままの名前ですね。次回は”Acoustic Contrast Optimization”について説明したいと思います。

 

 

今回はこのへんでGood luck

スポンサーリンク

clear all;close all;clc


% % % % % % 初期設定
x=-2:0.1:2;
y=-2:0.1:2;
rho=1.293;
A=1;
c=340; %音速
freq=[10:10:100 200:100:1000 2000:1000:10000]; %周波数
w=2*pi*freq;


% % % % % % 音源位置の設定
x_sp=-1:0.02:1;
y_sp=ones(1,length(x_sp));
input=[x_sp x_sp
    y_sp*(-1.5) y_sp*(+1.5)].';

% % % % % % ブライトゾーンとダークゾーンの設定
bz_geo=[];
dz_geo=[];
bz_x=-0.3:0.1:0.3;
bz_y=-0.3:0.1:0.3;

for ii1=1:length(bz_x)
    for ii2=1:length(bz_y)
        bz_geo=[bz_geo;bz_x(ii1)+0.4 bz_y(ii2)+0.7];
        dz_geo=[dz_geo;bz_x(ii1)-0.4 bz_y(ii2)-0.7];
    end
end
zone=[bz_geo;dz_geo];

Gb=zeros(size(bz_geo,1),size(input,1),length(freq));
Gd=zeros(size(dz_geo,1),size(input,1),length(freq));

for ii1=1:size(bz_geo,1)
    for ii2=1:size(input,1)
        r_bz=sqrt( (bz_geo(ii1,1)-input(ii2,1))^2 + (bz_geo(ii1,2)-input(ii2,2))^2 );
        r_dz=sqrt( (dz_geo(ii1,1)-input(ii2,1))^2 + (dz_geo(ii1,2)-input(ii2,2))^2 );
        if r_bz<0.1
            r_bz=0.1; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
        end
        if r_dz<0.1
            r_dz=0.1; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
        end
        for ii3=1:length(w)
            k=w(ii3)/c;
        Gb(ii1,ii2,ii3)=j*w(ii3)*rho*A./(4*pi*r_bz).*exp(j*(k*r_bz));
        Gd(ii1,ii2,ii3)=j*w(ii3)*rho*A./(4*pi*r_dz).*exp(j*(k*r_dz));
        end
    end
end
ArrayEffort=zeros(2,length(freq));
AcousticContrast=zeros(2,length(freq));
for ii1=1:length(freq)
    G=[Gb(:,:,ii1);Gd(:,:,ii1)];
    q=ones(length(input),1);
    p=G*q;
    ave_p=mean(abs(p));
    
    p_target_def=[ones(size(bz_geo,1),1);zeros(size(dz_geo,1),1)];
    q_def=pinv(G)*p_target_def;
    p_def=G*q_def;
    ave_p_def=mean(abs(p_def));
    r_def=ave_p_def/ave_p;
    
    p_target_ls=G*ones(size(input,1),1);
    p_target_ls(size(bz_geo,1)+1:end,1)=zeros(size(dz_geo,1),1);
    q_ls=inv(G'*G+1*eye(size(G'*G)))*G'*p_target_ls;
    p_ls=G*q_ls;
    ave_p_ls=mean(abs(p_ls));
    r_ls=ave_p_ls/ave_p;
    
    ArrayEffort(1,ii1)=q_def'*q_def/(q'*q*r_def^2);
    ArrayEffort(2,ii1)=q_ls'*q_ls/(q'*q*r_def^2);
    
    AcousticContrast(1,ii1)=size(dz_geo,1)*q_def'*Gb(:,:,ii1)'*Gb(:,:,ii1)*q_def/........
        ( size(bz_geo,1)*q_def'*Gd(:,:,ii1)'*Gd(:,:,ii1)*q_def );
    AcousticContrast(2,ii1)=size(dz_geo,1)*q_ls'*Gb(:,:,ii1)'*Gb(:,:,ii1)*q_ls/........
        ( size(bz_geo,1)*q_ls'*Gd(:,:,ii1)'*Gd(:,:,ii1)*q_ls );
end

figure
subplot(211)
loglog(freq,ArrayEffort)
legend('ただの逆行列計算','最小二乗法を用いた最適化手法')
xlabel('周波数')
ylabel('Array Effort')
subplot(212)
loglog(freq,AcousticContrast)
legend('ただの逆行列計算','最小二乗法を用いた最適化手法')
xlabel('周波数')
ylabel('Acoustic Contrast')

コメント