MATLAB.rar_matlab 球型波_wider1g_平面波MATLAB_球面波_球面波matlab

clc clear Qs=2.5*10^-5; f=1000; %声源频率 c=343; %声速 g=1.29; %空气密度 w=2*pi*f; k=w/c; %后面正则化参数的选择会用到 SNR=40; %后面正则化参数的选择会用到 %参考坐标系为（0 0 0） Zs=0.06; %全息面离参考坐标系的距离 Zh=0.03; %重构面离参考坐标系的距离 d1=Zs; %全息面和重建面之间的距离，后面正则化参数的选择会用到 %设声源1中心坐标为（x1,y1,z1） x1=0; y1=0; z1=0; %开始设置全息面坐标，参考坐标系设在与实验室8*8手持阵列的中心等高的地方 s=0.03; %麦克风间距 D=[-0.105 -0.105 Zs -0.075 -0.105 Zs -0.045 -0.105 Zs -0.015 -0.105 Zs 0.015 -0.105 Zs 0.045 -0.105 Zs 0.075 -0.105 Zs 0.105 -0.105 Zs -0.105 -0.075 Zs -0.075 -0.075 Zs -0.045 -0.075 Zs -0.015 -0.075 Zs 0.015 -0.075 Zs 0.045 -0.075 Zs 0.075 -0.075 Zs 0.105 -0.075 Zs -0.105 -0.045 Zs -0.075 -0.045 Zs -0.045 -0.045 Zs -0.015 -0.045 Zs 0.015 -0.045 Zs 0.045 -0.045 Zs 0.075 -0.045 Zs 0.105 -0.045 Zs -0.105 -0.015 Zs -0.075 -0.015 Zs -0.045 -0.015 Zs -0.015 -0.015 Zs 0.015 -0.015 Zs 0.045 -0.015 Zs 0.075 -0.015 Zs 0.105 -0.015 Zs -0.105 0.015 Zs -0.075 0.015 Zs -0.045 0.015 Zs -0.015 0.015 Zs 0.015 0.015 Zs 0.045 0.015 Zs 0.075 0.015 Zs 0.105 0.015 Zs -0.105 0.045 Zs -0.075 0.045 Zs -0.045 0.045 Zs -0.015 0.045 Zs 0.015 0.045 Zs 0.045 0.045 Zs 0.075 0.045 Zs 0.105 0.045 Zs -0.105 0.075 Zs -0.075 0.075 Zs -0.045 0.075 Zs -0.015 0.075 Zs 0.015 0.075 Zs 0.045 0.075 Zs 0.075 0.075 Zs 0.105 0.075 Zs -0.105 0.105 Zs -0.075 0.105 Zs -0.045 0.105 Zs -0.015 0.105 Zs 0.015 0.105 Zs 0.045 0.105 Zs 0.075 0.105 Zs 0.105 0.105 Zs]; %对于给定的麦克风阵列和麦克风放置位置，则（A+）*A确定，重构公式为p=PT*inv((（A+）*A)+q)*((A+)*a) k1=1; for j1=1:size(D,1) for j2=1:size(D,1) r(k1)=sqrt((D(j1,1)-D(j2,1))^2+(D(j1,2)-D(j2,2))^2); if r(k1)==0 B(k1)=1/2+1/(2*k*d1)^2; else J1(k1)=sph_bessel(k,1,r(k1)); %开始积分 syms t f11(k1)=r(k1)*sqrt(1+((t^2)/(2*k*d1)^2)); B11(k1)=IntGaussLager('sph_bessel(k,0,f11(k1))',5); digits(10); B1(k1)=vpa(B11(k1)); B(k1)=(J1(k1)/(k*r(k1)))+(1/(2*k*d1)^2)*B1(k1); end k1=k1+1; end end for k2=1:size(D,1) BA(k2,:)=B(1+(k2-1)*64:k2*64); end %单极子球的声场公式为P=-(i*w*g*a^Qs)*exp(i*k*r)/(r*4*pi) %a为脉动球半径，r为空间任意一点到声源中心的距离 %开始计算全息面的理论声压 for i5=1:size(D,1) r1(i5)=sqrt((D(i5,1)-x1)^2+(D(i5,2)-y1)^2+(D(i5,3)-z1)^2); P1(i5,1)=-(i*w*g*Qs)*exp(i*k*r1(i5))/(r1(i5)*4*pi); end %开始设置重建面尺寸,22*22个点麦克风间距0.01米，也就是1厘米 for i0=1:8 for i10=1:8 cell0{i0,i10}=[-0.105+(i10-1)*0.03 -0.105+(i0-1)*0.03]; %将22*22点的x,y坐标存在cell0里，cell0的维数22*22 end %end i9=1; for i11=1:8 for i12=1:8 cell11{i9,1}=cell0{i11,i12}; %将22*22点的坐标按从左到右从上到下的顺序展开成484维数的列cell11 i9=i9+1; end end for i11=1:8*8 C1(i11,:)=cell11{i11,1}; %将cell转换成矩阵维数484*2,484代表点数，2代表x,y end for i12=1:8*8 C2(i12,1)=Zh; %将重构距离作为重构面的第3列 end C=[C1 C2]; %下面开始计算（A+）*a,对于不同的重构点则有不同的（A+）*a k3=1; for j3=1:size(C,1) for j4=1:size(D,1) r(k3)=sqrt((D(j4,1)-C(j3,1))^2+(D(j4,2)-C(j3,2))^2); if r(k3)==0 B4(k3)=1/(k^2*(C(j3,3)+d1)^2); else B22(k3)=IntGauss('exp(-i*k*cos(m)*(C(j3,3)-d1))*besselj(0,k*(r(k3)*sin(m)))*sin(m)*cos(m)',0,pi/2,3); B2(k3)=eval(B22(k3)); %开始无穷的定积分 B33(k3)=IntGaussLager('besselj(0,k*(r(k3)*sqrt(1+((t^2)/(k^2*(C(j3,3)+d1)^2))))))*t',4); B3(k3)=eval( B33(k3)); B4(k3)=B2(k3)+(1/(k^2*(C(j3,3)+d1)^2))*B3(k3); end k3=k3+1; end end for k4=1:size(D,1) BB(k4,:)=B4(1+(k4-1)*size(C,1):k4*size(C,1)); end %下面开始计算正则化参数q q=(1/2+(1/((2*k*d1)^2)))*10^(-(SNR/10)); %d1为全息面和重构面之间的距离，SNR为包含所有误差和噪声的信噪比 %重构面重构声压 I=eye(size(D,1)); %I为单位对焦矩阵 P=(P1.'*inv(BA+q*I)*BB).'; %重构面的理论声压计算 for i14=1:size(C,1) r2(i14)=sqrt((C(i14,1)-x1)^2+(C(i14,2)-y1)^2+(C(i14,3)-z1)^2); P2(i14,1)=-(i*w*g*Qs)*exp(i*k*r2(i14))/(r2(i14)*4*pi); end %P_cha=abs(P2)-abs(P); P_cha=abs(P2-P); P_cha_2norm=norm(P_cha); P2_2norm=norm(P2); Error1=100*abs(P_cha_2norm/P2_2norm) %开始画图 P21=(abs(P2))'; for ii13=1:8 X(ii13)=-0.105+(ii13-1)*0.03; end for ii14=1:8 Y(ii14)=-0.105+(ii14-1)*0.03; end for ii8=1:8 Z(ii8,:)=P21(1+(ii8-1)*8:8*ii8); end figure surfc(X,Y,Z)%三维曲面 title('重建面理论声压') PH1=(abs(P))'; for iii13=1:8 X(iii13)=-0.105+(iii13-1)*0.03; end for iii14=1:8 Y(iii14)=-0.105+(iii14-1)*0.03; end for iii8=1:8 Z(iii8,:)=PH1(1+(iii8-1)*8:8*iii8); end figure surfc(X,Y,Z)%三维曲面 title('重建面重构声压')