CN110333504B - 一种空时二维滤波的快速宽带波束形成方法 - Google Patents

一种空时二维滤波的快速宽带波束形成方法 Download PDF

Info

Publication number
CN110333504B
CN110333504B CN201910640904.3A CN201910640904A CN110333504B CN 110333504 B CN110333504 B CN 110333504B CN 201910640904 A CN201910640904 A CN 201910640904A CN 110333504 B CN110333504 B CN 110333504B
Authority
CN
China
Prior art keywords
integer
row
column
points
sampling
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910640904.3A
Other languages
English (en)
Other versions
CN110333504A (zh
Inventor
李秀坤
王集
徐天杨
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Engineering University
Original Assignee
Harbin Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN201910640904.3A priority Critical patent/CN110333504B/zh
Publication of CN110333504A publication Critical patent/CN110333504A/zh
Application granted granted Critical
Publication of CN110333504B publication Critical patent/CN110333504B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明属于水下声传感器阵列信号处理领域，具体涉及一种空时二维滤波的快速宽带波束形成方法，包括以下步骤：以符合Niquest定理的采样率进行采样，每个通道采集的数据按时间顺序排列成列向量，组成一个T×N维的数据矩阵，列的编号依次为j＝0,1,…,N‑1，行的编号依次为i＝0,1,…,T‑1；确定采样倍数M以及相邻通道时延补偿点数T0；计算坐标矩阵主值区间P(0)(T0)，规定坐标矩阵的列编号从左至右依次为j＝0,1,…,N‑1，行编号从上至下依次为i＝0,1,…,M‑1，则P(0)(T0)的第i行第j列元素为
Figure DDA0002131838330000011
本发明提出了坐标矩阵计算方法，并利用坐标矩阵直接对有效数据进行计算，略去了大部分的冗余计算量，从而达到快速计算的目的。

Description

一种空时二维滤波的快速宽带波束形成方法
技术领域
本发明属于水下声传感器阵列信号处理领域，具体涉及一种空时二维滤波的快速宽带波束形成方法。
背景技术
在水下被动定位中，目标声源往往为宽带声信号。快速、准确的宽带波束形成技术对提升水下武器装备的性能具有重要意义。
常规的宽带数字波束形成建立在量化时延的基础上，因而输出波束的扫描范围、角度步进值等参量均受到采样率的限制。为了保证输出波束的质量，往往需要较高的采样率。当接收阵为标准半波长均匀线阵时，
Figure BDA0002131838310000011
描述了宽带时域波束形成的过程。其中， N表示线阵中一共存在N个阵元，xi(n)表示第i个通道的接收信号时间序列。τs为时延补偿，该补偿必须为采样周期的整数倍。即τs＝kTs,k＝0,±1,±2,…，其中Ts为采样周期，与采样率fs的关系为
Figure BDA0002131838310000012
由于k只能取离散值，所以τs只能为离散值。若希望使用更小数值的τs进行精确度更高的扫描，则只能使用更小的Ts值，即需要提高系统的采样率。根据采样定理，为了能够不失真地恢复信号，采样率应当为信号最高频率的2倍，即萘奎斯特(Niquest)采样定理。在数字多波束中，为了保证增益不明显下降，通常要求采样率为信号最高频率的3-5倍。而这种要求只能保证增益下降不多，若希望波束指向性图失真很小，则一般要求采样率为信号最高频率的10倍以上，由此导致高采样率。高采样率对系统整体性能提出了更高的要求。在数据采集阶段，高采样率要求使用高转换速率的AD芯片，造成了硬件实现上的困难，提高了系统成本。在数据存储和传输阶段，高采样率导致了大的数据吞吐量，对数据存储设备、数据传输过程都提出了更高的要求。在数据处理阶段，更大的数据量要求处理器具有更强的计算能力。因此，需要寻求一种方法，既能够保证输出波束的质量，又能够保持较低的采样率。
采样定理表明，只要按照萘奎斯特采样定律对模拟信号进行采样，就可以无失真地将其恢复。所以，可以按照较低的采样率进行采样，在处理阶段提高采样率，这样便可以在使用低样率的同时保证输出波束的质量。通常的方法是，对每个通道采集的信号进行升采样操作，然后在高采样率下进行时延波束形成。对于存在N个通道的声呐接收阵，该方法需要进行N 次插值操作。考虑到插值是一种计算量较大的数值计算方法，所以这种方法虽然降低了对硬件设备的要求，但具有较大的计算量。利用时延波束形成和升采样均为线性过程这一特点，还可以先对接收信号进行补零，然后进行时延波束形成，最后再通过插值的方法获得高采样率下的波束。对于存在N个通道的接收阵，这种方法只需要进行1次插值操作，所以降低了计算量。但由于补零过程的存在，使得在波束形成阶段的计算量成倍地增加。当希望进行M 倍升采样时，计算量则提升M倍。这限制了宽带时延波束形成的快速实现。
发明内容
针对升采样导致计算量成倍地增加的问题，本发明提出了一种基于空时二维滤波的快速宽带时延波束形成方法。一方面，该方法能够得到任意原采样率整数倍数的采样率下的输出波束，从而能够保证输出波束的精度。另一方面，该方法在保证输出波束质量的前提下几乎没有增大计算量。实际上，对于传统方法而言，升采样M倍时计算量也提升了M倍。而本文所述方法为传统方法的
Figure BDA0002131838310000021
即本发明所述方法在不增加计算量的前提下完成了升采样M倍的任务，计算量与升采样倍数无关且保持恒定。
下面对本发明方法的原理进行阐述。
一种空时二维滤波的快速宽带波束形成方法，包括以下步骤：
(1)以符合Niquest定理的采样率进行采样，每个通道采集的数据按时间顺序排列成列向量，组成一个T×N维的数据矩阵，列的编号依次为j＝0,1,…,N-1，行的编号依次为i＝ 0,1,…,T-1；
(2)确定采样倍数M以及相邻通道时延补偿点数T0；
(3)计算坐标矩阵主值区间P(0)(T0)，规定坐标矩阵的列编号从左至右依次为 j＝0,1,…,N-1，行编号从上至下依次为i＝0,1,…,M-1，则P(0)(T0)的第i行第j列元素为
Figure BDA0002131838310000022
(4)去除P(0)(T0)中非整数点，只保留整数点；P(0)(T0)的每列有且只有1个整数点元素，设定循环因子k＝0；
(5)依次取出P(k)(T0)中第0～M-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束
Figure BDA0002131838310000023
的第kM,kM+1,…,kM+M-1个点；
(6)对P(k)(T0)中所有整数点数值加1，得到下一个周期的坐标矩阵P(k+1)(T0)；更新循环因子：k←k+1；
(7)判断kM≥T是否成立，若成立，则说明已经遍历数据矩阵中所有数值点，执行步骤(8)，否则执行步骤(5)；
(8)计算滤波器通带范围；
(9)根据步骤(8)中计算的参数设计滤波器对
Figure BDA0002131838310000031
进行滤波，得到输出波束y(i)。
所述计算坐标矩阵主值区间P(0)(T0)，规定坐标矩阵的列编号从左至右依次为 j＝0,1,…,N-1，行编号从上至下依次为i＝0,1,…,M-1，则P(0)(T0)的第i行第j列元素为
Figure BDA0002131838310000032
包括：
由步骤(2)确定的采样倍数M以及相邻通道时延补偿点数T0得输出波束表达式为：
Figure BDA0002131838310000033
其中，
Figure BDA0002131838310000034
表示第j个通道中序号为
Figure BDA0002131838310000035
的采样点参与了累加，得到了输出波束的第 i个数据点，
Figure BDA0002131838310000036
显然，
Figure BDA0002131838310000037
必为整数，由Niquest采样定理可知，只要采样率大于信号最高频率的2倍，则通过插值的方式恢复信号而无信息丢失；阵列信号中，目标信号最高频率
Figure BDA0002131838310000038
对接收数据进行M倍降采样处理，利用补零和时延求和得到波束
Figure BDA0002131838310000039
然后再通过低通滤波无失真地恢复输出波束，具体过程表示为：
Figure BDA00021318383100000310
其中，
Figure BDA00021318383100000311
Figure BDA00021318383100000312
为降采样M倍的信号，
Figure BDA00021318383100000313
为第j个通道参与波束形成的数据序号；
Figure BDA00021318383100000314
与pij(T0)关系为：
Figure BDA00021318383100000315
其中，非整数的pij(T0)表示序号为
Figure BDA00021318383100000316
的数据点不参与降采样后的累加过程，即在降采样后，这部分数据点不参与运算；
定义P(T0)为：
pij(T0)＝{P(T0)}ij,(i＝0,1,…,T-1,j＝0,1,…,N-1)
P(T0)为T×N的矩阵，称为坐标矩阵，将坐标矩阵的行从上至下依次称为第0行，第1 行，…，第T-1行，同时将坐标矩阵的列从左至右依次称为第0列，第1列，…，第N-1列，于是，P(T0)的第i行，第j列元素为pij(T0)。
所述去除P(0)(T0)中非整数点，只保留整数点；P(0)(T0)的每列有且只有1个整数点元素，设定循环因子k＝0，包括：
由于坐标矩阵表示降采样后数据点的序号，且只有当坐标矩阵为整数的元素表示的数据才能参与波束形成运算，所以只需要找出坐标矩阵中的整数元素，然后找到该元素表示的数据点，对其进行累加即完成了降采样后的波束形成；
下面考虑坐标矩阵的第0行，其元素依次为：
Figure BDA0002131838310000041
p00(T0)必为整数，且值为0；设下一个为整数的元素为p0j(T0),则
Figure BDA0002131838310000042
其中，
Figure BDA0002131838310000043
表示整数组成的集合，k表示某个整数；若S0为T0和M的最小公倍数，则
Figure BDA0002131838310000044
Figure BDA0002131838310000045
均为整数，此时，p0j(T0)为整数；若j'＝kj0,k＝0,1,2,…，则p0j'(T0)亦为整数；所以，坐标矩阵的第0行中存在整数的列编号呈周期性变化，周期为
Figure BDA0002131838310000046
所述依次取出P(k)(T0)中第0～M-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束
Figure BDA0002131838310000047
的第kM,kM+1,…,kM+M-1个点，包括：
P(T0)的第0行第j列出现了整数坐标，则考察P(T0)的第(j+c)列第i行的整数坐标的分布情况，由周期性可知，只需考察第j列到第2j列之间的范围即可，故j＜j+c＜2j，pi(j+c)(T0)为：
Figure BDA0002131838310000048
由于
Figure BDA0002131838310000049
若希望
Figure BDA00021318383100000410
则令
Figure BDA00021318383100000411
i＝kM-cT0，又由于 0≤i≤M-1，得到：
Figure BDA0002131838310000051
若P(T0)的第0行第j列出现了一个整数坐标点，则P(T0)的第(j+c)列第i＝kM-cT0行也将出现一个整数坐标点；
若P(T0)的第(j+c)列第i行出现了一个整数点坐标，下面考察第(j+c)列中整数点坐标出现的情况，设以第(j+c)列第i行为起点，向上或者向下移动m行出现下一个整数点，则该整数点的坐标为第(i+m)行，第(j+c)列，于是此处数值为：
Figure BDA0002131838310000052
由于
Figure BDA0002131838310000053
为整数，故只需要求
Figure BDA0002131838310000054
为整数；所以：
m＝kM,k＝0,±1,±2,…
m的取值说明，在第(j+c)列中整数点出现的位置是周期的，周期为M，则：
定义坐标矩阵P(T0)的第kM～kM+M-1行为坐标矩阵的子矩阵，记为P(k)(T0)，称k＝0时的子矩阵P(0)(T0)为坐标矩阵的主值区间。
所述对P(k)(T0)中所有整数点数值加1，得到下一个周期的坐标矩阵P(k+1)(T0)；更新循环因子：k←k+1，包括：
主值区间内有且只有N个整数点，可求得P(T0)中相邻整数点的数值为：
Figure BDA0002131838310000055
即同一列中相邻整数点之间的数值总是相差1；
在坐标矩阵中，每行出现整数点是周期的，周期为
Figure BDA0002131838310000056
每列出现整数点的位置是周期的，周期为M，主值区间内每列有且必有1个整数点；由于坐标矩阵是周期的，并且相邻整数点数值相差1，故只需要计算主值区间内的数值，即能得到坐标矩阵全部区域的值。
所述计算滤波器通带范围，包括：
滤波器的下截止频率为wfL，上截止频率为wfH，则wfL和wfH应当满足以下约束：对于实信号和复信号wfL应当满足：
Figure BDA0002131838310000057
对于实信号的上截止频率应当满足：
Figure BDA0002131838310000058
对于复解析信号wfH应当满足：
Figure BDA0002131838310000059
其中，w0L和w0H是采样率为fs时信号通带的上截止频率和下截止频率。
本发明的有益效果在于：
本发明提出了坐标矩阵计算方法，并利用坐标矩阵直接对有效数据进行计算，略去了大部分的冗余计算量，从而达到快速计算的目的。
附图说明
图1是本发明的流程图；
具体实施方式
下面结合附图对本发明做进一步描述。
针对升采样导致计算量成倍地增加的问题，本发明提出了一种基于空时二维滤波的快速宽带时延波束形成方法。一方面，该方法能够得到任意原采样率整数倍数的采样率下的输出波束，从而能够保证输出波束的精度。另一方面，该方法在保证输出波束质量的前提下几乎没有增大计算量。实际上，对于传统方法而言，升采样M倍时计算量也提升了M倍。而本文所述方法为传统方法的
Figure BDA0002131838310000061
即本发明所述方法在不增加计算量的前提下完成了升采样M倍的任务，计算量与升采样倍数无关且保持恒定。
下面对本发明方法的原理进行阐述。
假设接收系统中共存在N个传感器(又称N个通道)，依次编号为0,1,…,N-1。第j个传感器的输出为
Figure BDA0002131838310000062
i＝0,1,…,T-1，j＝0,1,…,N-1。其中，i表示接收数据的序号，每个通道中共采集了T个数据点。时延波束形成通过在相邻通道之间插入时延补偿的方式完成波束的扫描。在数字系统中，插入的时延只能为1个采样周期的整数倍。设系统采样率为Mfs，其中M为整数。设相邻通道之间需要插入
Figure BDA0002131838310000063
个采样点，则输出波束表达式为
Figure BDA0002131838310000064
其中，
Figure BDA0002131838310000065
表示第j个通道中序号为
Figure BDA0002131838310000066
的采样点参与了累加，得到了输出波束的第i个数据点。
Figure BDA0002131838310000067
的表达式为
Figure BDA0002131838310000068
显然，
Figure BDA0002131838310000069
必为整数。由萘奎斯特采样定理可知，只要采样率大于信号最高频率的2倍，则可以通过插值的方式恢复信号而无信息丢失。若阵列信号中，目标信号最高频率
Figure BDA00021318383100000610
则可以对接收数据进行M倍降采样处理，利用补零和时延求和得到波束
Figure BDA0002131838310000071
然后再通过低通滤波无失真地恢复输出波束，具体过程表示为
Figure BDA0002131838310000072
此处，
Figure BDA0002131838310000073
为
Figure BDA0002131838310000074
其中，
Figure BDA0002131838310000075
为降采样M倍的信号，
Figure BDA0002131838310000076
为第j个通道参与波束形成的数据序号。
Figure BDA0002131838310000077
与pij(T0)关系为
Figure BDA0002131838310000078
显然，pij(T0)可能不为整数。非整数的pij(T0)表示序号为
Figure BDA0002131838310000079
的数据点不参与降采样后的累加过程，即在降采样后，这部分数据点可以不参与运算。
定义P(T0)为
pij(T0)＝{P(T0)}ij,(i＝0,1,…,T-1,j＝0,1,…,N-1)
则P(T0)为T×N的矩阵，称为坐标矩阵。为描述方便，将坐标矩阵的行从上至下依次称为第0行，第1行，…，第T-1行，同时将坐标矩阵的列从左至右依次称为第0列，第1列，…第N-1列。于是，P(T0)的第i行，第j列元素为pij(T0)。
由于坐标矩阵表示降采样后数据点的序号，且只有当坐标矩阵为整数的元素表示的数据才能参与波束形成运算，所以只需要找出坐标矩阵中的整数元素，然后找到该元素表示的数据点，对其进行累加即可完成降采样后的波束形成。由于坐标矩阵中只有部分元素为整数，故只有部分数据点参与了运算，从而节约了计算量，提高了计算速度。
下面考虑坐标矩阵的第0行。其元素依次为
Figure BDA00021318383100000710
显然p00(T0)必为整数，值为0。设下一个为整数的元素为p0j(T0),则
Figure BDA00021318383100000711
其中，
Figure BDA00021318383100000712
表示整数组成的集合，k表示某个整数。若S0为T0和M的最小公倍数，则
Figure BDA00021318383100000713
Figure BDA0002131838310000081
均为整数。此时，p0j(T0)为整数。若j'＝kj0,k＝0,1,2,…，则p0j'(T0)亦为整数。所以，坐标矩阵的第0行中存在整数的列编号呈周期性变化，周期为
Figure BDA0002131838310000082
下面利用反证法证明在第0列和第
Figure BDA0002131838310000083
列之间不存在另一个整数坐标。设p0j'(T0)为整数，且0＜j'＜j0，则
Figure BDA0002131838310000084
于是有
S＝j'T0＝k'M
所以，S为T0和M的公倍数，于是应当有
Figure BDA0002131838310000085
这与假设矛盾。故P(T0)矩阵的第0行的整数坐标周期性地出现，周期为
Figure BDA0002131838310000086
设P(T0)的第0行第j列出现了整数坐标，下面考察P(T0)的第(j+c)列第i行的整数坐标的分布情况。由周期性可知，只需考察第j列到第2j列之间的范围即可，故j＜j+c＜2j。pi(j+c)(T0)为
Figure BDA0002131838310000087
由于
Figure BDA0002131838310000088
若希望
Figure BDA0002131838310000089
则当令：
Figure BDA00021318383100000810
于是有
i＝kM-cT0
由于0≤i≤M-1，所以得到
Figure BDA00021318383100000811
这表示，若P(T0)的第0行第j列出现了一个整数坐标点，则P(T0)的第(j+c)列第 i＝kM-cT0行也将出现一个整数坐标点。
下面证明，若p0j(T0)为整数，则对于j＜j+c＜2j，0≤i≤M-1，一定存在pi(j+c)(T0)为整数。
证明：
Figure BDA0002131838310000091
由于
Figure BDA0002131838310000092
则只需证
Figure BDA0002131838310000093
设cT0＝kM+r，其中
Figure BDA00021318383100000910
0≤r≤M-1。取i＝M-r，则有
cT0+i＝(k+1)M
证毕。
若P(T0)的第(j+c)列第i行出现了一个整数点坐标，下面考察第(j+c)列中整数点坐标出现的情况。设以第(j+c)列第i行为起点，向上或者向下移动m行出现下一个整数点，则该整数点的坐标为第(i+m)行，第(j+c)列。于是此处数值为
Figure BDA0002131838310000095
由于
Figure BDA0002131838310000096
为整数，故只需要求
Figure BDA0002131838310000097
为整数。所以，可以得到
m＝kM,k＝0,±1,±2,…
m的取值说明，在第(j+c)列中整数点出现的位置是周期的，周期为M。结合上一结论可知，在0～M-1行中，第(j+c)列有且只有1个整数点。为了叙述方便，作出如下定义：
定义坐标矩阵P(T0)的第kM～kM+M-1行为坐标矩阵的子矩阵，记为P(k)(T0)。称k＝0时的子矩阵P(0)(T0)为坐标矩阵的主值区间。
于是，主值区间内有且只有N个整数点。可以求得P(T0)中相邻整数点的数值为
Figure BDA0002131838310000098
即，同一列中相邻整数点之间的数值总是相差1。
综上所述，在坐标矩阵中，每行出现整数点是周期的，周期为
Figure BDA0002131838310000099
每列出现整数点的位置是周期的，周期为M。主值区间内每列有且必有1个整数点。
利用此结论，可以计算出利用坐标矩阵进行波束形成的计算量。由于坐标矩阵是周期的，并且相邻整数点数值相差1，故只需要计算主值区间内的数值，即可得到坐标矩阵全部区域的值。
对于经典方法需要在采样率为Mfs情况下对所有数据点进行求和，所以主值区间内共需要计算M×N次加法。若按照本文方法，只需要在采样率为fs时对坐标矩阵中的整数点进行求和，则共需要计算N次加法。于是，利用坐标矩阵，可以将时延波束形成的计算量降低至原方法的
Figure BDA0002131838310000101
并且可以看到，利用坐标矩阵进行波束形成的方法，其计算量与M无关。这个特点允许将采样率提升至任意倍数而不增加计算量。
所以，只要确定了M、T0与N的值，便可以计算得到完整的坐标矩阵中的整数点。在实现过程中，考虑到数据存储的成本，可以利用坐标矩阵的周期性与相邻整数点数值相差1的特点，只需要存储主值区间内的数值，便可推知坐标矩阵任意位置的数值。据此，便可以设计基于空时二维联合滤波的快速时延波束形成算法。
下面考察从
Figure BDA0002131838310000102
得到y(i)的方法。本文所述方法虽然没有对数据补零的操作，但本文方法实际上利用坐标矩阵进行波束形成相当于同时完成了对数据补零并对其累加的操作。所以，也需要设计一个滤波器对波束进行恢复。下面将通过计算得出滤波器应当满足的设计指标。
对于离散信号，当采样率为时fs，设接收信号为x(n)，于是在z域有
Figure BDA0002131838310000103
频谱为
Figure BDA0002131838310000104
其中，j为复数单位，满足关系为j2＝-1。
利用坐标矩阵进行时延波束形成实际上相当于在x(n)的相邻采样点中插入(M-1)个零，得到序列
Figure BDA0002131838310000105
Figure BDA0002131838310000106
可以表示为
Figure BDA0002131838310000107
于是，序列
Figure BDA0002131838310000108
的z变换为
Figure BDA0002131838310000111
所以，
Figure BDA0002131838310000112
的频谱
Figure BDA0002131838310000113
为
Figure BDA0002131838310000114
结合上述公式，可以得到补零前后信号频谱的关系为
Figure BDA0002131838310000115
观察到式中存在关系：
Figure BDA0002131838310000116
由于X(w)为周期函数，最小正周期为2π，则
Figure BDA0002131838310000117
Figure BDA0002131838310000118
也是周期的，最小正周期为
Figure BDA0002131838310000119
上式表明了补零前后信号的频谱变化关系以及补零后信号频谱的周期性质。利用这两个等式可以实现对补零后频谱的清晰刻画，同时能够为滤波器设计提供指导。
若x(n)存在于频带w0L～w0H之间，则补零后频带变为
Figure BDA00021318383100001110
同时，由于
Figure BDA00021318383100001111
具有周期性，最小正周期为
Figure BDA00021318383100001112
则
Figure BDA00021318383100001113
中也将出现周期性延拓，频谱范围表达式为
Figure BDA00021318383100001114
对于实信号，由于频谱存在对称性，所以x(n)在频率-w0H～(-w0L)之间也存在能量分布，且与w0L～w0H之间关于w＝0对称。此时，
Figure BDA00021318383100001115
中也将出现另一种周期性延拓，频谱范围表达式为
Figure 1
Figure BDA00021318383100001117
的频谱分布范围为上式的叠加。不妨称k＝0时表示的信号为基频信号，|k|＞0时表示的信号称为谐波。可以看到，补零操作使信号出现了谐波，且谐波的频率取决于参数M。在宽带时延波束形成中，基频信号为希望得到的信号，而谐波信号则为干扰，应当避免。
本发明中，采用滤波器将谐波滤除。下面设计滤波器的指标参数。
由于目标信号存在于k＝0表示的区间内，而谐波信号存在于|k|＞0表示的区间内。所以，可以设计一个低通或者带通类型的滤波器，将|k|＞0的部分都滤除。设滤波器的通带范围是 wfL～wfH,不失一般性，可以限制wfL≥0。对于wfL，为了使信号不失真，则应当要求
Figure BDA0002131838310000121
下面考察wfH的选择。下表列出了信号类型与
Figure BDA0002131838310000122
通带的关系。
Figure BDA0002131838310000123
从表中可以看出，当信号为实信号时，wfH应当满足的限制以及滤波器带宽为
Figure BDA0002131838310000124
Figure BDA0002131838310000125
当信号为复解析信号时，wfH应当满足的限制以及滤波器带宽为
Figure BDA0002131838310000126
Figure BDA0002131838310000127
显然，复解析信号可以允许更大的过渡带。更宽的过渡带意味着对滤波器要求更低，允许使用更低阶次的滤波器，从而尽可能降低滤波器引入的信号畸变。由于希尔伯特变换能够将实信号变换为复解析信号，则可以借用希尔伯特变换得到复解析信号，从而降低对滤波器的要求。
同时，从式(35)和(36)中也可以看到，对于同一个信号，升采样倍数不能太高。因为随着 M的增大，Δwr和Δwi将逐渐变小。这意味着，随着M的增大，也需要逐渐提高滤波器阶数以获得更窄的过渡带才能够得到基频信号。而随着滤波器阶数的增加，滤波器引起的信号畸变也变得更加严重。所以，升采样倍数M不能设置得任意大，否则将因为滤波器过渡带过窄导致信号严重畸变，或者谐波进入滤波器通带，导致波束形成失败。
一种空时二维滤波的快速宽带波束形成方法，提出了步骤1～10，使得系统可以任意提升采样率而几乎不增加计算量。
滤波器设计指标须满足
Figure BDA0002131838310000131
和
Figure BDA0002131838310000132
其中，M为升采样倍数。 w0L,w0H分别为X(w)的下截止角频率和上截止角频率，wfL,wfH分别为多项滤波器的下截止角频率和上截止角频率。
提出了坐标矩阵计算方法，并利用坐标矩阵直接对有效数据进行计算，略去了大部分的冗余计算量，从而达到快速计算的目的。
步骤1：以符合萘奎斯特(Niquest)定理的采样率进行采样，每个通道采集的数据按时间顺序排列成列向量，组成一个T×N维的数据矩阵，列的编号依次为j＝0,1,…,N-1，行的编号依次为i＝0,1,…,T-1。
步骤2：确定升采样倍数M以及相邻通道时延补偿点数T0。
步骤3：计算坐标矩阵主值区间P(0)(T0)。规定坐标矩阵的列编号从左至右依次为j＝0,1,…,N-1，行编号从上至下依次为i＝0,1,…,M-1。则P(0)(T0)的第i行第j列元素为
Figure BDA0002131838310000133
步骤4：去除P(0)(T0)中非整数点，只保留整数点。P(0)(T0)的每列有且只有1个整数点元素。设定循环因子k＝0。
步骤5：依次取出P(k)(T0)中第0～M-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束
Figure BDA0002131838310000134
的第kM,kM+1,…,kM+M-1个点。由于P(k)(T0)中共有M×N个元素，但只存在N个元素为整数，故利用P(k)(T0)将加法次数从M×N次降至N次，计算量为传统方法的
Figure BDA0002131838310000135
倍(对于传统方法而言，若将采样率提升M倍，则计算量亦提升M倍)。于是，本方法的在时延波束形成阶段的计算量与未升采样时的计算量相同，与升采样倍数无关。故本方法在保持计算量不变的情况下可以任意提升采样率。
步骤6：对P(k)(T0)中所有整数点数值加1，得到下一个周期的坐标矩阵P(k+1)(T0)。更新循环因子：k←k+1。
步骤7：判断kM≥T是否成立。若成立，则说明已经遍历数据矩阵中所有数值点，执行步骤8，否则执行步骤5。
步骤8：计算滤波器通带范围。设滤波器的下截止频率为wfL，上截止频率为wfH。则wfL和wfH应当满足以下约束：对于实信号和复信号wfL应当满足：
Figure BDA0002131838310000141
对于实信号的上截止频率应当满足
Figure BDA0002131838310000142
对于复解析信号wfH应当满足，
Figure BDA0002131838310000143
其中，w0L和w0H是采样率为fs时信号通带的上截止频率和下截止频率。
步骤9：根据步骤8中计算的参数设计滤波器对
Figure BDA0002131838310000144
进行滤波，得到输出波束y(i)。
步骤10：结束。

Claims (5)

1.一种空时二维滤波的快速宽带波束形成方法，其特征在于，包括以下步骤：
(1)假设接收系统中共存在N个传感器，又称N个通道，依次编号为0,1,…,N-1；第j个传感器的输出为
Figure FDA0003809275740000011
其中，i表示接收数据的序号，每个通道中共采集了T个数据点，时延波束形成通过在相邻通道之间插入时延补偿的方式完成波束的扫描，在数字系统中，插入的时延只能为1个采样周期的整数倍，设系统采样率为Mfs，其中M为整数，设相邻通道之间需要插入
Figure FDA0003809275740000012
个采样点，则输出波束表达式为
Figure FDA0003809275740000013
其中，
Figure FDA0003809275740000014
表示第j个通道中序号为
Figure FDA0003809275740000015
的采样点参与了累加，得到了输出波束的第i个数据点，
Figure FDA0003809275740000016
的表达式为：
Figure FDA0003809275740000017
显然，
Figure FDA0003809275740000018
必为整数，由萘奎斯特采样定理可知，只要采样率大于信号最高频率的2倍，则可以通过插值的方式恢复信号而无信息丢失，若阵列信号中，目标信号最高频率
Figure FDA0003809275740000019
则可以对接收数据进行M倍降采样处理，利用补零和时延求和得到波束
Figure FDA00038092757400000110
然后再通过低通滤波无失真地恢复输出波束，具体过程表示为
Figure FDA00038092757400000111
此处，
Figure FDA00038092757400000112
为
Figure FDA00038092757400000113
其中，xj(i)为降采样M倍的信号，表示第j个通道中序号为
Figure FDA00038092757400000114
的采样点参与了累加，
Figure FDA00038092757400000115
与pij(T0)关系为
Figure FDA00038092757400000116
显然，pij(T0)可能不为整数，非整数的pij(T0)表示序号为
Figure FDA00038092757400000117
的数据点不参与降采样后的累加过程，即在降采样后，这部分数据点可以不参与运算；
定义P(T0)为pij(T0)＝{P(T0)}ij,(i＝0,1,…,T-1,j＝0,1,…,N-1)
则P(T0)为T×N的矩阵，称为坐标矩阵，为描述方便，将坐标矩阵的行从上至下依次称为第0行，第1行，…，第T-1行，同时将坐标矩阵的列从左至右依次称为第0列，第1列，…第N-1列，于是，P(T0)的第i行，第j列元素为pij(T0)；
(2)确定采样倍数M以及相邻通道时延补偿点数T0；
(3)计算坐标矩阵主值区间P(0)(T0)，规定坐标矩阵的列编号从左至右依次为j＝0,1,…,N-1，行编号从上至下依次为i＝0,1,…,M-1，则P(0)(T0)的第i行第j列元素为
Figure FDA0003809275740000021
(4)去除P(0)(T0)中非整数点，只保留整数点；P(0)(T0)的每列有且只有1个整数点元素，设定循环因子k＝0；
(5)依次取出P(k)(T0)中第0～M-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束
Figure FDA0003809275740000022
的第kM,kM+1,…,kM+M-1个点；
(6)对P(k)(T0)中所有整数点数值加1，得到下一个周期的坐标矩阵P(k+1)(T0)；更新循环因子：k←k+1；
(7)判断kM≥T是否成立，若成立，则说明已经遍历数据矩阵中所有数值点，执行步骤(8)，否则执行步骤(5)；
(8)计算滤波器通带范围；
(9)根据步骤(8)中滤波器通带范围设计滤波器对
Figure FDA0003809275740000023
进行滤波，得到输出波束y(i)。
2.根据权利要求1所述的一种空时二维滤波的快速宽带波束形成方法，其特征在于，所述去除P(0)(T0)中非整数点，只保留整数点；P(0)(T0)的每列有且只有1个整数点元素，设定循环因子k＝0，包括：
由于坐标矩阵表示降采样后数据点的序号，且只有当坐标矩阵为整数的元素表示的数据才能参与波束形成运算，所以只需要找出坐标矩阵中的整数元素，然后找到该元素表示的数据点，对其进行累加即完成了降采样后的波束形成；
下面考虑坐标矩阵的第0行，其元素依次为：
Figure FDA0003809275740000031
p00(T0)必为整数，且值为0；设下一个为整数的元素为p0j(T0),则
Figure FDA0003809275740000032
其中，
Figure FDA0003809275740000033
表示整数组成的集合，q表示某个整数；若S0为T0和M的最小公倍数，则
Figure FDA0003809275740000034
Figure FDA0003809275740000035
均为整数，此时，p0j(T0)为整数；若j′＝qj0,q＝0,1,2,…，则p0j′(T0)亦为整数；所以，坐标矩阵的第0行中存在整数的列编号呈周期性变化，周期为
Figure FDA0003809275740000036
3.根据权利要求1所述的一种空时二维滤波的快速宽带波束形成方法，其特征在于，所述依次取出P(k)(T0)中第0～M-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束
Figure FDA0003809275740000037
的第kM,kM+1,…,kM+M-1个点，包括：
P(T0)的第0行第j列出现了整数坐标，则考察P(T0)的第(j+c)列第i行的整数坐标的分布情况，由周期性可知，只需考察第j列到第2j列之间的范围即可，故j＜j+c＜2j，pi(j+c)(T0)为：
Figure FDA0003809275740000038
由于
Figure FDA0003809275740000039
若希望
Figure FDA00038092757400000310
则令
Figure FDA00038092757400000311
又由于0≤i≤M-1，得到：
Figure FDA00038092757400000312
若P(T0)的第0行第j列出现了一个整数坐标点，则P(T0)的第(j+c)列第i＝kM-cT0行也将出现一个整数坐标点；
若P(T0)的第(j+c)列第i行出现了一个整数点坐标，下面考察第(j+c)列中整数点坐标出现的情况，设以第(j+c)列第i行为起点，向上或者向下移动m行出现下一个整数点，则该整数点的坐标为第(i+m)行，第(j+c)列，于是此处数值为：
Figure FDA00038092757400000313
由于
Figure FDA00038092757400000314
为整数，故只需要求
Figure FDA00038092757400000315
为整数；所以：
m＝kM,k＝0,±1,±2,…
m的取值说明，在第(j+c)列中整数点出现的位置是周期的，周期为M，则：
定义坐标矩阵P(T0)的第kM～kM+M-1行为坐标矩阵的子矩阵，记为P(k)(T0)，称k＝0时的子矩阵P(0)(T0)为坐标矩阵的主值区间。
4.根据权利要求1所述的一种空时二维滤波的快速宽带波束形成方法，其特征在于，所述对P(k)(T0)中所有整数点数值加1，得到下一个周期的坐标矩阵P(k+1)(T0)；更新循环因子：k←k+1，包括：
主值区间内有且只有N个整数点，可求得P(T0)中相邻整数点的数值为：
Figure FDA0003809275740000041
即同一列中相邻整数点之间的数值总是相差1；
在坐标矩阵中，每行出现整数点是周期的，周期为
Figure FDA0003809275740000042
每列出现整数点的位置是周期的，周期为M，主值区间内每列有且必有1个整数点；由于坐标矩阵是周期的，并且相邻整数点数值相差1，故只需要计算主值区间内的数值，即能得到坐标矩阵全部区域的值。
5.根据权利要求1所述的一种空时二维滤波的快速宽带波束形成方法，其特征在于，所述计算滤波器通带范围，包括：
滤波器的下截止频率为wfL，上截止频率为wfH，则wfL和wfH应当满足以下约束：对于实信号和复信号wfL应当满足：
Figure FDA0003809275740000043
对于实信号的上截止频率应当满足：
Figure FDA0003809275740000044
对于复解析信号wfH应当满足：
Figure FDA0003809275740000045
其中，w0L和w0H是采样率为fs时信号通带的上截止频率和下截止频率。
CN201910640904.3A 2019-07-16 2019-07-16 一种空时二维滤波的快速宽带波束形成方法 Active CN110333504B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910640904.3A CN110333504B (zh) 2019-07-16 2019-07-16 一种空时二维滤波的快速宽带波束形成方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910640904.3A CN110333504B (zh) 2019-07-16 2019-07-16 一种空时二维滤波的快速宽带波束形成方法

Publications (2)

Publication Number Publication Date
CN110333504A CN110333504A (zh) 2019-10-15
CN110333504B true CN110333504B (zh) 2022-11-18

Family

ID=68146791

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910640904.3A Active CN110333504B (zh) 2019-07-16 2019-07-16 一种空时二维滤波的快速宽带波束形成方法

Country Status (1)

Country Link
CN (1) CN110333504B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111260782B (zh) * 2020-01-16 2023-10-20 广东工业大学 离散数据集的插值数据计算与读取方法、系统及存储介质

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1866356A (zh) * 2005-08-15 2006-11-22 华为技术有限公司 一种宽带波束形成方法和装置

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6110115A (en) * 1998-12-17 2000-08-29 Acuson Corporation Medical diagnostic ultrasonic imaging system and method
US6952460B1 (en) * 2001-09-26 2005-10-04 L-3 Communications Corporation Efficient space-time adaptive processing (STAP) filter for global positioning system (GPS) receivers
KR100898082B1 (ko) * 2003-12-24 2009-05-18 노키아 코포레이션 상보적 노이즈 분리 필터를 이용한 효율적 빔포밍 방법
US8085871B2 (en) * 2005-04-21 2011-12-27 Broadcom Corporation Adaptive modulation in a multiple input multiple output wireless communication system with optional beamforming
CN101609150B (zh) * 2009-07-07 2011-09-14 哈尔滨工程大学 一种提高基阵分辨力和增益的快速波束形成方法
CN102608588B (zh) * 2012-03-14 2014-04-16 西安电子科技大学 基于子带分解的宽带子阵自适应波束形成方法
CN109004970A (zh) * 2018-04-18 2018-12-14 哈尔滨工程大学 一种零范数约束的自适应稀疏阵列波束形成方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1866356A (zh) * 2005-08-15 2006-11-22 华为技术有限公司 一种宽带波束形成方法和装置

Also Published As

Publication number Publication date
CN110333504A (zh) 2019-10-15

Similar Documents

Publication Publication Date Title
CN107315160B (zh) 基于内插虚拟阵列信号原子范数最小化的互质阵列波达方向估计方法
CN107102291B (zh) 基于虚拟阵列内插的无网格化互质阵列波达方向估计方法
CN107290709B (zh) 基于范德蒙分解的互质阵列波达方向估计方法
CN108802705B (zh) 一种基于稀疏的空时自适应处理方法及系统
CN110333504B (zh) 一种空时二维滤波的快速宽带波束形成方法
CN109901101A (zh) 基于电磁矢量传感器互质阵列相干信号到达角估计方法
CN101902258A (zh) 一种数字预失真处理参数求取的方法及装置
CN107561484A (zh) 基于内插互质阵列协方差矩阵重建的波达方向估计方法
JP2812679B2 (ja) トランルスバーサルフイルタの形式のデジタル信号処理方法及び装置
CN105490665B (zh) 一种最优指数幂多项式插值滤波器系数的计算方法
CN106603036A (zh) 一种基于低阶内插滤波器的自适应时延估计方法
CN110361733B (zh) 一种基于时频联合重采样的中轨sar大斜视成像方法
Veligosha et al. Adjustment of adaptive digital filter coefficients in modular codes
CN110749856B (zh) 一种基于零化去噪技术的互质阵欠定测向方法
CN113221313A (zh) 一种音频均衡器的设计方法
CN115801508A (zh) 基于自适应分段小数滤波与频域均衡的多通道校正方法
CN109061564B (zh) 基于高阶累积量的简化近场定位方法
CN113054948A (zh) 并行卡尔曼滤波系统和方法
CN106154245A (zh) 基于等效阵列方向图的集中式mimo雷达阵列设计方法
US5233549A (en) Reduced quantization error FIR filter
CN114189229B (zh) 一种基于自适应分段算法的小数延时滤波方法
Han et al. A new architecture for ultrasound sigma-delta modulation beamformer
CN114036975B (zh) 基于频域-波数域解卷积的目标信号提取方法
KR101581686B1 (ko) 다목적 빔 집속 시스템
CN110504948B (zh) 信号处理方法及等间隔分数时延滤波器

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant