Sonar_Exp_beamform.zip_时域_时域波束形成_相移_频域_频域 beamform

% clc; clear;close all; %% 参数 c=1500; f0=1.5e3; fs=150e3; %100倍采样 T=50e-3; SNR=20; lamda=c./f0; D=lamda./2; %间距 0.5m N=21; maxdelay=(N-1).*D./c; %端射方向 最大延时时间 Nmaxdelay=round(maxdelay.*fs); %最大延时点 取近似 为1000 alpha=(-90:90)*pi./180; %扇面 不会出现左右模糊 px=(0:N-1).*D; %阵的位置 theta0=10*pi./180; %来向方向 N_delay0=D.*sin(theta0)./c.*fs; % 不同角度 相邻阵元的延时点数 Wn=2*[1000 2000]/fs; %滤波器 学 1000-2000 bb=fir1(128,Wn,'bandpass'); %FIR带通滤波器1-2k系数 带限白噪声 % noise=filter(bb,1,NN); % figure();freqz(bb,1); % [h,w]=freqz(bb,1);% 看滤波器 增益图 % figure();plot(w/pi*fs/2,20*log10(abs(h))); %% 信号模型 t=0:1./fs:T; %时间 sig00=sqrt(2).*cos(2.*pi.*f0.*t); for mm=1:N; N0=round(N_delay0.*(mm-1)); %延时点 30度时 最大为500个点 sig0(mm,:)=[zeros(1,Nmaxdelay.*2-N0) sig00 zeros(1,Nmaxdelay.*1+N0)]; %形成每一路信号 延时3N个点 cw信号 end len=length(sig0(1,:)); %信号长度 for mm=1:N; NN=normrnd(0,1,1,len); %% 加噪声 Normal random numbers noise=filter(bb,1,NN); %滤波 带限白噪声 noise=noise/std(noise); %噪声能量归一化 std： 标准差 noise=noise/10^(SNR/20); %噪声与信噪比的转化 sig(mm,:)=sig0(mm,:)+noise; %加噪声 end %% 时域波束形成 tt=(0:len-1)./fs; %tt最大为 0.07 tt1=(Nmaxdelay+1:len-1-Nmaxdelay)./fs; %0.0067到0.063 8500个点 有用的点做波束形成 for mm=1:length(alpha); %扫描 for nn=1:N; tc0=(nn-1).*D.*sin(alpha(mm))./c; %延时差 ss(nn,:)=interp1(tt,sig0(nn,:),tt1-tc0); %interp1 一维数据插入 相当于延时 先到延时 end %插值还插小了 ，补偿回来 总共N个 是对的 sss(mm,:)=sum(ss,1); %求和 对21路数据求和 % sss1(mm,:)=sum(ss(1:20,:),1); % sss2(mm,:)=sum(ss(2:21,:),1); sout(mm)=std(sss(mm,:));%std： 标准差 x=开方（1/N求和【（xi-u）^2）】 相当于平方+积分 % mm end figure plot(alpha.*180./pi,20.*log10(sout./max(sout))) % plot(alpha.*180./pi,sout./max(sout)) xlabel('方位角/度') ylabel('波束输出/分贝') title('时域波束形成') grid on %% 频域波束形成 vscan=exp(-j*2*pi*(f0)/c*px'*sin(alpha)); %补偿相位 for ii=1:N Sf=fft(sig(ii,:),fs); %步骤1 fft R1(ii)=Sf(f0); %对同一序号k的谱线 抽取 end R0=R1*vscan; yy=abs(real(R0)); %取实部 yy=yy./max(yy); figure plot(alpha.*180./pi,20.*log10(yy)) xlabel('方位角/度') ylabel('波束输出/分贝') title('频域波束形成') grid on %% DFT波束形成 for ii=1:N Sf=fft(sig(ii,:),fs); R1(ii)=Sf(f0+1); end nfft=1024; Beam_DFT=fftshift(abs(fft(R1,nfft))); % fft Beam_DFT=20*log10(Beam_DFT./max(Beam_DFT)); %归一化 nn1=[nfft/2+1:nfft,1:nfft/2]; angle = asin(lamda/D*[-nfft/2:nfft/2-1]/nfft)*180/pi; %理解了 -90-90度 figure plot(angle,Beam_DFT) xlabel('方位角/度') ylabel('波束输出/分贝') title('DFT波束形成') grid on % 画在一张图上 figure ;subplot(131); plot(alpha.*180./pi,20.*log10(sout./max(sout))) xlabel('方位角/度');ylabel('波束输出/分贝') title('时域波束形成');grid on subplot(132);plot(alpha.*180./pi,20.*log10(yy)) xlabel('方位角/度');ylabel('波束输出/分贝');title('频域波束形成');grid on subplot(133);plot(angle,Beam_DFT);xlabel('方位角/度') ylabel('波束输出/分贝');title('DFT波束形成');grid on