shujuchuli.m

双十字阵矢量水听器阵的指向性研究

clear all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%读数据%%%%%%%%%%%%%%%%%%%%%%%%
[filename, pathname] = uigetfile('*.wav', 'Pick an M-file');
fid=fopen([pathname filename]);
N1=24;%采集时通道数；
frewind(fid);
fseek(fid,42,-1);
sn=fread(fid,12.8*2.56*1000*N1,'int16');
len=length(sn);t1=(0:(len-1)/N1);%/25.6/2.56/1000;
x1=reshape(sn,N1,len/N1)*0.00043524;
c=1500;      %声速 
r1=0.226;    %大子阵半径
r2=0.113;    %小子阵半径
f0=8000;     %信号频率
w=2*pi*f0; 
N=256;     
fs=12.8*2.56*1000;   
wn1=7000*2/fs;wn2=9000*2/fs;    
b=fir1(N,[wn1,wn2]);%Window-based finite impulse response filter design 阶数　上下截止频率  %上下截止频率
nn=1:len/N1;
len1=length(nn);
t=0:(len1-1);
for i=1:24
 y(i,:)=filter(b,1,x1(i,:));%过滤ｘ１，ｂ和１是滤波器的分子和分母　　得到实数的接收数据
end
b=15500;
g=1500;
for i=1:24
    x(i,:)=hilbert(y(i,:));   
end
% x11=hilbert(x(1,:));
% x=[x1;x11];
p1=[x(2,b:b+g);x(5,b:b+g);x(8,b:b+g);x(11,b:b+g);x(14,b:b+g);x(17,b:b+g);x(20,b:b+g);x(23,b:b+g)];
vx=[1/w*x(3,b:b+g);1/w*x(6,b:b+g);1/w*x(9,b:b+g);1/w*x(12,b:b+g);1/w*x(15,b:b+g);1/w*x(18,b:b+g);1/w*x(21,b:b+g);1/w*x(24,b:b+g)];
vy=[1/w*x(4,b:b+g);1/w*x(7,b:b+g);1/w*x(10,b:b+g);1/w*x(13,b:b+g);1/w*x(16,b:b+g);1/w*x(19,b:b+g);1/w*x(22,b:b+g);1/w*x(1,b:b+g)];
for n=1:8;   
      if n==1
          p1(n,:)=p1(n,:)*exp(j*(0)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(26.9446)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(66.8017)*pi/180);
      elseif n==2
          p1(n,:)=p1(n,:)*exp(j*(9.5206)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(9.5206+14.3652)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(9.5206+0.5682)*pi/180);
      elseif n==3
          p1(n,:)=p1(n,:)*exp(j*(29.5568)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(29.5568+39.8564)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(29.5568-1.1615)*pi/180);
      elseif n==4
          p1(n,:)=p1(n,:)*exp(j*(47.659)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(47.659+68.4397)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(47.659+70.8330)*pi/180);
      elseif n==5
          p1(n,:)=p1(n,:)*exp(j*(0)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(25.5726)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(30.2295)*pi/180);
      elseif n==6
          p1(n,:)=p1(n,:)*exp(j*(-13.881)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(86.4930-13.881)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(54.1803-13.881)*pi/180);
      elseif n==7
          p1(n,:)=p1(n,:)*exp(j*(-15.7086)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(-15.7086+67.5710)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(-15.7086+65.3769)*pi/180);
      elseif n==8
          p1(n,:)=p1(n,:)*exp(j*(9.525)*pi/180);
          vx(n,:)=vx(n,:)*exp(j*(9.525-55.0714)*pi/180);
          vy(n,:)=vy(n,:)*exp(j*(13.5488+9.525)*pi/180);
      end
end
Ns=g+1;%采样点数 
N=8;
ph1=-180;                   
ph2=180;                   
phn=1;  
ax=(ph2-ph1)/phn+1;
     Xp1=zeros(ax,Ns);  
     Xp0=zeros(ax,Ns); 
     Xvx01=zeros(ax,Ns);
     Xvy01=zeros(ax,Ns);
     Xvx02=zeros(ax,Ns);
     Xvy02=zeros(ax,Ns);
     Xvx1=zeros(ax,Ns);
     Xvy1=zeros(ax,Ns);
     Xvx2=zeros(ax,Ns);
     Xvy2=zeros(ax,Ns);
     Xvc1=zeros(ax,Ns);
for i=1:ax;     
    num=ph1+(i-1)*phn;
    fain=num*pi/180;  
    for n=1:N;
        a=2;
      if mod(n,a)~=0 
        Xp0(i,:)=p1(n,:)*exp(j*w*r1*cos(fain-(n-1)*2*pi/N)/c); 
        Xvx01(i,:)=vx(n,:)*exp(j*w*r1*cos(fain-(n-1)*2*pi/N)/c);
        Xvy01(i,:)=vy(n,:)*exp(j*w*r1*cos(fain-(n-1)*2*pi/N)/c);
      else
        Xp0(i,:)=p1(n,:)*exp(j*w*r2*cos(fain-(n-1)*2*pi/N)/c); 
        Xvx02(i,:)=vx(n,:)*exp(j*w*r2*cos(fain-(n-1)*2*pi/N)/c);
        Xvy02(i,:)=vy(n,:)*exp(j*w*r2*cos(fain-(n-1)*2*pi/N)/c);
      end
      Xp1(i,:)=Xp1(i,:)+Xp0(i,:);
      Xvx1(i,:)=Xvx1(i,:)+Xvx01(i,:);
      Xvy1(i,:)=Xvy1(i,:)+Xvy01(i,:);
      Xvx2(i,:)=Xvx2(i,:)+Xvx02(i,:);
      Xvy2(i,:)=Xvy2(i,:)+Xvy02(i,:);
      Xvc01(i,:)=Xvx1(i,:)*cos(fain)+Xvy1(i,:)*sin(fain);
      Xvc02(i,:)=Xvx2(i,:)*sin(fain-pi/4)+Xvy2(i,:)*cos(fain-pi/4);
      Xvc1(i,:)=Xvc01(i,:)+Xvc02(i,:);
    end
   Xp(i)=(sum(abs(Xp1(i,:)').^2));
   Xvc(i)=(sum(abs(Xvc1(i,:)').^2));
   Rp(i)=(Xp(i)+Xvc(i))*Xvc(i);
end 
i=1:ax;
  num=ph1+(i-1)*phn;
  fain=num*pi/180;
maxRp=max(Rp);
Rpp=10*log10(Rp(i)./maxRp);
figure
plot(ph1+(i-1)*phn,Rpp);
legend('交叉矢量阵(p+vc)*vc',4);
axis([-180 180 -40 0])
grid on
hold on
xlabel('角度')
ylabel('dB')