Python信号处理:波束形成及目标方位估计,CBF、MVDR

2023-11-09 17:20

本文主要是介绍Python信号处理:波束形成及目标方位估计,CBF、MVDR,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

摘要:一直以来都是用MATLAB做信号处理,得到预处理的特征后再用Python进一步应用神经网络之类的方法。这里利用Python实现了常规波束形成(CBF)、MVDR波束形成以及波束扫描方位估计。本文实现的都是窄带波束形成,一种简单的宽带波束形成方法时直接将不同频率的窄带输出相加,比较容易扩展。

系列目录

  1. Python信号处理:快速傅里叶变换(FFT),短时傅里叶变换(STFT),窗函数,以及滤波
  2. Python信号处理:自相关函数(对标MATLAB中的autocorr)
  3. Python信号处理:波束形成及目标方位估计,CBF、MVDR
  4. Python信号处理:cvxpy工具包求解稀疏约束优化问题

目录

  1. 波束形成和方位估计
  2. 常规波束形成器
  3. MVDR波束形成器

1. 波束形成和方位估计

波束形成就是通过一个加权向量 w w w对阵列接收到的数据 x x x加权相加得到输出 y = w H x y = w^Hx y=wHx,所以确定 w w w是波束形成中非常重要的步骤,知道 w w w就知道了波束输出。波束形成器的设计问题其实就是设计加权向量 w w w

已知阵列的方向响应向量为 a ( θ ) = [ 1 , e − i 2 π λ d sin ⁡ θ , . . . , e − i 2 π λ ( N − 1 ) d sin ⁡ θ ] H a(\theta)=[1, e^{-\rm{i}\frac{2\pi}{\lambda}d\sin{\theta}}, ..., e^{-\rm{i}\frac{2\pi}{\lambda}(N-1)d\sin{\theta}}]^H a(θ)=[1,eiλ2πdsinθ,...,eiλ2π(N1)dsinθ]H
,则波束响应为
B ( θ ) = w H a ( θ ) . B(\theta)=w^Ha(\theta). B(θ)=wHa(θ).
波束响应只与加权向量和方向响应向量有关,而与是否存在信号无关。对于确定方向的加权向量 w w w,改变方向响应向量就可以得到波束响应能量相对于不同方向的函数图,也就是波束图。波束图反映了波束形成器的空域滤波性能。

加权向量 w w w是波束观察方向 θ 0 \theta_0 θ0的函数,因此在波束观察方向的波束输出功率可以表示为
σ y 2 ( θ 0 ) = y ( θ 0 ) y H ( θ 0 ) = w H ( θ 0 ) x x H w ( θ 0 ) = w ( θ 0 ) H R x w ( θ 0 ) , \sigma_y^2(\theta_0)=y(\theta_0)y^H(\theta_0)=w^H(\theta_0)xx^Hw(\theta_0)=w(\theta_0)^HR_xw(\theta_0), σy2(θ0)=y(θ0)yH(θ0)=wH(θ0)xxHw(θ0)=w(θ0)HRxw(θ0),
其中 R x R_x Rx是接收数据的协方差矩阵。改变波束观察方向 θ 0 \theta_0 θ0,就可以获得在整个观察区内的波束输出功率,此时的到的就是波束扫描方位谱

波束形成和方位估计,分别可以得到波束图和波束扫描方位谱,两者十分相似,在常规波束形成中,两者甚至在形式上一样的。但两者的意义差别很大。

  • 波束图用于表示波束形成器的空域滤波性能,即给定波束形成器,看它对各方位信号的响应大小。计算波束图时,波束加权向量是固定的,通过扫描阵列的方向响应向量获得不同方向的波束响应。从波束图可以观察到各方位信号对感兴趣的方位(波束观察方向)波束输出的影响。波束图与是否存在信号无关。

  • 方位谱用于估计各方位到达信号功率的大小。利用波束扫描法计算方位谱时,加权向量是观察方向的函数,方位谱的峰值所在可以用于估计信号的方位。

2. 常规波束形成器

常规波束形成器的原理是延时加权求和,加权向量正好就是阵列方向相应向量相位的变化,即 w = a ( θ 0 ) w=a(\theta_0) w=a(θ0),因此波束响应为
B C B F ( θ ) = a H ( θ 0 ) a ( θ ) , B_{CBF}(\theta)=a^H(\theta_0)a(\theta), BCBF(θ)=aH(θ0)a(θ),
其中,波束形成时的 θ 0 \theta_0 θ0是不变的。

用于方位估计时,信号方位 θ \theta θ时固定的,波束观察方向 θ 0 \theta_0 θ0在观察区中变化,得到方位谱为
σ y 2 ( θ 0 ) = y ( θ 0 ) y H ( θ 0 ) = w H ( θ 0 ) x x H w ( θ 0 ) = σ s 2 w H ( θ 0 ) a ( θ ) a H ( θ ) w ( θ 0 ) \begin{aligned} \sigma_y^2(\theta_0)&=y(\theta_0)y^H(\theta_0)=w^H(\theta_0)xx^Hw(\theta_0)\\ &=\sigma_s^2w^H(\theta_0)a(\theta)a^H(\theta)w(\theta_0)\\ \end{aligned} σy2(θ0)=y(θ0)yH(θ0)=wH(θ0)xxHw(θ0)=σs2wH(θ0)a(θ)aH(θ)w(θ0)

python实现如下,代码中变量命名和公式中相同,波束图和波束扫描方位谱如图所示。在白噪声背景下,CBF的两幅图是相同的。

import matplotlib.pyplot as plt
import numpy as npd = np.array([15]) * np.ones([18, 1])
f = 50
lmbda = 1500 / f
N = 18
n = np.arange(0, N, 1).reshape([-1, 1])# 波束形成
theta = np.arange(-90, 90, 0.1).reshape([1, -1])
theta_0 = 10
theta = theta / 180 * np.pi
theta_0 = theta_0 / 180 * np.pia = np.exp(1j * 2 * np.pi * np.sin(theta) * n * d / lmbda)
w = np.exp(1j * 2 * np.pi * np.sin(theta_0) * n * d / lmbda)
B = np.dot(a.transpose(), np.conj(w))
B = np.abs(B) / np.max(np.abs(B))fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(7, 2), dpi=300)
plt.subplots_adjust(wspace=0.1, hspace=0.05)ax[0].plot(theta.reshape([-1,1])*180/np.pi, 20 * np.log10(B.reshape([-1,1])), '--')# 方位估计
theta_0 = np.arange(-90, 90, 0.1).reshape([1, -1])
theta = 10
theta = theta / 180 * np.pi
theta_0 = theta_0 / 180 * np.pia = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta) * n * d / lmbda))
w = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta_0) * n * d / lmbda))sigma = np.transpose(np.conj(w)) * a * np.transpose(np.conj(a)) * w
sigma = np.diag(sigma)
sigma = sigma / np.max(sigma)ax[1].plot(theta_0.reshape([-1,1])*180/np.pi, 10 * np.log10(sigma))

在这里插入图片描述

3. MVDR波束形成器

最小方差无失真响应波束形成器,顾名思义有两个成分:1. 最小方差,2. 无失真,因此可以表示为
min ⁡ w w H R n w , s . t . w H a = 1 , \min_{w}\quad w^HR_nw, \quad \quad s.t.\quad w^Ha=1, wminwHRnw,s.t.wHa=1,
其中,优化目标表明为最小输出噪声方差,约束为无失真约束。实际中通常使用接收信号的协方差矩阵代替噪声协方差矩阵,此时方法也可以叫做最小功率无失真响应波束形成器(MPDR),但叫MVDR也可以。

通过求解该约束优化问题,就可以得到波束形成器的加权向量,进而得到波束输出响应。有两种思路求解该优化问题:

  1. 常规解法,通过lagrange乘子法,手动求解计算加权向量以及波束相应输出。
  2. 工具包求解,利用一些优化工具包,直接求解该优化问题得到答案——这种方法实际上解决了大部分的波束形成问题,因为很多波束形成问题都是一个目标函数,再加上很多的约束,如果工具包能求解对应的优化问题,波束输出就有了。

这里直接给出第一种方法求出来的加权向量。第二种解法在后面内容中会详细介绍,作为一种通用的波束形成方法。因为都是优化问题,波束形成在算法上几乎没有什么难题,只有工程中的问题。

scipy.optimize.minimize失败了,因为只能处理实数的问题。MATLAB的CVX工具箱是能处理复数的,不知道cvxpy怎么样。。

cvxpy好像解不出复数最优解,目标函数里可以有复数操作,但最终的值必须是实数的,而且解出来的最优解也一直是实数,肯定哪里不太对。。

w M V D R = R x − 1 a ( θ 0 ) a H ( θ 0 ) R x − 1 a ( θ 0 ) w_{MVDR}=\frac{R_x^{-1}a(\theta_0)}{a^H(\theta_0)R_x^{-1}a(\theta_0)} wMVDR=aH(θ0)Rx1a(θ0)Rx1a(θ0)
波束响应和波束输出功率,和第一部分的公式相同,即
B M V D R = w H a ( θ ) B_{MVDR}=w^Ha(\theta) BMVDR=wHa(θ)
σ M V D R ( θ 0 ) = w H ( θ 0 ) R x w ( θ 0 ) = 1 a H ( θ 0 ) R x − 1 a ( θ 0 ) \sigma_{MVDR}(\theta_0)=w^H(\theta_0)R_xw(\theta_0)=\frac{1}{a^H(\theta_0)R_x^{-1}a(\theta_0)} σMVDR(θ0)=wH(θ0)Rxw(θ0)=aH(θ0)Rx1a(θ0)1

python实现如下,代码中变量命名和公式中相同,波束图和波束扫描方位谱如图所示。从波束图可以看出,当信号方向(-10°)和波束观察方向(10°)不同时,MVDR在信号方向会形成一个凹陷,降低了信号(此时信号算是干扰)的能量。

import matplotlib.pyplot as plt
import numpy as npd = np.array([15]) * np.ones([18, 1])
f = 50
lmbda = 1500 / f
N = 18
n = np.arange(0, N, 1).reshape([-1, 1])# 信号
theta = -10
theta = theta / 180 * np.pi
x = np.exp(1j * 2 * np.pi * np.sin(theta) * n * d / lmbda)
Rx = np.mat(1000 * np.dot(x, np.transpose(np.conj(x))) + 1 * np.eye(N))# # 波束形成
theta = np.arange(-90, 90, 0.1).reshape([1, -1])
theta_0 = 10
theta = theta / 180 * np.pi
theta_0 = theta_0 / 180 * np.pia = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta) * n * d / lmbda))a_theta_0 = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta_0) * n * d / lmbda))
w = Rx.I * a_theta_0 / (a_theta_0.H * Rx.I * a_theta_0)B = w.H * a
B = np.abs(B) / np.max(np.abs(B))fig, ax = plt.subplots(1, 2, sharex=True, sharey=False, figsize=(7, 2), dpi=300)
plt.subplots_adjust(wspace=0.35, hspace=0.05)ax[0].plot(theta.reshape([-1,1])*180/np.pi, 20 * np.log10(B.reshape([-1,1])), '--')# 方位估计
theta_0 = np.arange(-90, 90, 0.1).reshape([1, -1])
theta_0 = theta_0 / 180 * np.pisigma = []
for i in range(theta_0.shape[1]):a_theta_0 = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta_0[:, i]) * n * d / lmbda))sigma.append((1 / (a_theta_0.H * Rx.I * a_theta_0)))sigma = np.array(sigma).reshape([-1, 1])ax[1].plot(theta_0.reshape([-1,1])*180/np.pi, 10 * np.log10(sigma))

在这里插入图片描述

这篇关于Python信号处理:波束形成及目标方位估计,CBF、MVDR的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/377516

相关文章

Python基于wxPython和FFmpeg开发一个视频标签工具

《Python基于wxPython和FFmpeg开发一个视频标签工具》在当今数字媒体时代,视频内容的管理和标记变得越来越重要,无论是研究人员需要对实验视频进行时间点标记,还是个人用户希望对家庭视频进行... 目录引言1. 应用概述2. 技术栈分析2.1 核心库和模块2.2 wxpython作为GUI选择的优

Python如何使用__slots__实现节省内存和性能优化

《Python如何使用__slots__实现节省内存和性能优化》你有想过,一个小小的__slots__能让你的Python类内存消耗直接减半吗,没错,今天咱们要聊的就是这个让人眼前一亮的技巧,感兴趣的... 目录背景:内存吃得满满的类__slots__:你的内存管理小助手举个大概的例子:看看效果如何?1.

Python+PyQt5实现多屏幕协同播放功能

《Python+PyQt5实现多屏幕协同播放功能》在现代会议展示、数字广告、展览展示等场景中,多屏幕协同播放已成为刚需,下面我们就来看看如何利用Python和PyQt5开发一套功能强大的跨屏播控系统吧... 目录一、项目概述:突破传统播放限制二、核心技术解析2.1 多屏管理机制2.2 播放引擎设计2.3 专

Python中随机休眠技术原理与应用详解

《Python中随机休眠技术原理与应用详解》在编程中,让程序暂停执行特定时间是常见需求,当需要引入不确定性时,随机休眠就成为关键技巧,下面我们就来看看Python中随机休眠技术的具体实现与应用吧... 目录引言一、实现原理与基础方法1.1 核心函数解析1.2 基础实现模板1.3 整数版实现二、典型应用场景2

Python实现无痛修改第三方库源码的方法详解

《Python实现无痛修改第三方库源码的方法详解》很多时候,我们下载的第三方库是不会有需求不满足的情况,但也有极少的情况,第三方库没有兼顾到需求,本文将介绍几个修改源码的操作,大家可以根据需求进行选择... 目录需求不符合模拟示例 1. 修改源文件2. 继承修改3. 猴子补丁4. 追踪局部变量需求不符合很

python+opencv处理颜色之将目标颜色转换实例代码

《python+opencv处理颜色之将目标颜色转换实例代码》OpenCV是一个的跨平台计算机视觉库,可以运行在Linux、Windows和MacOS操作系统上,:本文主要介绍python+ope... 目录下面是代码+ 效果 + 解释转HSV: 关于颜色总是要转HSV的掩膜再标注总结 目标:将红色的部分滤

Python 中的异步与同步深度解析(实践记录)

《Python中的异步与同步深度解析(实践记录)》在Python编程世界里,异步和同步的概念是理解程序执行流程和性能优化的关键,这篇文章将带你深入了解它们的差异,以及阻塞和非阻塞的特性,同时通过实际... 目录python中的异步与同步:深度解析与实践异步与同步的定义异步同步阻塞与非阻塞的概念阻塞非阻塞同步

Python Dash框架在数据可视化仪表板中的应用与实践记录

《PythonDash框架在数据可视化仪表板中的应用与实践记录》Python的PlotlyDash库提供了一种简便且强大的方式来构建和展示互动式数据仪表板,本篇文章将深入探讨如何使用Dash设计一... 目录python Dash框架在数据可视化仪表板中的应用与实践1. 什么是Plotly Dash?1.1

在C#中调用Python代码的两种实现方式

《在C#中调用Python代码的两种实现方式》:本文主要介绍在C#中调用Python代码的两种实现方式,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录C#调用python代码的方式1. 使用 Python.NET2. 使用外部进程调用 Python 脚本总结C#调

Python下载Pandas包的步骤

《Python下载Pandas包的步骤》:本文主要介绍Python下载Pandas包的步骤,在python中安装pandas库,我采取的方法是用PIP的方法在Python目标位置进行安装,本文给大... 目录安装步骤1、首先找到我们安装python的目录2、使用命令行到Python安装目录下3、我们回到Py