全息面声压用平面波展开遇到了问题

txl

在自由场中，声波符合亥姆霍兹齐次微分方程，在直角坐标中的解是平面波，分为传播波和倏逝波两种，传播规律比较简单，其线性组合必然满足齐次亥姆霍兹微分方程。但是我在用平面波展开全息面声压是遇到了问题，用MATLAB进行仿真时发现预测面的理论声压数据和预测面的计算声压数据相差很大，请教高手，问题到底出现在哪？计算的程序和结果图片都贴上来请高手解答，非常感谢！ 请教1 - 副本1.docx (42.17 KB, 下载次数: 6)

2017-5-1 22:01 上传

点击文件名下载附件

txl

这是程序：
clear ;
%基本参数设置
m0=1.29;     %空气密度1.29
c0=340;     %空气中声速
v0=2.5;     %点声源表面振动速度
r0=0.001;    %点声源半径
Q0=4*pi*r0.^2*v0;%点声源强度
complex=sqrt(-1);
PI=3.1415926;
f=200;k=2*pi*f/343;
      %分析频率和波数
Zh=0.05;      %全息测量面距离
Z=0.1;        %预测面距离

x0=0;
y0=0;
z0=0;%点声源坐标

%全息面的离散化
Lx=1;Ly=1; %全息面大小
M=8; N=8;      %x方向和y方向的阵元数

delta_X=Lx/(M-1);
delta_Y=Ly/(N-1);定义采样间隔
xl=-0.5:delta_X:0.5;
yl=-1*(-0.5:delta_Y:0.5);
[Xp,Yp]=meshgrid(xl,yl);
%采样点的坐标值

%全息面声压数据和预测面理论声压数据
RZh=sqrt((Xp-x0).^2+(Yp-y0).^2+Zh.^2);%全息测量面半径
RZ=sqrt((Xp-x0).^2+(Yp-y0).^2+Z.^2);%预测面半径
ph=complex*m0*c0*k*Q0*exp(-complex*k*RZh)./(4*pi.*RZh);%全息面声压
pz=complex*m0*c0*k*Q0*exp(-complex*k*RZ)./(4*pi.*RZ);%预测面理论声压
%%%%%%%%全息面声压数据用平面传播波和平面倏逝波的线性组合表示，下面求线性组合的系数，矩阵法求得。
pro_1=zeros(64,64);%用来存放传播波和倏逝波的基向量
z=Zh;
x=-0.5:1/7:0.5;y=-(-0.5:1/7:0.5);%全息数据点的坐标
kx=-3:1:4;ky=-3:1:4;%kx和ky的取值范围，应该是整个实数范围，为了计算方便，只取了部分整数
for nx=1:8
    for ny=1:8
        for lx=1:8
            for ly=1:8
                if kx(lx)^2+ky(ly)^2<k^2
                    kz=sqrt(k^2-kx(lx)^2-ky(ly)^2);
                    pro_1((nx-1)*8+ny,(lx-1)*8+ly)=exp(i*(kx(lx)*x(nx)+ky(ly)*y(ny)+kz*z)); %传播波
                else
                    kz1=sqrt(kx(lx)^2+ky(ly)^2-k^2);
                    pro_1((nx-1)*8+ny,(lx-1)*8+ly)=exp(i*(kx(lx)*x(nx)+ky(ly)*y(ny)))*exp(-1*kz1*z); %倏逝波
                end

            end

        end

    end

end   
pro_1_real=real(pro_1);
pro_1_imag=imag(pro_1); %分别提取平面波矩阵的实数部分矩阵和虚数部分矩阵
pro=[pro_1_real;pro_1_imag];%实数部分矩阵和虚数部分矩阵重新构成矩阵
[U,S,V]=svd(pro);%采用svd方法分解新构成的矩阵
ph=ph';
p=reshape(ph,64,1);
p_real=real(p);
p_imag=imag(p);
p=[p_real;p_imag];%全息面数据构成列向量，并分开实数部分和虚数部分
C=[(pro'*pro)^(-1)]*pro'*p;%采用最小二乘法求传播波和倏逝波的系数
%预测
z_p=0.1;
for nx=1:8
    for ny=1:8
        for lx=1:8
            for ly=1:8
                for nn=1:64
                if kx(lx)^2+ky(ly)^2<k^2
                    kz(nn)=sqrt(k^2-kx(lx)^2-ky(ly)^2);
                    pro_p((nx-1)*8+ny,(lx-1)*8+ly)=exp(i*(kx(lx)*x(nx)+ky(ly)*y(ny)+kz(nn)*z_p)); %预测面，传播波传播规律：幅值不变，相位变化
                else
                    kz1(nn)=sqrt(kx(lx)^2+ky(ly)^2-k^2);
                    pro_p((nx-1)*8+ny,(lx-1)*8+ly)=exp(i*(kx(lx)*x(nx)+ky(ly)*y(ny)))*exp(-1*kz1(nn)*z_p);%预测面，倏逝波传播规律：相位不变，幅值呈指数减小
                end
                end
            end

        end

    end

end

p_p=pro_p*C;  %预测面数据列向量

p_p=reshape(p_p,8,8);
p_p=p_p';    %预测面列向量数据转换成矩阵，与预测点一一对应
figure;
surf(abs(ph));
figure
surf(abs(p_p));
figure
surf(abs(pz));

txl

结果图片在《请教1-副本1》中

txl

怎么没人回答呢，自己顶

eastar

程序运行没问题吧？

txl

eastar 发表于 2017-5-3 08:58
程序运行没问题吧？

程序可以运行出结果，但是运行的结果和理论值对不上去，不知道什么问题

txl

怎么没有人来回答呢

txl

期待解答

Ljq

应该是倏逝波的问题，当kx(lx)^2+ky(ly)^2>k^2的时候，那个语句不能再这样写了，否则会因为指数的递增而掩盖掉真实的重构。详情可以参阅《傅里叶声学》【美】Earl G.Wiliams