weixin_51363190
唯一1202
2021-04-01 16:06
采纳率: 0%
浏览 64

波束形成法做宽带doa

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是无穷矩阵

如何解决

  • 点赞
  • 写回答
  • 关注问题
  • 收藏
  • 邀请回答

2条回答 默认 最新

  • QA_Assistant
    有问必答小助手 2021-04-02 10:23

    你好,我是有问必答小助手。为了技术专家团更好地为您解答问题,烦请您补充下(1)问题背景详情,(2)您想解决的具体问题,(3)问题相关图片。便于技术专家团更好地理解问题,并给出解决方案。

    您可以点击问题下方的【编辑】,进行补充修改问题。

    点赞 评论
  • superdont
    superdont 2021-04-02 13:58

    这个是MATLAB代码,应该在MATLAB的模块内进行提问。 建议把代码分成模块来写,不要一气呵成写下来。不然排错还是相当困难的,因为很难看清楚你的思路到底是什么。

    点赞 评论