%----------------------阵列和接收信号参数------------------------
clc;
clear all;
close all;
c=1500;%海水声速设定最高和最低频率
flow=200;fhigh=400;%设定最高和最低频率
wavelength=c/fhigh;
d=wavelength/2;%声源间隔
M=12;%振源个数
L=d*(M-1);%线阵长度
Bs=floor(c/L)/10;%由窄带条件确定子带间带宽间隔
f=[flow:Bs:fhigh];
%----------------------常规宽带波束图------------------------
theta=-90:1:90-0.3; %搜索范围
theta0=10; %信号源的来波方向
fcbf=zeros(length(f),length(theta));
wc=zeros(M,1);
for fu=1:length(f)
steer=exp(1i*2*pi*f(fu)/c*d*sin(theta0*pi/180)*[0:M-1]');
wc(:,fu)=steer/M;
for nu=1:length(theta)
a=exp(1i*2*pi*f(fu)/c*d*sin(theta(nu)*pi/180)*[0:M-1]');
fcbf(fu,nu)=wc(:,fu)'*a;
end;
end;
fcbf=10*log10(abs(fcbf)/(max(max(abs(fcbf)))));
figure(1)
mesh(theta,f,fcbf);
xlabel('theta');
ylabel('frequency');
title('cbf response for wide band signals');
fs=5*fhigh;%采样率
N=512;%接收数据的长度
T=N/fs;%接收数据的脉宽
L=N/2;%缓存的存储长度
k=21:1:56;% 取出频域数据中k = 21,...56的子带数据,对应频率范围为[0.4102f0,1.0938f0]
fk=zeros(1,length(k));
for nk=1:length(k)%采样频率点的个数
fk(nk)=fs*k(nk)/L;
end
theta1=-30*pi/180;%信号来向
%------------------------接收到的数据--------------------------
delaytime=-[0:M-1]*d*sin(theta1)/c;
delaypoints=floor(delaytime*fs);%延时采样间隔取整数
%仿真经过采样后的线性调频信号sin(2*pi*(f1+(fu-f1)/(2*T)*t(i))*t(i))
te=0:1:N-1;
t=(0:1:N+delaypoints(M)-1)/fs;
receiveLFMtmp=zeros(M,length(t));
receiveLFM=zeros(M,N);
for m=1:M %阵元的循环
for i=1:length(t) %时间的循环
receiveLFMtmp(m,i)=sin(2*pi*(flow+(fhigh-flow)/(2*T)*(t(i)+delaytime(m)))*(t(i)+delaytime(m)));%阵元m接收到的信号sm(t)=s(t-τm)
if t(i)-delaytime(m)<0
receiveLFMtmp(m,i)=0; %采样时间 0<t<T,如不满足作补零处理,即信号没到达时接收到的数据为0
elseif t(i)-delaytime(m)>=T
receiveLFMtmp(m,i)=0;
end
if i<=N
receiveLFM(m,i)=receiveLFMtmp(m,i);%从receiveLMFtmp中取前512个点作为采样数据
end
end
end
figure(2);
for m=1:M
t = 1:1:(N+delaypoints(M));
plot(t,receiveLFMtmp(m,:)/2+m);
hold all;
end
xlabel('\iti');
ylabel('阵元号\itm');
axis([0 N+delaypoints(M)-1 0 13]);
ktheta1=zeros(1,length(k));
atheta1=zeros(M,length(k));
for nf1=1:length(k)
ktheta1(nf1)=2*pi*fk(nf1)/c*sin(theta1);
for m=1:M
atheta1(m,nf1)=exp(-1i*ktheta1(nf1)*(m-1)*d);
end
end
w_c=atheta1/M;%每一个窄带常规波束形成器的加权向量wc(fk)
%------------------------接收到的数据的功率谱分析---------------------------------------
figure(3);
Hs = spectrum.periodogram;
Hpsd = psd(Hs,receiveLFM(1,:),'Fs',fs/fhigh); % fs/fhigh归一化处理
F_2 = Hpsd.Frequencies;
PSD_2 = Hpsd.Data;
plot(F_2,10*log10(PSD_2),'LineWidth',2);
xlabel('频率({\itf}_0)');
ylabel('功率谱/(dB/Hz)')
title('Periodogram周期图法估计功率谱密度');
axis([0 fs/fhigh/2 -50 10]);
grid on;
%-----------------------------DFT-Beamforming no overlap— 无重叠---------------
%无重叠分段
xnseg1= zeros(L,M,2);%(时间点,阵元编号,分段号)
Xkseg1= zeros(L,M,2);
for m=1:M
for i=1:L
xnseg1(i,m,1) = receiveLFM(m,i);
xnseg1(i,m,2) = receiveLFM(m,i+L);
end
Xkseg1(:,m,1) = fft(xnseg1(:,m,1));
Xkseg1(:,m,2) = fft(xnseg1(:,m,2));
end
Ykseg1=zeros(2,length(k));%频域个窄带内的波束形成
for i=1:length(k)
Ykseg1(1,i)=w_c(:,i)'*Xkseg1(k(i)+1,:,1).'; % k(i)+1 才对应X(21)→X(56)
Ykseg1(2,i)=w_c(:,i)'*Xkseg1(k(i)+1,:,2).';
end
%2段波束输出频域数据进行IDFT,得到时域输出
Ykseg1_full = zeros(L,2);
for i=1:length(k)
Ykseg1_full(k(i)+1,1) = Ykseg1(1,i); % k(i)+1 才对应Y(21)→Y(56)
Ykseg1_full(L-k(i)+1,1) = conj(Ykseg1(1,i)); %L为存储的数据长度,这里L-k(i)+1对应的负频率
Ykseg1_full(k(i)+1,2) = Ykseg1(2,i);
Ykseg1_full(L-k(i)+1,2) = conj(Ykseg1(2,i));
end
%其他频段直接置零,进行截断
ynseg1_full = zeros(L,2);
ynseg1_full(:,1) = ifft(Ykseg1_full(:,1));%作傅里叶反变换
ynseg1_full(:,2) = ifft(Ykseg1_full(:,2));
%数据无重叠进行拼接
ynseg1 = zeros(1,N);
for i=1:L
ynseg1(i)=ynseg1_full(i,1);
ynseg1(i+L)=ynseg1_full(i,2);
end
figure(4)
realsignals=receiveLFM(1,:);
plot(te,real(ynseg1),'r',te,real(realsignals),'b');
legend('received signals using CBF','real signals');
title('无重叠下的常规波数形成接收信号');
%-----------------------------DFT-Beamforming 有重叠—---------------------------------
%阵元有效数据分段进行DFT:50%重叠分为3段
xnseg2 = zeros(L,M,3);
Xkseg2 = zeros(L,M,3);
for m=1:M
for i=1:L
xnseg2(i,m,1) = receiveLFM(m,i);
xnseg2(i,m,2) = receiveLFM(m,i+L/2);
xnseg2(i,m,3) = receiveLFM(m,i+L);
end
Xkseg2 (:,m,1) = fft(xnseg2(:,m,1));
Xkseg2 (:,m,2) = fft(xnseg2(:,m,2));
Xkseg2 (:,m,3) = fft(xnseg2(:,m,3));
end
Ykseg2 = zeros(3,length(k));
for i=1:length(k)
Ykseg2(1,i) = w_c(:,i)'*Xkseg2 (k(i)+1,:,1).';
Ykseg2(2,i) = w_c(:,i)'*Xkseg2 (k(i)+1,:,2).';
Ykseg2(3,i) = w_c(:,i)'*Xkseg2 (k(i)+1,:,3).';
end
Ykseg2_full = zeros(L,3);
for i=1:length(k)
Ykseg2_full(k(i)+1,1) = Ykseg2(1,i);
Ykseg2_full(L-k(i)+1,1) = conj(Ykseg2(1,i));
Ykseg2_full(k(i)+1,2) = Ykseg2(2,i);
Ykseg2_full(L-k(i)+1,2) = conj(Ykseg2(2,i));
Ykseg2_full(k(i)+1,3) = Ykseg2(3,i);
Ykseg2_full(L-k(i)+1,3) = conj(Ykseg2(3,i));
end
ynseg2_full = zeros(L,3);
ynseg2_full(:,1) = ifft(Ykseg2_full(:,1));
ynseg2_full(:,2) = ifft(Ykseg2_full(:,2));
ynseg2_full(:,3) = ifft(Ykseg2_full(:,3));
% 数据按照50%重叠进行拼接,顺序不能变动,而且不能在一个循环里完成
ynseg2 = zeros(1,N);
for i=1:L
ynseg2(i) = ynseg2_full(i,1);
end
for i=1:L
ynseg2(i+L/2) = ynseg2_full(i,2);
end
for i=1:L
ynseg2(i+L) = ynseg2_full(i,3);
end
figure(5)
plot(te,real(ynseg2),'r',te,real(realsignals),'b');
legend('received signals using CBF','real signals');
title('50%重叠下的常规波数形成接收信号');
figure(6)
plot(te,real(ynseg1)-real(realsignals),'--k','LineWidth',2);
hold on;
plot(te,real(ynseg2)-real(realsignals),'-b','LineWidth',2);
xlabel('\iti');
ylabel('{\ity}({\iti})-{\its}({\iti})');
legend('不重叠','重叠50%','Location','NorthEast');
axis([0 N-1 -1.5 1.5]);
title('真实信号与接收信号的误差');
%分析由于频域的阶段效应带来的误差???有问题
评论0