宽带非相干信号DOA估计仿真

119 阅读 0 评论 79 点赞

我是靠谱客的博主朴实服饰，这篇文章主要介绍宽带非相干信号DOA估计仿真，现在分享给大家，希望可以做个参考。

宽带非相干信号DOA估计仿真
产生两个相同的宽带LFM信号，采用非相干信号子空间法进行DOA估计。

复制代码

clear
 
Bsignal = 20e6;                       %调频带宽
Tsignal = 10e-6;                      %脉冲宽度
Ksignal = Bsignal/Tsignal;                        %调频斜率
Fssignal = 4*Bsignal;                       %采样频率
F0signal = 1e6;                       %发射信号时的瞬时频率，也是信号有效区间发射信号的中心频率
Tssignal = 1/Fssignal;                      %采样周期
Nsignal = ceil(Tsignal/Tssignal);                 %采样点数
tsignal = linspace(-Tsignal/2,Tsignal/2,Nsignal);

N=8;                                     %接收阵元个数
L=Nsignal;                                   %快拍数目
K=2;                                     %
c=3e8;                                   %传播速度
lamda=c/(F0signal+Bsignal/2);                               %波长
d=lamda/2;                            %发射阵元间距取半个波长 
theat0=[-10 20];                          %目标与阵列的夹角

%.................................................................................
%产生线性调频信号
LFM1 = ones(1,L).*exp(j *(2 * pi * F0signal * tsignal + pi * Ksignal * tsignal.^2));     %线性调频信号复数表达式
LFM2 = ones(1,L).*exp(j *(2 * pi * F0signal * tsignal + pi * Ksignal * tsignal.^2));     %线性调频信号复数表达式
Q=length(theat0);

LFM_s(1,:)=LFM1;
LFM_s(2,:)=LFM2;

for lamdaindex=1:10
    lfm_dived=LFM_s(1:2,1+40*(lamdaindex-1):40*lamdaindex);
    LFMfft(1,:)=fftshift(fft(lfm_dived(1,:)));
    [maxx maxxindex]=max(LFMfft(1,:));
    LFMfft(2,:)=fftshift(fft(lfm_dived(2,:)));
    lamda=c/((maxxindex-20)*1e6);

for q=1:Q
        A_shiliang(:,q)=exp(-1j*2*pi*d/lamda*sin(theat0(q)*pi/180)*[0:N-1]');%接收阵列引导矢量
    end
        rsignaldoa=LFMfft;

rsignal=A_shiliang*rsignaldoa;
    % 计算协方差矩阵
    Rxx=rsignal*rsignal'/40;
    % 特征值分解
    [EV,D]=eig(Rxx);                   %特征值分解
    EVA=diag(D)';                      %将特征值矩阵对角线提取并转为一行
    [EVA,I]=sort(EVA);                 %将特征值排序 从小到大
    EV=fliplr(EV(:,I));                % 对应特征矢量排序
    % 遍历每个角度，计算空间谱
    i=1;
    for iang = -90:0.1:90;
        phim=pi/180*iang;
        AMusic=exp(-1j*2*pi/lamda*d*sin(phim)*[0:N-1]).';    
        Pbf(i)=AMusic'*Rxx*AMusic;
        En=EV(:,K+1:N);                   % 取矩阵的第M+1到N列组成噪声子空间
        Pmusicj(i)=1/(AMusic'*En*En'*AMusic);
        i=i+1;
    end
    Pmusicj=abs(Pmusicj);
    Pmmaxj=max(Pmusicj);
    Pmusic(lamdaindex,:)=Pmusicj/Pmmaxj;            % 归一化处理
end

iang = -90:0.1:90;
figure
plot(iang,mean(Pmusic,1))
hold on

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
clear
 
Bsignal = 20e6;                       %调频带宽
Tsignal = 10e-6;                      %脉冲宽度
Ksignal = Bsignal/Tsignal;                        %调频斜率
Fssignal = 4*Bsignal;                       %采样频率
F0signal = 1e6;                       %发射信号时的瞬时频率，也是信号有效区间发射信号的中心频率
Tssignal = 1/Fssignal;                      %采样周期
Nsignal = ceil(Tsignal/Tssignal);                 %采样点数
tsignal = linspace(-Tsignal/2,Tsignal/2,Nsignal);

N=8;                                     %接收阵元个数
L=Nsignal;                                   %快拍数目
K=2;                                     %
c=3e8;                                   %传播速度
lamda=c/(F0signal+Bsignal/2);                               %波长
d=lamda/2;                            %发射阵元间距取半个波长 
theat0=[-10 20];                          %目标与阵列的夹角


%.................................................................................
%产生线性调频信号
LFM1 = ones(1,L).*exp(j *(2 * pi * F0signal * tsignal + pi * Ksignal * tsignal.^2));     %线性调频信号复数表达式
LFM2 = ones(1,L).*exp(j *(2 * pi * F0signal * tsignal + pi * Ksignal * tsignal.^2));     %线性调频信号复数表达式
Q=length(theat0);

LFM_s(1,:)=LFM1;
LFM_s(2,:)=LFM2;

for lamdaindex=1:10
    lfm_dived=LFM_s(1:2,1+40*(lamdaindex-1):40*lamdaindex);
    LFMfft(1,:)=fftshift(fft(lfm_dived(1,:)));
    [maxx maxxindex]=max(LFMfft(1,:));
    LFMfft(2,:)=fftshift(fft(lfm_dived(2,:)));
    lamda=c/((maxxindex-20)*1e6);

    for q=1:Q
        A_shiliang(:,q)=exp(-1j*2*pi*d/lamda*sin(theat0(q)*pi/180)*[0:N-1]');%接收阵列引导矢量
    end
        rsignaldoa=LFMfft;

    rsignal=A_shiliang*rsignaldoa;
    % 计算协方差矩阵
    Rxx=rsignal*rsignal'/40;
    % 特征值分解
    [EV,D]=eig(Rxx);                   %特征值分解
    EVA=diag(D)';                      %将特征值矩阵对角线提取并转为一行
    [EVA,I]=sort(EVA);                 %将特征值排序 从小到大
    EV=fliplr(EV(:,I));                % 对应特征矢量排序
    % 遍历每个角度，计算空间谱
    i=1;
    for iang = -90:0.1:90;
        phim=pi/180*iang;
        AMusic=exp(-1j*2*pi/lamda*d*sin(phim)*[0:N-1]).';    
        Pbf(i)=AMusic'*Rxx*AMusic;
        En=EV(:,K+1:N);                   % 取矩阵的第M+1到N列组成噪声子空间
        Pmusicj(i)=1/(AMusic'*En*En'*AMusic);
        i=i+1;
    end
    Pmusicj=abs(Pmusicj);
    Pmmaxj=max(Pmusicj);
    Pmusic(lamdaindex,:)=Pmusicj/Pmmaxj;            % 归一化处理
end

iang = -90:0.1:90;
figure
plot(iang,mean(Pmusic,1))
hold on