奇异值分解波束形成MATLAB程序

close all

clear all

clc

c=1500;

fs=10000;

f=500;

w=2*pi*f;

k=w/c;

M=13; %阵元个数

Nmid=1; %参考点

d=3;%阵元间距

m=[0:1:M-1];

yi=zeros(M,1);% 返回一个M*1维的零矩阵

zi=zeros(M,1);

xi=m*d;

xi=xi.';

%各阵元坐标

y1=12;

x1=10;z1=10;% 声源位置， y轴指向声源平面

Ric1=sqrt((x1-xi).^2+(y1-yi).^2+(z1-zi).^2);%声源至各阵元的距离 M*1维

Rn1=Ric1-Ric1(Nmid);%声源至各阵元与参考阵元的声程差矢量 M*1维

T=1;

t=1/fs:1/fs:T;

L=length(t);%时间长度

s1=exp(j*w*t);%参考阵元接收到的信号 1*L维

p1=zeros(M,L);%M*L维

for k1=1:M

tao1(k1)=Rn1(k1)/c;

p1(k1,:)=Ric1(Nmid)/Ric1(k1)*s1.*exp(-j*w*tao1(k1)); %各阵元经过幅度衰减和相位延迟后接收到的信号，M*L维

end

SNR=20;%信噪比

Am=10.^(-SNR/20);

n1=Am*(randn(M,L)+j*randn(M,L));%各阵元接收的噪声矢量 M*L维

p=p1+n1;%各阵元接收的声压信号矩阵 M*L

NFFT=2048;

for k1=1:M

P(k1,:)=fft(p(k1,:),NFFT);

end

% ----------------------------------------------------------

% 扫描范围

step_x=0.5;

step_z=0.5;

y=y1;

x=[0:step_x:20];

z=[0:step_z:20];

for k1=1:length(z) % 纵坐标

for k2=1:length(x)

Ri=sqrt((x(k2)-xi).^2+(y-yi).^2+(z(k1)-zi).^2); %该扫描点至各阵元的聚焦距离矢量

Rn=Ri-Ri(Nmid);%扫描点到各阵元与参考阵元的程差矢量 M*1

tao=Rn/c;

for k3=1:M

PP(k3,:)=P(k3,:).*exp(j*w*tao(k3,:));

end

for i = 2:M

A(i-1,:) = 2*[real(PP(i,:)-PP(i-1,:)),imag(PP(i,:)-PP(i-1,:))];

end

E(k1,k2) = 1/trace(A*A.');

end

end

% 归一化

E=E/max(max(E));

figure;

mesh(x,z,E);

xlabel('x(m)'),ylabel('z(m)')

title('奇异值分解单亮点声图')

figure;

pcolor(x,z,E);shading interp;

xlabel('x(m)'),ylabel('z(m)')

title('奇异值分解单亮点声图')

colorbar

[a,b] = find(E==max(max(E)));

xm1=b/length(x)*20, zm1=a/length(z)*20,

[a b]=find(E > 0.5*max(max(E)));

deta_x=(max(b)-min(b))/length(x)*20

deta_z=(max(a)-min(a))/length(z)*20