clc
clear all
close all
%% 常数定义
M=15; %阵元数
snapshots=1e5; %快拍数（时间点）时域采样点数
fc=40e6; %中心频率
snrdB=0;snr=10^(snrdB/20); %信噪比
c=3*10^8; %声速
fc1=60e6;
T = 1e-3; %信号脉宽
fs = 1e8; %采样频率
ts=1/fs; %采样周期
%t = -T/2:1/fs:T/2-1/fs; %时间
t=1/fs:1/fs:1/fs*snapshots;
B = 10e6; %带宽
K = B/T;
ff= fc1+1*K*t; %信号频率
fmax=max(ff); %lfm信号最高频率
fmin=min(ff); %lfm信号最低频率
lambda=c/fmax; %波长
d=0.5*lambda; %阵元间距
angle=[40 140]; %角度
P=length(angle); %声源数
tau=d/c*sin(angle*pi/180); % 相邻阵元间延迟时间，为-0.002，即-T
NFFT=2e5; %快速傅里叶变换点数
n=0:(NFFT-1);
s=zeros(P,snapshots);
s(1,:)=exp(1j*2*pi*(fc*t+0.5*K*t.^2)); %lfm信号1 虚信号
s(2,:)=exp(1j*2*pi*(fc*t-0.5*K*t.^2)); %lfm信号2
%% 阵元上的接收信号波形
ss1 = zeros(M,snapshots); %信号1
for j = 1:M % 生成15个阵元上的接收信号波形
ss1(j,:) = s(1,:).*exp(1j*2*pi*fc*(j-1)*d/c*sin(angle(1)*pi/180));
end
ss2 = zeros(M,snapshots); %信号2
for j = 1:M % 生成15个阵元上的接收信号波形
ss2(j,:) = s(2,:).*exp(1j*2*pi*fc*(j-1)*d/c*sin(angle(2)*pi/180));
end
%% 对阵元接收波形做fft变换
for i=1:M
ss1_f1(i,:) = fft(ss1(i,:),NFFT);
end
for i=1:M
ss2_f1(i,:)=fft(ss2(i,:),NFFT);
end
%% 对三段数据各子带进行波束形成
bf1 = zeros(1,NFFT);
bf2 = zeros(1,NFFT);
for i=floor(35e6/fs*NFFT):ceil(45e6/fs*NFFT) % 取100个频点，确保覆盖35e6-45e6HZ范围
f = i*fs/NFFT; % 每个频点处频率
w = exp(j*2*pi*f*d/c*sin(angle(1)*pi/180)*[0:M-1]')./M; % 常规波束形成，加权向量即为每个频点处阵列响应向量
bf1(i) = w'* ss1_f1(:,i);
end
for i=floor(35e6/fs*NFFT):ceil(45e6/fs*NFFT)
f = i*fs/NFFT; % 每个频点处频率
w = exp(j*2*pi*f*d/c*sin(angle(2))*[0:M-1]')./M; % 常规波束形成，加权向量即为每个频点处阵列响应向量
bf2(i) = w'* ss2_f1(:,i);
end
bf1_out = real(ifft(bf1,1e5));
bf2_out = real(ifft(bf2,1e5));
sr(1,:)=bf1_out;
sr(2,:)=bf2_out;
figure(1)
subplot(211)
plot(1:length(t),sr(1,:),'r'); % LFM信号时域波形
subplot(212)
plot(1:length(t),sr(2,:),'b');
xlim([1 1e5]);
xlabel('(fs = 1e8Hz,1e5个样点)');
ylabel('');
title('LFM信号时域波形')
%% CBF
for kk=1:M
signal(kk,:) = sr(1,:).*exp(1j*2*pi./(c./(fc+K*t))*d.*(kk-1)*cos(angle(1)*pi/180))+sr(2,:).*exp(1j*2*pi./(c./(fc-K*t))*d.*(kk-1)*cos(angle(2)*pi/180));
end;
n=exp(j*2*pi*randn(M,snapshots))/snr;%噪声
data=signal;
SS=data+n; %阵列接收数据
Nseg=1000; %时域分为1000段
duan=100; %每段100个时间点
NFFT1=2000; %做2000个点的快速傅里叶变换
SS_FFT = zeros(M,Nseg*NFFT1);
for seg=1:Nseg %数据段
for m = 1:M
SS_FFT(m,NFFT1*(seg-1)+1:NFFT1*seg)=fft(signal(m,duan*(seg-1)+1:duan*seg),NFFT1);%时域分成1000段 每段100个点 对每段的100个点做2000个点的fft变换 产生2000个频率数据 存入SS_FFT中
end
end
%%
P_temp=zeros(1,180);
Pmusic=zeros(1,180);
ii=0;
theta=1:1:180;
for fn=700:900
f=fs/NFFT1*fn;
fXX=SS_FFT(:,fn:NFFT1:(Nseg-1)*NFFT1+fn); %对每个频点建立类似时域的阵列接收数据
RfXX=fXX*fXX'/size(fXX,2); %频域的协方差阵
[EigenVectors,EigenValues]=eig(RfXX);
Lemda=diag(EigenValues); %计算矩阵特征值
[SortedLemda,Index]=sort(Lemda); %特征值升序排列
Index=flipud(Index); %将特征值降序排列
NoiseSubspace(:,1:M-P)=EigenVectors(:,Index(P+1:M));%把噪声空间分成两部分，分别为加入信号空间的一部分和噪声
for i=1:1:180
a=exp(1j*2*pi*f*d*(0:M-1)'*cos(i*pi/180)/c);
MUSIC_Spec(i)=a'*a/(a'*NoiseSubspace*(NoiseSubspace)'*a);
end
%P=1./sum(MUSIC_Spec);%每列相加
P_temp = 10*log10(MUSIC_Spec);
Pmusic = Pmusic+P_temp;
ii=ii+1;
end
Pmusic=Pmusic/ii;
figure(2)
plot(theta,(Pmusic))
grid on;
num=zeros(1,100);
loc=zeros(1,100);
[num loc]=findpeaks(abs(Pmusic));%查找峰值 num为峰值 loc为峰值所在位置
[a_num a_loc]=max(num); %查找峰值中的最大值 a_loc最大值所在位置
location_1=loc(a_loc);%提取最大峰值所在位置
num(a_loc)=min(abs(Pmusic));%在找出全部峰值后，把最大峰值赋成最小值
num_deal_max=num;
[b_num b_loc]=max(num);%找到剩下峰值的最大值
location_2=loc(b_loc);%次最大峰值的位置
angle=[theta(location_1) theta(location_2)]
如何用波束形成法做宽带信号DOA估计
RfXX是无穷矩阵
如何解决
