阵列信号处理及matlab实现_球形麦克风阵列设计
問題背景
要設計球形麥克風陣列(SMA),并保證麥克風在球面等距離分布(可以理解為,球面上任一麥克風與相鄰麥克風的距離相等),需要給出每個麥克風的具體坐標位置。
麥克風坐標計算
球面上的麥克風最終組成一個規則的凸包多面體(參見Platonic solid),麥克風坐標數值解法其實是一個Thomson problem。
此類問題有多種求解思路(可參考cofludy的博客),其中一種是基于力學原理,大致思路如下:
1)首先建立單位球體,然后將待分布的N個麥克風作為帶電粒子,隨機分布在該單位球體的表面;
2) 計算每個帶電粒子在各個方向上受到的合力大小,然后計算出合力在對應帶電粒子上的切線分量;
3) 根據切線分量計算對應帶電粒子沿切向運動飛出該單位球體表面的坐標,然后對每一帶電粒子的坐標沿徑向進行歸一化,使所有帶電粒子再次回到該單位球體的表面;
4) 步驟2)、3)循環若干次,當各帶電粒子所受合力均小于一設定值時,得到各帶電粒子的球面均勾分布,即為N個麥克風在該球形的陣列分布。
詳細的算法步驟可參考專利描述:一種3D錄音系統球面麥克風陣列分布方法【北京大學】
幸運的是,網上有一段現成的MATLAB代碼,可通過鏈接查看,源代碼如下:
function []=main() N=12; %點數量 a=rand(N,1)*2*pi;%根據隨機求面均勻分布,先生成一個初始狀態 b=asin(rand(N,1)*2-1); r0=[cos(a).*cos(b),sin(a).*cos(b),sin(b)]; v0=zeros(size(r0)); G=1e-2;%斥力常數,試驗這個值比較不錯 for ii=1:200%模擬200步,一般已經收斂,其實可以在之下退出[rn,vn]=countnext(r0,v0,G);%更新狀態r0=rn;v0=vn; end plot3(rn(:,1),rn(:,2),rn(:,3),'.');hold on;%畫結果 [xx,yy,zz]=sphere(50); h2=surf(xx,yy,zz); %畫一個單位球做參考 set(h2,'edgecolor','none','facecolor','r','facealpha',0.7); axis equal; axis([-1 1 -1 1 -1 1]); hold off; endfunction [rn vn]=countnext(r,v,G) %更新狀態的函數 %r存放每點的x,y,z數據,v存放每點的速度數據 num=size(r,1); dd=zeros(3,num,num); %各點間的矢量差 for m=1:num-1for n=m+1:numdd(:,m,n)=(r(m,:)-r(n,:))';dd(:,n,m)=-dd(:,m,n);end end L=sqrt(sum(dd.^2,1));%各點間的距離 L(L<1e-2)=1e-2; %距離過小限定 F=sum(dd./repmat(L.^3,[3 1 1]),3)';%計算合力Fr=r.*repmat(dot(F,r,2),[1 3]); %計算合力徑向分量 Fv=F-Fr; %切向分量rn=r+v; %更新坐標 rn=rn./repmat(sqrt(sum(rn.^2,2)),[1 3]); vn=v+G*Fv;%跟新速度 end一開始因為計算機沒有安裝MATLAB只有Python,代碼看著也不多,也有那么一點點Python使用經驗,想著用Python重寫一遍好了。
結果,困難程度遠遠超出預期(請原諒我這只菜鳥~)。一邊查MATLAB中函數的用法,一邊查Numpy中函數的用法,折騰幾天還是沒能完成。不得已,安裝MATLAB后,一步一步確認Python代碼的輸出結果,又折騰幾天最終完成。
Numpy雖然說很多函數和MATLAB非常相似,但有幾個地方需要注意:
- MATLAB下標索引從1開始,Numpy下標索引從0開始;
- MATLAB中多維數組下標依次表示維度 行:列:頁,而Numpy中多維數組下標依次表示維度 頁:行:列;
計算得到球面麥克風點的坐標后,可利用Scipy.spatial中的ConvexHull函數繪制3D凸包,參考以下鏈接,源代碼如下
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy.spatial import ConvexHull# 8 points defining the cube corners pts = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],[0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1], ])hull = ConvexHull(pts)fig = plt.figure() ax = fig.add_subplot(111, projection="3d")# Plot defining corner points ax.plot(pts.T[0], pts.T[1], pts.T[2], "ko")# 12 = 2 * 6 faces are the simplices (2 simplices per square face) for s in hull.simplices:s = np.append(s, s[0]) # Here we cycle back to the first coordinateax.plot(pts[s, 0], pts[s, 1], pts[s, 2], "r-")# Make axis label for i in ["x", "y", "z"]:eval("ax.set_{:s}label('{:s}')".format(i, i))plt.show()計算結果
好了,讓我們一起來看看最終的實現效果吧~
球面均勻分布4點球面均勻分布36點總結
以上是生活随笔為你收集整理的阵列信号处理及matlab实现_球形麦克风阵列设计的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 开发常用工具
- 下一篇: vue获取div中的值_一篇文章看懂Vu