关于阵列信号处理的一些东西（一）

于 2023-11-17 17:10:06 发布
阅读量276 收藏 5
点赞数 2
分类专栏： 雷达信号处理 文章标签： 信号处理 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_45946324/article/details/134439950
版权

空间匹配滤波

阵列方向图函数
F ( θ ) = W H ∗ α ( θ ) F(\theta) = W^{H}*\alpha(\theta) F(θ)=WH∗α(θ)
当权重矢量 W = α ( θ 0 ) W=\alpha(\theta_0) W=α(θ0​)时，阵列方向图增益G(θ)在θ0处增益达到最大。

f0 = 16e9;  % 雷达工作频率16GHz
c = 3e8;
lambda = c/f0;
d = lambda/2;   % 阵元距离
sita = -5/9*pi:0.01:5/9*pi;
sita0 = pi/6;
n = 16;
N = 0:16-1;     % 16个阵元
move_sita = exp(1j*N'*2*pi*d*sin(sita0)/lambda);    % 权重
alpha_sita = exp(1j*N'*2*pi*d*sin(sita)/lambda);    % 导向矢量
F = move_sita'*alpha_sita;
refer = sin(pi*d*n*(sin(sita))/lambda)/n./sin(pi*d*(sin(sita))/lambda);
figure('Name','空域匹配滤波');
plot(sita/2/pi*360,20*log10(abs(F)/max(abs(F))));
hold on
plot(sita/2/pi*360,20*log10(abs(refer)/max(abs(refer))),'--');
xticks(-100:20:100);
xlabel('角度');
ylabel('归一化功率增益');
hold off

在30°获得最大增益

最优波束形成

最大信干噪比(MINR)准则

在这里插入图片描述
在这里插入图片描述

%% MSINR准则下使用带功率的协方差公式
clc;clear;close;
N = 16;
n = 0:N-1;
f0 = 10e9;  % 雷达工作频率16GHz
c = 3e8;
lambda = c/f0;
d = lambda/2;   % 阵元距离
sita0 = 0;      % 期望信号方向
signal_e = 10^5;    % 信号能量
sita1 = deg2rad(-35);
sita2 = deg2rad(-19);
sita3 = deg2rad(45);
sita = [sita1 sita2 sita3];     % 干扰方向
INR1 = 10^(40/10);
INR2 = 10^(35/10);
INR3 = 10^(50/10);
INR = [INR1 INR2 INR3];         % 干扰强度
A = exp(1j*2*pi*n'*sin(sita)*d/lambda);     % 干扰信号的导向矢量
alpha0 = exp(1j*2*pi*n'*sin(sita0)*d/lambda);   % 期望信号的导向矢量值
Rin = zeros(N,N);
for i = 1:length(sita)
    Rin = Rin + INR(i)*A(:,i)*A(:,i)';          % 干噪协方差矩阵
end
Rin = Rin + eye(N);                         % 加入噪声
RS = signal_e*(alpha0*alpha0');
[V,D] = eig(RS,Rin);        % 得到广义特征值与特征向量
[~,I] = max(diag(D));
w = V(:,I);              % 最大特征值所对应的特征向量就是最优权矢量
angle = -5/9*pi:0.01:5/9*pi;
A1 = exp(1j*2*pi*n'*sin(angle)*d/lambda);
CBF = db(abs(sum(A1)).^2/N/N)/2;          % 权重全为1时的增益方向图
MSINR = abs(w'*A1);
MSINR = db(MSINR/max(MSINR));

CBF(CBF < -80) = -80;
MSINR(MSINR < -80) = -80;

[~,ind1] = min(abs(angle-sita1));
[~,ind2] = min(abs(angle-sita2));
[~,ind3] = min(abs(angle-sita3));

figure('Name','最优化波束形成');
plot(rad2deg(angle),CBF);
hold on
plot(rad2deg(angle),MSINR,'--');
plot(rad2deg(sita1), MSINR(ind1),'ro', 'MarkerFaceColor', 'r');
text(rad2deg(sita1), MSINR(ind1), '\leftarrow干扰1');
plot(rad2deg(sita2), MSINR(ind2),'ro', 'MarkerFaceColor', 'r');
text(rad2deg(sita2), MSINR(ind2), '\leftarrow干扰2');
plot(rad2deg(sita3), MSINR(ind3),'ro', 'MarkerFaceColor', 'r');
text(rad2deg(sita3), MSINR(ind3), '\leftarrow干扰3');
xlabel('角度');
ylabel('归一化功率增益');
legend('CBF','MSINR');
hold off

在这里插入图片描述
MSINR的图像相比于匹配滤波在干扰处确实有塌陷，但没有书上给的图效果好。

%% MSINR准则下使用快拍的方式求协方差矩阵
clc;clear;close;
N = 16;
n = 0:N-1;
f0 = 16e9;  % 雷达工作频率16GHz
c = 3e8;
lambda = c/f0;
d = lambda/2;   % 阵元距离
INR = [40 35 50];
SNR = 40;
theta = [-40 -15 35];
theta0 = 20;
snap = 1024;        % 快拍数
signal = exp(1j*2*pi*f0*(0:snap-1)/(2*snap));
INR = 10.^(INR/10)/2;
A = exp(1j*2*pi*n'*sin(deg2rad(theta))*d/lambda);
alpha0 = exp(1j*2*pi*n'*sin(deg2rad(theta0))*d/lambda);
xs = sqrt(10^(SNR/10))*alpha0*signal;
kj = length(theta);     % 干扰数
rs = zeros(kj,snap);
for i = 1:kj
    for j = 1:snap
        rs(i,j) = sqrt(INR(i))*(randn(1)+1j*randn(1))*signal(j);
    end
end
noise = (randn(N,snap)+1j*randn(N,snap))/sqrt(2);
Rs = 1/snap*(xs*xs');
J = A*rs;
Rin = 1/snap*(J*J'+noise*noise');
[V,D] = eig(Rs,Rin);
[~,I] = max(diag(D));
w = V(:,I);
w0 = exp(1j*2*pi*n'*sin(deg2rad(theta0))*d/lambda);
angle = -5/9*pi:0.001:5/9*pi;
F = w'*exp(1j*2*pi*n'*sin(angle)*d/lambda);
F0 = w0'*exp(1j*2*pi*n'*sin(angle)*d/lambda);
figure('Name','MSINR');
F1 = db(abs(F)/max(abs(F)));
plot(rad2deg(angle),F1);
hold on
plot(rad2deg(angle),db(abs(F0)/max(abs(F0))),'--');

angle = rad2deg(angle);
[~,ind1] = min(abs(angle-theta(1)));
[~,ind2] = min(abs(angle-theta(2)));
[~,ind3] = min(abs(angle-theta(3)));

text(theta(1),F1(ind1),'\leftarrow干扰1');
plot(theta(1),F1(ind1),'ro','MarkerFaceColor','r');
text(theta(2),F1(ind2),'\leftarrow干扰2');
plot(theta(2),F1(ind2),'ro','MarkerFaceColor','r');
text(theta(3),F1(ind3),'\leftarrow干扰3');
plot(theta(3),F1(ind3),'ro','MarkerFaceColor','r');

legend('最优化','匹配滤波');
title('MSINR,干扰方向[-40 -15 35]');
hold off

在这里插入图片描述
效果明显更好，但目前不知道这两种方式为什么会有这么大的差异。

线性约束最小方差准则(LCMV)与最小方差无失真响应(MVDR)

在这里插入图片描述
在这里插入图片描述

%% MVDR
clc;clear;close;
N = 16;
n = 0:N-1;
f0 = 16e9;  % 雷达工作频率16GHz
c = 3e8;
lambda = c/f0;
d = lambda/2;   % 阵元距离
INR = [40 35 50];
SNR = 40;
theta = [-40 -15 35];
theta0 = 20;
snap = 1024;        % 快拍数
signal = exp(1j*2*pi*f0*(0:snap-1)/(2*snap));
INR = 10.^(INR/10)/2;
A = exp(1j*2*pi*n'*sin(deg2rad(theta))*d/lambda);
alpha0 = exp(1j*2*pi*n'*sin(deg2rad(theta0))*d/lambda);
xs = sqrt(10^(SNR/10))*alpha0*signal;
kj = length(theta);     % 干扰数
rs = zeros(kj,snap);
for i = 1:kj
    for j = 1:snap
        rs(i,j) = sqrt(INR(i))*(randn(1)+1j*randn(1))*signal(j);
    end
end
noise = (randn(N,snap)+1j*randn(N,snap))/sqrt(2);
X = A*rs + noise;          
Rx_pinv = pinv(1/snap*(X*X'));
% J = A*rs;
% Rx_pinv = 1/snap*pinv(J*J'+noise*noise');
w = Rx_pinv*alpha0/(alpha0'*Rx_pinv*alpha0);
w0 = exp(1j*2*pi*n'*sin(deg2rad(theta0))*d/lambda);
angle = -5/9*pi:0.001:5/9*pi;
F = w'*exp(1j*2*pi*n'*sin(angle)*d/lambda);
F0 = w0'*exp(1j*2*pi*n'*sin(angle)*d/lambda);
figure('Name','MVDR');
F1 = db(abs(F)/max(abs(F)));
plot(rad2deg(angle),F1);
hold on
plot(rad2deg(angle),db(abs(F0)/max(abs(F0))),'--');

angle = rad2deg(angle);
[~,ind1] = min(abs(angle-theta(1)));
[~,ind2] = min(abs(angle-theta(2)));
[~,ind3] = min(abs(angle-theta(3)));

text(theta(1),F1(ind1),'\leftarrow干扰1');
plot(theta(1),F1(ind1),'ro','MarkerFaceColor','r');
text(theta(2),F1(ind2),'\leftarrow干扰2');
plot(theta(2),F1(ind2),'ro','MarkerFaceColor','r');
text(theta(3),F1(ind3),'\leftarrow干扰3');
plot(theta(3),F1(ind3),'ro','MarkerFaceColor','r');

legend('最优化','匹配滤波');
title('MVDR,干扰方向[-40 -15 35]');
hold off

在这里插入图片描述
比较MSINR和MVDR的图像可以推断两者在Ri+n和θ0已知并相等的情况下确实是等价的。
在这里插入图片描述

%% LCMV
clc;clear;close;
N = 16;
n = 0:N-1;
f0 = 16e9;  % 雷达工作频率16GHz
c = 3e8;
lambda = c/f0;
d = lambda/2;   % 阵元距离
INR = [40 35 50];
SNR = 40;
theta = [-40 -15 35];
theta0 = 20;
theta1 = -10;
snap = 1024;        % 快拍数
signal = exp(1j*2*pi*f0*(0:snap-1)/(2*snap));
INR = 10.^(INR/10)/2;
A = exp(1j*2*pi*n'*sin(deg2rad(theta))*d/lambda);
alpha0 = exp(1j*2*pi*n'*sin(deg2rad(theta0))*d/lambda);
alpha1 = exp(1j*2*pi*n'*sin(deg2rad(theta1))*d/lambda);
C = [alpha0 alpha1];			% 约束矩阵
[~,L] = size(C);
f = ones(L,1);
kj = length(theta);     % 干扰数
rs = zeros(kj,snap);
for i = 1:kj
    for j = 1:snap
        rs(i,j) = sqrt(INR(i))*(randn(1)+1j*randn(1))*signal(j);
    end
end
noise = (randn(N,snap)+1j*randn(N,snap))/sqrt(2);
X = A*rs + noise;          
% X = A*rs + noise + C*[signal;signal]; 	%将Rin替换成Rx
Rx_pinv = pinv(1/snap*(X*X'));
w = Rx_pinv*C/(C'*Rx_pinv*C)*f;
angle = -5/9*pi:0.001:5/9*pi;
F = w'*exp(1j*2*pi*n'*sin(angle)*d/lambda);
figure('Name','LCMV');
F1 = db(abs(F)/max(abs(F)));
plot(rad2deg(angle),F1);
hold on
angle = rad2deg(angle);
[~,ind1] = min(abs(angle-theta(1)));
[~,ind2] = min(abs(angle-theta(2)));
[~,ind3] = min(abs(angle-theta(3)));
[~,ind4] = min(abs(angle-theta0));
[~,ind5] = min(abs(angle-theta1));

text(theta(1),F1(ind1),'\leftarrow干扰1');
plot(theta(1),F1(ind1),'ro','MarkerFaceColor','r');
text(theta(2),F1(ind2),'\leftarrow干扰2');
plot(theta(2),F1(ind2),'ro','MarkerFaceColor','r');
text(theta(3),F1(ind3),'\leftarrow干扰3');
plot(theta(3),F1(ind3),'ro','MarkerFaceColor','r');
text(theta0,F1(ind4),'目标1');
plot(theta0,F1(ind4),'ro','MarkerFaceColor','r');
text(theta1,F1(ind5),'目标2');
plot(theta1,F1(ind5),'ro','MarkerFaceColor','r');
title('LCMV,干扰方向[-40 -15 35] 目标方向[-15 20]');
hold off

在这里插入图片描述

若将上述代码的 f 改为[1 0]'，则图像如下图所示。
在这里插入图片描述
将代码中的Ri+n替换成Rx，方向图旁瓣上升。在这里插入图片描述

均方误差最小准则(MSE)

在这里插入图片描述
在这里插入图片描述

三种准则比较

在这里插入图片描述

SMI(采样协方差矩阵求逆)

在这里插入图片描述
有限次快拍下的方向图副瓣较高，且快拍数越少，副瓣越高。

%% SMI
clc;clear;close;
N = 16;
n = 0:N-1; 
f0 = 16e9;  % 雷达工作频率16GHz
c = 3e8;
lambda = c/f0;
d = lambda/2;   % 阵元距离
INR = [40 35 50];
SNR = 40;
theta = [-40 -15 35];
theta0 = 20;
INR = 10.^(INR/10)/2;
A = exp(1j*2*pi*n'*sin(deg2rad(theta))*d/lambda);
alpha0 = exp(1j*2*pi*n'*sin(deg2rad(theta0))*d/lambda);
kj = length(theta);     % 干扰数
snaps = [32 128 1024];
times = 50;         % 统计次数
figure('Name','SMI');
hold on
for k = 1:length(snaps)
    sumF = 0;
    for t = 1:times
        snap = snaps(k);        % 快拍数
        signal = exp(1j*2*pi*f0*(0:snap-1)/(2*snap));
        rs = zeros(kj,snap);
        for i = 1:kj
            for j = 1:snap
                rs(i,j) = sqrt(INR(i))*(randn(1)+1j*randn(1))*signal(j);
            end
        end
        noise = (randn(N,snap)+1j*randn(N,snap))/sqrt(2);
        X = A*rs + noise;          
        Rx_pinv = pinv(1/snap*(X*X'));
        w = Rx_pinv*alpha0/(alpha0'*Rx_pinv*alpha0);
        angle = -5/9*pi:0.001:5/9*pi;
        F = w'*exp(1j*2*pi*n'*sin(angle)*d/lambda);
        F1 = db(abs(F)/max(abs(F)));
        sumF = sumF + F1;
    end
    sumF = sumF/times;
    plot(rad2deg(angle),sumF);
end

angle = rad2deg(angle);
[~,ind1] = min(abs(angle-theta(1)));
[~,ind2] = min(abs(angle-theta(2)));
[~,ind3] = min(abs(angle-theta(3)));

text(theta(1),sumF(ind1),'\leftarrow干扰1');
plot(theta(1),sumF(ind1),'ro','MarkerFaceColor','r');
text(theta(2),sumF(ind2),'\leftarrow干扰2');
plot(theta(2),sumF(ind2),'ro','MarkerFaceColor','r');
text(theta(3),sumF(ind3),'\leftarrow干扰3');
plot(theta(3),sumF(ind3),'ro','MarkerFaceColor','r');
legend('snap = 32','snap = 128','snap = 1024');
title('SMI,干扰方向[-40 -15 35]');
hold off

在这里插入图片描述

对角加载技术

在这里插入图片描述

%% SMI
clc;clear;close;
N = 16;
n = 0:N-1; 
f0 = 16e9;  % 雷达工作频率16GHz
c = 3e8;
lambda = c/f0;
d = lambda/2;   % 阵元距离
DL = [0 10];        % 10db对角加载，这里是真值
INR = [40 35 50];
SNR = 40;
theta = [-40 -15 35];
theta0 = 20;
INR = 10.^(INR/10)/2;
A = exp(1j*2*pi*n'*sin(deg2rad(theta))*d/lambda);
alpha0 = exp(1j*2*pi*n'*sin(deg2rad(theta0))*d/lambda);
kj = length(theta);     % 干扰数
snaps = [32,1024];
times = 50;         % 统计次数
figure('Name','SMI');
hold on
sumF = 0;
for s = 1:length(snaps)
    snap = snaps(s);
    if s == 2
        DL = [0];
    end
    for k = 1:length(DL)
        for t = 1:times
            signal = exp(1j*2*pi*f0*(0:snap-1)/(2*snap));
            rs = zeros(kj,snap);
            for i = 1:kj
                for j = 1:snap
                    rs(i,j) = sqrt(INR(i))*(randn(1)+1j*randn(1))*signal(j);
                end
            end
            noise = (randn(N,snap)+1j*randn(N,snap))/sqrt(2);
            X = A*rs + noise;          
            Rx_pinv = pinv(1/snap*(X*X') + DL(k)*eye(N));			%对角加载
            w = Rx_pinv*alpha0/(alpha0'*Rx_pinv*alpha0);
            angle = -5/9*pi:0.001:5/9*pi;
            F = w'*exp(1j*2*pi*n'*sin(angle)*d/lambda);
            F1 = db(abs(F)/max(abs(F)));
            sumF = sumF + F1;
        end
        sumF = sumF/times;
        plot(rad2deg(angle),sumF);
    end
end
angle = rad2deg(angle);
[~,ind1] = min(abs(angle-theta(1)));
[~,ind2] = min(abs(angle-theta(2)));
[~,ind3] = min(abs(angle-theta(3)));

text(theta(1),sumF(ind1),'\leftarrow干扰1');
plot(theta(1),sumF(ind1),'ro','MarkerFaceColor','r');
text(theta(2),sumF(ind2),'\leftarrow干扰2');
plot(theta(2),sumF(ind2),'ro','MarkerFaceColor','r');
text(theta(3),sumF(ind3),'\leftarrow干扰3');
plot(theta(3),sumF(ind3),'ro','MarkerFaceColor','r');
legend('snap = 32','对角加载DL','snap = 1024');
title('SMI,干扰方向[-40 -15 35]');
hold off

在这里插入图片描述
在这里插入图片描述

广义旁瓣相消器(GSC)

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
最终得到计算公式如下所示
在这里插入图片描述

参考

线性最优离散滤波器——维纳滤波器及LCMV MVDR GSC （自适应滤波）

波束形成算法综述

b站视频

码乐乐
关注
专栏目录
阵列信号处理
m0_50086718的博客
07-23 4527
一、阵列信号模型 在阵列信号处理中，空间电磁环境、噪声特性、传播媒介等因素的变化，将 会引起信号模型的变化。为简化分析，做出以下假设。 信源：点源，位于远场区域。 噪声：加性零均值高斯白噪声。 传播媒介：各向同性、无色散，信号在媒介中传播时满足叠加原理。 其原理如图 1所示。 图 1阵列原理图 假设阵元为窄带信号，到达阵列的第n个信号则可以表示为 ...
阵列信号处理学习笔记（一）--阵列信号处理定义
最新发布
小豆包的博客
07-20 1489
MOOC 阵列信号处理 国防科技大学电子科学学院视频地址
阵列信号处理知识点合集
梦 醉了谁的博客
08-06 3836
无
（初级）数字信号处理目录（不只是目录）
Reborn Lee
07-29 3043
目录 前言 复指数信号与傅里叶分析 LTI系统对复指数信号的响应 连续周期时间信号的傅里叶级数（FS） 离散时间复指数信号的周期性质 离散周期信号的傅里叶级数（DFS） 傅里叶级数与线性时不变系统 滤波器简介 周期方波的傅里叶级数系数 连续时间傅里叶变换（FT） 周期信号的傅里叶变换 几个判断线性时不变系统的精彩例子 连续时间傅里叶变换的共轭与共轭对称性 连续时间傅里叶...
小波变换滤波器组在一维信号矢量的完美恢复
沙特国王大学学报小波变换滤波器组在一维信号矢量Fahima Tabassuma，Md.Imdadul Islama，M.R.阿明ba孟加拉国达卡贾汉吉尔纳加尔大学计算机科学与工程系b孟加拉国达卡东西大学电子和通信工程系阿提奇莱因福奥文章历史...
"2015年雷恩第一大学博士论文：高效CAD流程中的运行时硬件重新配置
2015年论文/雷恩第一大学由布列塔尼欧洲大学盖章为的等级雷恩第一大学博士提及：信号处理和电信马蒂斯博士学校提交人ChristopheHURI aU x在IRISA研究单位（UMR 6074）计算机科学和随机系统研究所编写国立应用科学...
动态场景视频去噪中的原始视频信号处理
2301基于Benchmark数据集的动态场景监督式原始视频去噪天津大学电气与信息工程学院{huanjing.yue，caocong 123，leolei，chu rh，yjy}@tju.edu.cnhttps://github.com/cao-cong/RViDeNet摘要近年来，有监督学习策略在...
YB菜菜的毫米波雷达自学之路——LFMCW雷达信号概述和信号模型初步建立
byang_18的博客
11-28 2949
雷达信号模型
阵列信号处理(全).ppt
12-12
阵列信号与波束形成的入门级资料，很不错。需要的可以下载了。 掌握空间传播波携带信号的获取与处理的基本理论和方法，特别是空时多维信号算法，熟悉参数估计和自适应波束形成的常用算法。
阵列信号处理试题及答案
03-15
中国科学院研究生课程阵列信号处理的期末考试题目及答案，仅供参考
第二章_优化阵列信号处理_
09-29
在学习过程中自行学习编写的《优化阵列信号处理》例题程序，有需要的下。
第三章(二)_优化阵列信号处理_
10-03
这是本人在学习过程中自行学习编写的《优化阵列信号处理》中第三章后几个例题的程序，有需要的下。
第三章(一)_优化阵列信号处理_
10-03
在学习过程中自行学习编写的《优化阵列信号处理》例题程序，有需要的下。
阵列信号处理（信源传播模型）
weixin_43195281的博客
07-26 2235
阵列信号处理的远场信源模型与近场信源模型介绍。
阵列信号处理（波束形成）
weixin_43195281的博客
11-15 1304
普通波束形成
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值