%傳統波束形成，CBF （Ps：這個程序是別人的，不是我寫的，但是具體是在哪里找到的已經忘了）
clear all; close all; clc;
%---------初始化常量----------%
c = 334; % 聲速c
fs = 1000; % 抽樣頻率fs
T = 0.1; %
t = 0:1/fs:T; % 時間 [0,0.1]
L = length(t); % 時間長度，采樣總點數：101個
f = 500; % 感興趣的頻率，信號自己的頻率
w = 2*pi*f; % 角頻率 k = w/c; % 波數 k
M = 18; % 陣元個數,實際只有17個
%接下來是陣元坐標，聲源位置，這里設置的是（12,10,12）點，同時設置陣元的高斯白噪聲。
%---------各陣元坐標---------% %
Nmid = 12; % 參考點 %
d = 3; % 陣元間距 %
m = (0:1:M-1) yi = zeros(M,1); % 生成一個M*1維的零矩陣
zi = [0;3;6;9;12;15;18;21;24;12;12;12;12;12;12;12;12;12];
xi = [12;12;12;12;12;12;12;12;12;0;3;6;9;12;15;18;21;24]; %形成柿子麥克風陣列坐標 %xi = xi.' % 列向量 m*d 陣元數*陣元間距
figure(1)
plot(xi,zi,'r*'); title('十字形麥克風陣列') %--------- 聲源位置----------% x1 = 12; y1 = 10; z1 = 12; %聲源位置 （12,10,12） x,z為水平面
x2 = 12; y2 = 0; z2 = 12; %參考點：（12，0，12）
Ric1 = sqrt((x1-xi).^2+(y1-yi).^2+(z1-zi).^2); % 聲源到各陣元的距離
Ric2 = sqrt((x1-x2).^2+(y1-y2).^2+(z1-z2).^2); %聲源到參考點的距離
Rn1 = Ric1 - Ric2; %聲源至各陣元與參考陣元的聲程差矢量
s1 = cos(2*w*t); % 參考陣元接收到的矢量
Am = 10^(-1); % 振幅
n1 = Am * (randn(M,L)+j*randn(M,L)); % 各陣元加上高斯白噪聲，考慮復數18*101的陣列
p1 = zeros(M,L);
%3、延遲求和----------------------------------------------------------------------------
%整個程序最關鍵的部分，延遲求和，同時得到各陣元接收的聲壓信號矩陣。以及協方差矩陣，這個還有疑問，要把論文讀懂來理解。
%--------------------------各陣元的延遲求和--------------------------------%
for k1 = 1:M
p1(k1,:) = Ric2/Ric1(k1)*s1.*exp(-j*w*Rn1(k1)/c); %Ric2/Ric1(k1)等於把每一個點的距離歸一化.exp()表示時延 18*101的矩陣 % 接收到的信號
%時間和陣元的矩陣
end
p = p1+n1; % 各陣元接收的聲壓信號矩陣 R = p*p'/L; % 接收數據的自協方差矩陣 A.'是一般轉置，A'是共軛轉置。 %%%R代表幅值嗎？ %18*18
%4、掃描整個聲源平面------------------------------------------------------------------------
%我們設置步長為0.1，掃描范圍是20x20的平面，雙重for循環得到M*1矢量矩陣，最后得到交叉譜矩陣（cross spectrum matrix）。由DSP理論，這個就是聲音的功率。 %-------掃描范圍------%
step_x = 0.1; % 步長設置為0.1 %和線陣里面掃描角度是一個道理
step_z = 0.1;
y = y1; %y=10
y1=10 x = (9:step_x:15); % 掃描范圍 9-15
z = (9:step_z:15);
for k1=1:length(z) %先掃x軸，再掃z軸，設掃描到的點為（x=k2，z=k1），y=10 k1=61，k2=61
for k2=1:length(x)
Ri = sqrt((x(k2)-xi).^2+(y-yi).^2+(z(k1)-zi).^2); %（x(k2),10,z(k1)）走到的點距離陣元的距離 %每走一步，產生一個18*1的矩陣
Ri2 = sqrt((x(k2)-x2).^2+(y-y2).^2+(z(k1)-z2).^2); %走到的點距離參考點的距離 %一個數 % 該掃描點到各陣元的聚焦距離矢量
Rn = Ri-Ri2; %每一步產生一個1*18的矩陣 % 掃描點到各陣元與參考陣元的程差矢量
b = exp(-j*w*Rn/c); % 聲壓聚焦方向矢量 %歸整到陣元，18*1的矩陣，消除時間的影響
Pcbf(k1,k2) = abs(b'*R*b); % CSM
end
end %聲音的功率
%5、歸一化處理
%歸一化處理的程序 %--------------------------------------歸一化------------------------------%
for k1 = 1:length(z);
pp(k1) = max(Pcbf(k1,:)); % Pcbf 的第k1行的最大元素的值
end
Pcbf = Pcbf/max(pp); % 所有元素除以其最大值 歸一化幅度
%6、作圖
%觀察得到的結果 %-------------------------------作圖展示-----------------------------------%
figure(2)；
surf(x,z,Pcbf); xlabel('x(m)'),ylabel('z(m)') ；title('CBF三維單聲源圖')； colorbar
figure(3) ；
pcolor(x,z,Pcbf); shading interp; xlabel('x(m)'); ylabel('z(m)'); title('CBF單聲源圖') ；colorbar
%7、結果
%十字型陣列最終得到的結果效果並不理想，沒有達到一個點聲源的理想結果
剛開始這個程序其實沒怎么看懂，於是在網上找了CBF的知識，找到一篇可以找到參考的論文。【水聲定位算法之CBF波束形成_360doc個人圖書館】http://www.360doc.cn/mip/558928696.html
自己照着這個程序寫了一個線陣的，可是寫了一天都沒有寫出來，后面的有些地方有點模糊，不知道怎么處理，以后再思考一下吧。