MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて4 （Acoustic Contrast Optimization）

2021.04.10

前回の復習

前回の記事の内容を確認してから本記事を読むことをお勧めいたいします。

MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて1
本記事ではフィードフォワードの音場制御方法の理論について、数式とMATLABのコードを交えながら説明いたします。理論としてはすごくシンプルなので、高校生でも理解できると思います。
mech-eng-uni.com
2021.04.01
MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて2 （最小二乗法を用いた最適化）
最小二乗法を用いた音場最適化制御方法について、理論式やMATLABプログラムを示しながら説明します。ここまで簡単に詳しく説明してくれている文献はほかにないので、ぜひ読んでみてください。
mech-eng-uni.com
2021.04.04
MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて3 （効果的な目標音場の設定）
最小二乗法を用いた音場最適化制御方法について、理論式やMATLABプログラムを示しながら説明します。この音場最適化手法では、目標音場の設定が重要となります。この目標音場の設定方法について詳しく説明いたします。ここまで簡単に詳しく説明してくれている文献はほかにないので、ぜひ読んでみてください。
mech-eng-uni.com
2021.04.06
MATLABで学ぶ信号処理・音響工学　音場制御の評価指標について1
音場制御では、制御がうまくいっているか、ブライトゾーンとダークゾーンの構築がうまくいっているかを評価するための指標があります。本記事では、Array Effort、Acoustic Contrast、ブライトゾーン内の音圧のばらつき、について説明し、MATLABのプログラムを示しながら、音場制御方法の違いと評価指標の差異を確認します。
mech-eng-uni.com
2021.04.08

スポンサーリンク

Acoustic Contrast Optimization

“Acoustic Contrast Optimization”はAcoustic Contrastを最大にする音場最適化手法です。前の記事の復習になりますが、Acoustic Contrastは下式になります。

$$
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}
$$

最適化を計算するためには、評価関数というものを定義します。その評価関数は下式になります。ちなみに評価関数J_ACはスカラーです。

$$
J_{AC} = q^H G_D^H G_D q + λ_{C}(q^H G_D^H G_D q -B) + λ_{M}(q^H q – E_m)
$$

BはブライトゾーンのΣパワースペクトルに対する制約条件です。また、Emはシステムが必要とする入力パワーを制約するための定数でしたね。λ_Mやλ_Cは重みづけをするためのパラメーターだと思ってください。

線形最適化では「評価関数を設計変数の偏微分　＝　0」を解けばよいので、入力ベクトルqで偏微分して下式を導きます。

$$
λ_{C}q= (G_D^H G_D)^{-1} (G_D^H G_D + λ_{M}I)q
$$

学生

学生

「q=」の形にできないやん！

という感じですよね。

ただ、この等式は解くことができるようです。
（詳細は割愛しますが、偉い先生方ががんばってといたっぽいです。）

最適な入力ベクトルは下式の最小固有値の固有ベクトルに”比例”することがわかっています。「”比例”すること」という意味は、フィルターを作成するのに必要であるスピーカー間の相対的な振幅比と位相差はわかるということです。

$$
(G_D^H G_D)^{-1} (G_D^H G_D + λ_{M}I)
$$

MATLABでの検証

では、プログラムを組んで理論を検証しましょう。
今回は３つの音場制御手法（単純な逆行列、LS手法、本手法）で、2つの評価指標（Array EffortとAcoustic Contrast）を比較します。また、今回は下のような音響マップは表示しないことにします。理由は、音響マップでは最適化手法を定量的に評価できないからです。（音響マップでは良さそう、悪そう、といった定性的な比較しかできない）

図1　音響マップの例

比較結果を図2に示します。

図2　比較結果

”Acoustic Contrast Optimization”が広帯域でArray Effortが低く、Acoustic Contrastが高いことがわかります。この3つの手法ですと、”Acoustic Contrast Optimization”が最も音場制御に適した最適化だと考えられます。

（Acoustic Contrastを最大にする手法なので、Acoustic Contrastが高くなるのはあたりまえですが。。。）ただ、10kHz近傍でLS法の方が高い帯域がありますね。これは、固有値計算が特異値に近いことが原因でうまく計算できていない可能性があります。ただ、この辺の課題も解決するテクニック（ノウハウ）があるので、機械があれば記事にしたいと思います。

 

次回はイギリスのUniversity of Surreyが研究していたPlanarityという最適化手法について説明できればと思っています。ただ、Planarityは説明するのが難しいので、違うことを記事にするかもしれません。（いずれにしても、いつかは記事にしたいと思います。）

今回はこのへんで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
lambdaM=0;
ArrayEffort=zeros(3,length(freq));
AcousticContrast=zeros(3,length(freq));
for ii1=1:length(freq)
    Gbtemp=Gb(:,:,ii1);Gdtemp=Gd(:,:,ii1);
    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;
    
    [q_ac d_ac]=eigs(inv(Gbtemp'*Gbtemp)*(Gdtemp'*Gdtemp+lambdaM),1,'sm') ;
    p_ac=G*q_ac;
    ave_p_ac=mean(abs(p_ac));
    r_ac=ave_p_ac/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);
    ArrayEffort(3,ii1)=q_ac'*q_ac/(q'*q*r_def^2);   
    
    AcousticContrast(1,ii1)=size(dz_geo,1)*q_def'*Gbtemp'*Gbtemp*q_def/........
        ( size(bz_geo,1)*q_def'*Gdtemp'*Gdtemp*q_def );
    AcousticContrast(2,ii1)=size(dz_geo,1)*q_ls'*Gbtemp'*Gbtemp*q_ls/........
        ( size(bz_geo,1)*q_ls'*Gdtemp'*Gdtemp*q_ls );
    AcousticContrast(3,ii1)=size(dz_geo,1)*q_ac'*Gbtemp'*Gbtemp*q_ac/........
        ( size(bz_geo,1)*q_ac'*Gdtemp'*Gdtemp*q_ac );
end

figure
subplot(211)
loglog(freq,ArrayEffort)
legend('ただの逆行列計算','Leaset Squares optimization','Acoustic Contrast Optimization')
xlabel('周波数')
ylabel('Array Effort')
subplot(212)
loglog(freq,AcousticContrast)
legend('ただの逆行列計算','Leaset Squares optimization','Acoustic Contrast Optimization')
xlabel('周波数')
ylabel('Acoustic Contrast')

コメント