CN103999151A - 计算上有效的宽带滤波和相加阵列聚焦 - Google Patents

计算上有效的宽带滤波和相加阵列聚焦 Download PDF

Info

Publication number
CN103999151A
CN103999151A CN201280053952.9A CN201280053952A CN103999151A CN 103999151 A CN103999151 A CN 103999151A CN 201280053952 A CN201280053952 A CN 201280053952A CN 103999151 A CN103999151 A CN 103999151A
Authority
CN
China
Prior art keywords
sensor
output signal
filter
interpolation
signal
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.)
Granted
Application number
CN201280053952.9A
Other languages
English (en)
Other versions
CN103999151B (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.)
Hottinger Bruel and Kjaer AS
Original Assignee
Bruel and Kjaer Sound and Vibration Measurement AS
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 Bruel and Kjaer Sound and Vibration Measurement AS filed Critical Bruel and Kjaer Sound and Vibration Measurement AS
Publication of CN103999151A publication Critical patent/CN103999151A/zh
Application granted granted Critical
Publication of CN103999151B publication Critical patent/CN103999151B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10KSOUND-PRODUCING DEVICES; METHODS OR DEVICES FOR PROTECTING AGAINST, OR FOR DAMPING, NOISE OR OTHER ACOUSTIC WAVES IN GENERAL; ACOUSTICS NOT OTHERWISE PROVIDED FOR
    • G10K11/00Methods or devices for transmitting, conducting or directing sound in general; Methods or devices for protecting against, or for damping, noise or other acoustic waves in general
    • G10K11/18Methods or devices for transmitting, conducting or directing sound
    • G10K11/26Sound-focusing or directing, e.g. scanning
    • G10K11/34Sound-focusing or directing, e.g. scanning using electrical steering of transducer arrays, e.g. beam steering
    • G10K11/341Circuits therefor
    • G10K11/346Circuits therefor using phase variation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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
    • G01S3/00Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/80Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
    • G01S3/802Systems for determining direction or deviation from predetermined direction
    • G01S3/808Systems for determining direction or deviation from predetermined direction using transducers spaced apart and measuring phase or time difference between signals therefrom, i.e. path-difference systems
    • G01S3/8083Systems for determining direction or deviation from predetermined direction using transducers spaced apart and measuring phase or time difference between signals therefrom, i.e. path-difference systems determining direction of source

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Acoustics & Sound (AREA)
  • Multimedia (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Circuit For Audible Band Transducer (AREA)

Abstract

本发明提供了一种确定包括多个传感器的传感器阵列的聚焦输出信号的方法,每个传感器可操作为输出响应于被测量的量的传感器输出信号,聚焦输出信号表示所计算的在聚焦点处的量;该方法包括:从每个传感器中接收相应测量的传感器输出信号;通过执行关于所测量的传感器信号的聚焦计算来计算聚焦的输出信号;其中,该方法进一步包括确定一组预定网格点的网格点子组,每个网格点具有与其相关的至少一个预先计算的滤波器参数;并且其中,计算聚焦输出信号包括执行关于网格点子组的插值,以便获得插值的聚焦输出信号。

Description

计算上有效的宽带滤波和相加阵列聚焦
技术领域
本发明涉及阵列信号处理,例如,声学信号的宽带波束形成。
背景技术
宽带波束形成是一种被广泛应用于定向接收诸如声学或无线电波的技术。已经在例如声源定位、声纳、雷达、无线通信等的背景下描述了波束形成技术。通常,在这种系统中,以使所产生的测量系统对来自特定方向的波特别敏感的方式放大和延迟传感器的信号。在这种测量系统中,能够将传感器阵列的敏感性引向特定方向-称为“波束形成”的处理。当同时记录所有信道时,这种系统仅需要非常少量的时间来用于单个测量。
例如,很多噪声消除问题涉及在复杂环境(诸如,汽车或飞机内部)中定位一个或多个噪声源。近年来,已经可以同时使用多个信道来执行测量。如今,存在具有安装在网格中的大量麦克风(例如,64或128)的测量系统。在其他测量系统中,麦克风通常安装在不规则的配置中。
由于麦克风(或其他传感器)和数据/信号采集硬件的成本,所以通常希望在波束形成系统中尽可能较少地使用传感器。另一方面,对系统的频率范围和空间精度的要求都趋向于增加在阵列中需要的传感器的数量。
在所谓的滤波和相加(Filter-And-Sum)(FAS)波束形成中,通过将单独的滤波器应用于传感器信号并随后将滤波的信号相加来计算在给定位置处的输出时间信号。Shefeng Yan等人的文章“Convex optimizationbased time-domain broadband beamforming with sidelobe control”(J.Acoust.Soc.Am.121(1),2007,1)描述了一种基于FIR滤波器的方法和一种用于例如为最小化波束形成器的旁瓣电平而优化FIR滤波器的方法。
优化的宽带滤波和相加(FAS)波束形成能够相对于延迟和相加(Delay-And-Sum)(DAS)波束形成大幅地减少旁瓣电平,并且对于球形阵列-球谐波束形成(SHB)的情况,例如,参照Shefeng Yan等人(同前)和Shefeng Yan等人的文章“Optimal Modal Beamforming for SphericalMicrophone Arrays”(IEEE Transactions on Audio,Speech and LanguageProcessing,Vol.19,No.2,2011,2,361-371)。
然而,如本文所述,计算最佳滤波器参数是非常繁重的计算任务。尤其地,在其中为多个聚焦点以及为大型传感器阵列执行波束形成操作的应用中,用于执行先有技术的优化宽带FAS波束形成的计算资源可能被禁止。而且,由Shefeng Yan等人(J.Acoust.Soc.Am,同前)概述的方法需要基于来自每个特定的测量中的信号的协方差矩阵进行最优化,引入与最小方差(或Capon)波束形成器的特性相似的特性:在由对两个源的测量所产生的输出并非分别对这两个源的每一个的测量所产生的的输出之和的意义上,输出并非源的线性函数。
发明内容
本文公开了一种确定包括多个传感器的传感器阵列的聚焦输出信号的方法,每个传感器可操作为输出响应于被测量的量的传感器输出信号,聚焦输出信号表示所计算的在聚焦点处的量。该方法的实施方式包括:
-从每个传感器中接收相应测量的传感器输出信号;
-通过执行关于所测量的传感器信号的聚焦计算,来计算聚焦输出信号。
聚焦输出信号可因此被视为表示为相应的一个传感器计算的多个被滤波的传感器信号的组合(例如,求和或其他线性组合);每个被滤波的传感器信号因此可表示通过与相应的传感器相关联的滤波器滤波的来自各个传感器的所测量的传感器信号,其中,从相应的传感器中接收测量的传感器信号,滤波器取决于依赖聚焦点的至少一个滤波器参数。
发明人已经认识到,通过定义跨越其中利用特定阵列执行聚焦计算的空间中的区域的聚焦点网格,利用相当低的计算资源便可获得具有最佳滤波器的波束形成和其他阵列聚焦技术(例如,声全息技术)的好处。由此定义的点也被称为网格点。为这些预定义的网格点的每一个,计算一组最佳滤波器参数,并且预先计算的滤波器参数可存储在与传感器阵列相关联的滤波器库(例如,文件、一组文件或数据库)中。例如,这种滤波器库可与阵列一起设置在合适的存储介质上。当已经进行了测量并且执行了在任意聚焦点r处的聚焦计算时,则随后识别例如作为周围最近的网格点的预定义的网格点子组,并且执行关于所选择的网格点子组的插值。在某些实施方式中,为这些网格点中的每一个执行聚焦计算,并且通过在来自周围网格点的聚焦输出值之间的插值逼近在聚焦点处的聚焦输出。在替代性实施方式中,对所识别的网格点子组的预先计算的滤波器参数插值,并且使用所产生的插值的滤波器参数来执行聚焦计算。因此,插值步骤和某些或所有聚焦计算步骤的顺序可互换。
因此,本文中公开的方法的实施方式进一步包括:确定一组预定网格点的网格点子组,每个网格点具有与其相关的至少一个预先计算的滤波器参数;并且为任意的聚焦点计算聚焦输出信号包括执行关于网格点子组的插值,以便获得插值的聚焦输出信号。尤其地,计算作为插值的聚焦输出信号的聚焦输出信号可包括将一个或多个辅助滤波器应用于相应测量的传感器输出信号的每一个;一个或多个辅助滤波器中的每一个与从中接收相应测量的传感器信号的相应的传感器相关联,并且一个或多个辅助滤波器的每一个取决于至少一个预先计算的滤波器参数。因此,对于与其相关联的每个传感器,每个网格点可具有用于与其相关联的每个传感器的至少一个预先计算的滤波器参数,即,每个网格点可将其与多个预先计算的滤波器参数相关联,每个预先计算的滤波器参数与一个传感器相关联。
本文中公开的方法的实施方式可应用于不同类型的阵列聚焦应用(诸如,波束形成和声全息技术)或者计算在期望的聚焦点处的来自不同传感器位置的贡献的类似应用中。在波束形成的情况下,聚焦输出信号可被称为波束形成的输出信号;它可表示在聚焦点处的声源对在传感器阵列的位置(例如,在传感器阵列的中心)处的被测量的量(例如,声压)的贡献。在这种实施方式中,聚焦计算因此是波束形成计算。
在声全息技术的情况下,聚焦输出信号表示在聚焦点处的量(例如,声压、粒子速度或不同的声学量)(对聚焦点处的量的估计)。在该情况下,不同聚焦点的输出信号可被用于重构(或“后向传播”)在期望区域(例如,物体的表面或体积)中的声场。
因此,在某些实施方式中,聚焦计算是波束形成计算;并且聚焦计算包括确定聚焦点对在传感器阵列的阵列位置处的声学量(例如,声压)的贡献的估计量。在替代性实施方式中,聚焦计算包括确定对在聚焦点处的声场的参数的估计。因此,在这种实施方式中,聚焦输出信号表示重构的声学量。通常,聚焦计算包括基于来自传感器阵列的传感器的测量的传感器信来号计算在聚焦点处的物理量。在某些实施方式中,聚焦计算是包括应用滤波器以获得被滤波的信号以及求和(或其他线性组合)以获得所产生的聚焦信号的滤波和相加计算。因此,在某些实施方式中,术语阵列聚焦旨在表示用于通过滤波和相加运算来估计在聚焦点处的量的处理。
因此,在某些实施方式中,术语聚焦计算是指滤波和相加计算,并且聚焦的输出信号表示滤波和相加计算的结果。滤波和相加计算包括将各个滤波器应用于来自各个传感器的传感器信号以及将被滤波的信号求和。
在某些实施方式中,计算聚焦输出信号包括:
-通过将各个辅助滤波器应用于所测量的传感器输出信号来为每个传感器计算多个辅助滤波传感器信号,每个辅助滤波器与从中接收测量的传感器信号的相应的传感器相关联,并且每个辅助滤波器取决于至少一个预先计算的滤波器参数;每个辅助滤波器因此与网格点子组中的一个相关联;
-为每个网格点组合为相应的传感器计算的多个辅助滤波传感器信号,以获得辅助聚焦输出信号,并且对为相应的网格点计算的辅助聚焦输出信号进行插值以获得所述插值的聚焦输出信号。
在某些实施方式中,计算聚焦输出信号包括:
-通过将各个辅助滤波器应用于所测量的传感器输出信号来为每个传感器计算多个辅助滤波传感器信号,每个辅助滤波器与从中接收测量的传感器信号的相应的传感器相关联,并且每个辅助滤波器取决于至少一个预先计算的滤波器参数;而且,每个辅助滤波器与网格点子组中的一个相关联;
-为每个传感器插值为各个网格点计算的多个辅助滤波传感器信号,以获得插值的滤波传感器信号,并且组合为各个传感器计算的插值的传感器信号,以获得所述插值的聚焦输出信号。
在其他实施方式中,计算聚焦输出信号包括:
-为每个传感器从与确定的网格点子组的相应一个相关联的预先计算的滤波器参数中计算至少一个插值的滤波器参数;
-通过将各个插值的滤波器应用于所测量的传感器输出信号来为每个传感器计算插值的滤波传感器信号;以及
-组合为各个传感器计算的插值的滤波传感器信号,以获得插值的聚焦输出信号。
将理解,可在不同的域(例如,时域和/或频域)中执行聚焦操作和/或插值操作。
在相对于应用的阵列的远场区域中,波束形成器没有距离(深度)分辨率,因此,仅可获得源的方向辨别。在中间距离,可实现某种程度的距离分辨率,因此在此可期望在3D中的特定点处的聚焦,例如,为了获得3D源分布图。本文中描述的方法的实施方式为聚焦点的网格使用预先计算的滤波器。如果波束形成器旨在仅用于远场区域内,那么仅需要方向(2D)网格,并且将进行方向上的聚焦。否则,可能需要3D网格,并且在特定的点上进行聚焦。相似的考虑适用于声全息技术。然而,为了本描述的目的,术语方向聚焦和点聚焦将被互换使用,并且聚焦方向也由在相关方向的位置向量表示。因此,本文中使用的术语网格点和聚焦点旨在包括在3D空间中的位置和在3D空间中的方向,例如,由在单位球面上的位置表示。因此,可由相对于合适的坐标系的三维空间坐标或者可由定义位于3D空间中的二维表面(例如,球面)上的位置的两维空间坐标来定义每个网格点和每个聚焦点。
传感器输出信号可为声学信号(即,表示被测量的声音),例如,噪音、可听见的声音、不可听见的声音(诸如,超音波或次声波等)或其组合。在其他实施方式中,传感器输出信号可表示任何其他声学或无线电信号,例如,声纳信号、雷达信号、无线通信信号等。
所测量的量可为声学量,例如,声压、声压梯度、粒子速度等。因此,每个传感器可为任何合适的听觉测量装置,例如,麦克风、水诊器、压力梯度换能器、粒子速度换能器、用于无线通信系统的接收器/换能器、雷达和/或声纳等或其组合。传感器阵列包括多个传感器,例如,布置在规则或不规则的网格(例如,二维或三维网格)中的一组传感器。
阵列的传感器可被布置在一组测量位置的相应一个上。这组测量位置可被布置在一个或多个测量平面内,例如,在单个平面内或者在两个或多个平行平面内。在每个平面内,测量位置可以不规则的图案或任意其它合适的方式布置在规则的网格上。而且,本文中描述的方法还可应用于非平面测量几何中,即,测量位置不是位于一个或多个平行平面内而是例如位于弯曲的表面上的布置。例如,本文中描述的方法可应用于球形阵列几何中。
术语插值旨在表示用于至少近似地计算在离散的一组已知的数据点附近的新数据点的任何合适的处理。在本描述的情况下,术语插值旨在表示从已知的滤波器参数计算与聚焦点相关联的聚焦输出信号,每个滤波器参数与离散的一组(子组)网格点相关联,其中,聚焦点在网格点的预定附近内。因此,选择网格点子组可包括选择在聚焦点附近的网格点,例如,选择在聚焦点的预定附近的所有网格点,或者选择最靠近聚焦点的预定数量的网格点。然而,将理解,可选择用于选择被用于执行插值的网格点子组的其他方案。
插值可为分段常数插值、线性插值、多项式插值、样条插值或任何其他合适的插值方法。插值可在角度和/或线性坐标(1D、2D或3D)中。将进一步理解,在某些实施方式中,插值和组合操作的顺序可互换。尤其地,这适用于其中插值是在每个网格点上计算的值的线性组合并且其中产生聚焦输出信号的多个滤波传感器信号的组合是所计算的滤波传感器信号的线性组合(例如,所计算的滤波传感器信号之和)的实施方式。
在某些实施方式中,可通过应用与各个网格点相关联的各个预先计算的滤波器参数来计算用于插值(以及用于每个传感器)的每个网格点的辅助滤波信号。随后,对于每个传感器,所选择的聚焦点需要的滤波传感器信号可被计算为所计算的辅助滤波传感器信号的插值。可替代地,可为每个网格点从该网格点和各个传感器的辅助滤波传感器信号中计算辅助聚焦输出信号,并且随后,可对各个网格点的辅助聚焦输出信号进行插值,以获得插值的聚焦输出信号。因此,可按照任何顺序执行对各个传感器的辅助滤波传感器信号的求和(或者另外组合)和对不同的网格点的插值,特别是当插值被计算为与各个网格点相关联的计算值的线性组合时。在例如为了获得波束形成图或重构的声场图而执行关于多个聚焦点的聚焦(例如,波束形成)的实施方式中,可按照以下方式进一步加快处理:首先,可横跨在获得波束形成或重构图中涉及的所有聚焦点上的插值所需要的至少网格点子组执行聚焦计算,即,通过应用与所述网格点子组的每一个相关联的各个预先计算的滤波器参数,以便计算各个辅助滤波信号,每个辅助滤波信号与一个网格点相关并且与一个传感器相关。对于每个网格点,计算多个相应的辅助滤波信号,为每个传感器计算至少一个辅助滤波信号。随后,可通过对所计算的辅助滤波信号的相应一个进行插值来计算用于每个期望的聚焦点的滤波传感器信号。可替代地,可通过插值从各个辅助滤波信号中计算的各个聚焦输出信号来计算用于每个期望的聚焦点的插值的聚焦输出信号。
在被滤波的传感器输出信号是(至少近似为)滤波器参数的线性函数时,例如与FIR滤波器的情况一样,可通过下列步骤计算与期望的聚焦点相关联的被滤波的传感器输出信号:首先通过从与网格点的相应一个相关联的预先计算的滤波器参数计算插值的滤波器参数,插值的滤波器参数定义插值的滤波器,并且随后通过将插值的滤波器应用于传感器输出信号来计算被滤波的传感器输出。可随后通过组合来自各个传感器的滤波传感器输出来计算插值的聚焦输出信号。在该实施方式中,在某些情况下,聚焦计算的数量将减少,从而降低阵列聚焦(例如,波束形成)操作的计算成本。
每个网格点可将其与用于每个传感器的至少一个预先计算的滤波器参数相关联。可通过最小化聚焦输出信号的功率同时要求在输出信号中完整地保留来自网格点的贡献来为给定的网格点计算预先计算的滤波器参数。
在某些实施方式中,可通过最小化最大旁瓣电平来计算与第一网格点相关联的每组预先计算的滤波器参数,最大的旁瓣电平表示来自不同于第一网格点的一组其他位置的并且位于一组预定的频率处的干扰抑制水平(level,电平)。已经认识到,该实施方式降低了在干扰源的某些位置处以及在某些频率处的高旁瓣灵敏度峰值的风险。
通常,每个网格点可具有与其相关联的一组预先计算的滤波器参数。
在某些实施方式中,通过下列步骤计算与第一网格点相关联的每个预先计算的滤波器参数组:
-通过最小化最大旁瓣电平来为一组预定频率中的每一个确定一组复杂的传感器权重,最大的旁瓣电平表示来自不同于第一网格点的一组其他位置的并且位于所述频率处的干扰的抑制水平(level,电平),每个传感器权重与各个传感器相关;
-通过将由所述预先计算的滤波器参数定义/建立的频率响应于由所述确定的传感器权重的子组构成的频率响应向量相配,来为每个传感器确定至少一个预先计算的滤波器参数,所述子组的传感器权重与特定传感器的各个频率相关联。已经证明,该实施方式大幅减少了用于计算预先计算的滤波器参数的所需要的计算资源,同时允许在不同的频率处实现截然不同的旁瓣抑制电平。
网格点可被选择为任何合适的规则的或不规则的点的网格,例如,二维或三维网格。网格点可布置在一个或多个平面中,例如,在单个平面中或在两个或多个平行平面中。在每个平面内,网格点可以不规则的图案或任何合适的方式布置在规则的网格上。可替代地,网格点可分布在一个或多个弯曲的表面(例如,具有相应半径的一个或多个球面)中。在某些实施方式中,网格点被布置为在相邻的网格点之间的距离小于所使用的波束形成器的本地波束宽度。
应注意,上述和以下描述的方法的特征可至少部分以软件或固件来实施并且在数据处理装置或由程序代码方法(诸如,计算机可执行的指令)的执行引起的其他处理器件上执行。在此处和下文中,术语处理器件包括适当地适配于执行以上功能的任何电路和/或装置。尤其地,以上术语包括通用或专用可编程微处理器、数字信号处理器(DSP)、专用集成电路(ASIC)、可编程逻辑阵列(PLA)、现场可编程门阵列(FPGA)、图形处理单元(GPU)、专用电子电路等或其组合。
可以不同的方式实施包括在上面和下面描述的方法、系统、装置和产品器件的本发明的实施方式,其每一个产生结合首次提及的方法描述的一个或多个好处和优点,并且每一个具有与结合首次提及的方法描述的和/或在从属权利要求中公开的实施方式对应的一个或多个实施方式。
尤其地,用于执行阵列聚焦(例如,波束形成)计算的处理设备的实施方式可包括:接口,用于从传感器阵列的各个传感器中接收响应于被测量的量的一组传感器输出信号;以及处理单元,被配置为执行在本文中定义的方法的实施方式的步骤,其中,处理单元包括用于存储一组预先计算的滤波器参数的存储介质。
一种阵列聚焦(例如,波束形成器或全息)系统可包括:如在上面和下面定义的处理设备;以及一组传感器,用于在一组测量位置处测量被测量的量,并且在通信连接中可连接至该设备,以便将被测量的第一声学量转发给处理设备。例如,这种系统可被用于在3D空间中定位声音(例如,噪声)源,例如,用于在围绕物中定位声源。
将理解,预先计算的滤波器参数可通过计算机程序生成,该计算机程序可与用于执行阵列聚焦(例如,波束形成)处理的计算机程序分开或者包含在该计算机程序内(例如,作为该计算机程序的一部分或位于一个或多个单独的文件内)或者作为其组合。例如,预先计算的滤波器参数可通过计算机程序或安装程序生成并且存储在数据库或利用传感器阵列递送的其他存储介质内。
一种计算机程序可包括程序代码方法,程序代码方法被适配为:当在数据处理系统上执行程序代码方法时,使数据处理系统执行在上面和下面公开的方法的步骤。计算机程序可存储在计算机可读存储介质上或者体现为数据信号。存储介质可包括用于存储数据的任何合适的电路或装置,诸如,RAM、ROM、EPROM、EEPROM、闪存、诸如CD ROM、DVD、硬盘等的磁或光学存储装置。
一种计算机可读存储介质具有存储在其上的一组预先计算的滤波器参数,每个滤波器参数与一组网格点中的一个相关联,并且当执行本文公开的方法的实施方式的步骤时,每个滤波器参数被本文定义的处理设备使用。
根据另一个一般方面,本文中公开的是一种确定包括多个传感器的传感器阵列的聚焦输出信号的方法,每个传感器可操作为输出响应于被测量的量的传感器输出信号,聚焦输出信号表示在所计算的在聚焦点处的量。该方法的实施方式包括:
-从每个传感器中接收相应测量的传感器输出信号;
-从接收的被测量的传感器输出信号的每一个中计算各个被滤波的传感器信号,被滤波的传感器信号表示由与相应的传感器相关联的滤波器滤波的被测量的传感器信号,其中,从相应的传感器中接收被测量的传感器信号,滤波器取决于依赖聚焦点的至少一个滤波器参数;
-组合多个所计算的滤波传感器信号,以获得聚焦输出信号。
附图说明
从下面参照附图描述的实施方式中,以上和其他方面将变得显而易并且被阐明,其中:
图1示出了波束形成器系统的示意性方框图;
图2示出了计算传感器阵列的波束形成的输出信号的处理的流程图;
图3示出了在平面阵列中具有麦克风的测量设置。
贯穿附图,在可行的情况下,相同的附图标记表示相同或相应的元件、特征或部件。
具体实施方式
在下文中,将参照波束形成系统更加详细地描述本发明的方面和实施方式。然而,将理解,本文中描述的方法、产品和系统的实施方式还可适用于声全息技术。
图1示出了用于执行声波的波束形成的波束形成器系统的示意性方框图。该系统包括一组声学接收器108和连接至声学接收器的分析单元103。
在下文中,声学接收器108将还被称为换能器。然而,将理解,声学接收器可为麦克风、水诊器或用于测量声学特性(诸如,声压、声压梯度、粒子速度或其他线性量)的任何其他合适的传感器。在图1的实施方式中,换能器108被实施为其中换能器108被布置成规则的网格(例如,一维、二维或三维网格)的换能器的阵列102。换能器可被布置在规则的网格或不规则的几何体上。在波束形成应用中,不规则的几何体表现得更好,并且因此被通常使用。可根据围绕物(在该围绕物内,定位声源)的尺寸和几何形状、感兴趣的频率范围、期望的空间分辨率和/或其他设计参数来选择换能器的数量和阵列的几何形状(例如,换能器间的间距)。
换能器阵列102连接至分析单元103,以便换能器108可将测量的信号例如经由线或无线信号连接转发给分析单元。由换能器测量的信号也被称为传感器信号。
分析单元103包括用于接收和处理来自换能器阵列102的所测量的信号的接口电路104、与接口电路104进行数据通信的处理单元105、存储介质112以及与处理单元105进行数据通信的输出单元106。即使在图1中示为单个单元,将理解,分析单元103也可物理地分成两个单独的装置(例如,采集前端和计算机)或者甚至两个以上的装置。同样,将理解,相对于分析单元的不同子块描述的功能可分成可替代的或另外的功能或硬件单元/模块。
接口电路包括信号处理电路,信号处理电路适合于从换能器108中接收输出信号并且处理所接收的信号以用于处理单元105的后续分析。接口单元执行同步时间数据采集,并且随后可通过处理单元105来进行所有进一步的处理,包括通常使用FFT数据至频域的转换。接口电路104可包括一个或多个以下部件:用于放大所接收的信号的一个或多个前置放大器、用于将所接收的信号转换成数字信号的一个或多个模数(A/D)转换器、一个或多个滤波器(例如,带宽滤波器)等。例如,接口电路可将作为每个换能器的频率函数的幅度和相位提供为输出数据。
处理单元105可为适当编程的微处理器、计算机的中央处理单元或者用于处理从接口电路104接收的信号的任何其他合适的装置(例如,ASIC、DSP、GPU等)。处理单元被适配于处理经由接口电路104接收的传感器信号,以便计算如本文所述的波束形成的输出信号。
存储介质112可包括用于存储表示一组预先计算的滤波器参数的数据的任何合适的电路或装置,诸如,RAM、ROM、EPROM、EEPROM、闪存、磁或光学存储装置(诸如CD ROM、DVD、硬盘等)。在图1中,存储介质被示为单独的但与处理单元通信连接。然而,将理解,存储介质112还可体现为处理单元105的一部分(例如,内部存储器)。
输出单元106可包括显示器或用于提供波束形成的信号的视觉表示(例如,关于不同聚焦点的波束形成的输出信号的图)的任何其他合适的装置或电路。合适的输出单元的实例包括打印机和/或打印机接口(用于提供打印的表示)。可替代地/另外,输出单元106可包括用于传递和/或存储表示波束形成的输出信号的数据的任何合适的电路或装置,例如,RAM、ROM、EPROM、EEPROM、闪存、磁或光学存储装置(诸如,CD ROM、DVD、硬盘)、有线或无线数据通信接口(例如,至计算机或诸如LAN、广域网和互联网的电信网络的接口)等。
分析单元103可被实施为适当编程的计算机,例如,包括合适的信号采集板或电路的PC。
在操作期间,换能器阵列102可位于在定位声源的周围的位置处或者绘制声源的位置处,例如,在包括发射声学辐射的声源的物体的表面附近或者在围绕物内部。可根据被分析的物体或环境的尺寸和几何复杂性、感兴趣的频率范围和期望的空间分辨率来选择换能器的数量、阵列的几何形状(例如,换能器间的间距)以及至可能的声源的距离。
阵列102的位置可通过例如位置检测装置来确定并且被馈送至分析单元103内。阵列102的换能器108可在其相应的位置处测量声压或另一种合适的声学量,并且所产生的传感器信号被发送至分析单元103。
例如,换能器阵列可为具有集成的位置检测装置的手持式阵列,从而允许在分布在物体周围的不同的可到达的位置处进行测量。另一种典型的应用可以是在车厢内部,其中,3D阵列网格可被用于能够识别在所有方向中的源,例如,可使用球面阵列或双层阵列(例如,包括8x8x2个传感器)。
分析单元103从通过换能器测量的信号中计算关于一个或多个聚焦点109和/或方向的波束形成的输出信号。分析单元可存储和/或输出波束形成的信号的表示,例如,作为方向和/或聚焦点的函数的在阵列位置处的声音强度或对声压的贡献的图。
将参照图2并继续参照图1来描述用于计算波束形成的输出信号的处理的实施方式。
通常,该处理的实施方式可计算在给定位置/方向r处的输出时间信号b(t,r)。例如,输出信号可为聚焦点/方向对传感器阵列(其中心)处的声压的贡献的估计量。如上所述,r可定义在3D中的位置或方向(例如,由在传感器阵列的中心或者另一个合适的坐标系的原点的周围的单位球面上的位置表示的)。FAS波束形成器可首先将单独的滤波器h应用于传感器信号pm(t),m=1,2,...,M,以获得被滤波的信号,然后将被滤波的信号相加:
b ( t , r ) = Σ m p m ( t ) ⊗ h ( t , α m ( r ) ) - - - ( 1 )
此处,符号表示时间卷积,并且向量αm(r)包含应用于换能器编号m的滤波器的参数,以帮助在位置r处聚焦波束形成器。因此,通过组合(在该实施方式中,为求和)多个所计算的被滤波的传感器信号,来获得波束形成的输出信号b(t,r)。通常,将理解,可通过对所计算的被滤波的传感器信号的不同的组合(尤其是线性组合)来计算波束形成的输出信号。滤波器可为FIR滤波器;然而,本文中描述的处理的实施方式还可应用其他类型的滤波器。S.Yan、C.Hou和X.Ma的文章“Convex optimizationbased time-domain broadband beamforming with sidelobe control”(J.Acoust.Soc.Am.121(1),2007,1,46-49)描述了基于FIR滤波器的方法(包括例如为了最小化波束形成器的旁瓣电平而优化FIR滤波器的方法)。技术人员将理解,在很多情况下,最佳滤波器参数向量αm(r)的计算是非常繁重的计算任务,解决该任务不适于在每次波束形成计算时必须为每个计算点执行。在下文中,将更加详细地描述执行最佳FAS波束形成的有效方法的实例。
在初始步骤S1中,该处理定义了合适的坐标系,并且获得形成网格110的一组N个网格点rn,n=1,2,...,N,可为该网格获得预先计算的最佳滤波器参数。位置rn的网格跨越其中利用特定的阵列执行波束形成计算的空间内的区域。通常,在网格的点rn之间的间距应略微小于波束形成器的本地波束宽度。该间距可取决于所使用的插值方案并且可在设计整个系统(滤波器、插值方案、网格生成)时确定。此时,间距可被选择地足够小以便获得期望的插值精度。该处理为这些网格点rn的每一个获得一组预先计算的最佳滤波器参数向量αm(rn),其中,n=1,2,...,N是与N个网格点相关联的指数,并且m=1,2,...,M是与该阵列的M个换能器相关的指数。滤波器参数αm(rn)因此定义了与各个网格点相关联的各个辅助滤波器。例如,网格点和相关联的滤波器参数可存储在滤波器库(例如,文件、一组文件或数据库)中,该滤波器库与该阵列相关联并且可与该阵列一起设置。下面将更加详细地描述用于计算最佳滤波器参数的方法的实例。通常,网格和滤波器参数通过特定的/单独的程序一次性定义/计算并且与传感器阵列一起递送。但是还可在安装期间或者在波束形成的软件内通过合适的初始化函数来进行(从而允许重新计算)。
在步骤S2中,该处理从阵列的各个换能器中接收测量,即,一组传感器信号pm(t),m=1,2,...,M。
在步骤S3中,该处理选择定义位置或方向(将计算关于其的波束形成的信号)的向量r。例如,该处理可自动选择一系列位置,或者该处理可接收表示期望的方向/位置的用户输入,或者用户可定义作为生成源分布图的基础的点r的网格。
在步骤S4中,该处理识别最靠近聚焦点r的周围最近的网格点。在图1中,最靠近点109的网格点表示为113。将理解,所识别的网格点的数量可取决于所使用的插值算法、坐标系的选择和网格的类型。例如,如果网格点被设置在立方体网格中,则该处理可将周围最近的网格点确定为聚焦点r位于其中的立方体的拐角。同样,如果网格点被设置在以坐标系的原点为中心的并且具有相应半径的一个或多个球面上,则该处理可确定最靠近位置r的球面,并且随后确定在每个球面上的最接近的网格点的预定数量。
在步骤S5中,该处理使用与各个网格点相关联的预先计算的滤波器参数来为每个确定的最接近的网格点执行波束形成计算:
b ( t , r n ) = Σ m p m ( t ) ⊗ h ( t , α m ( r n ) )
因此,通过组合(在该示例下,为求和)将辅助滤波器应用于每个传感器信号而产生的辅助滤波传感器信号来计算波束形成器输出。辅助滤波器反过来由预先计算的滤波器参数αm(rn)来定义。在步骤S6中,该处理执行对由此计算的在位置rn处的波束形成器输出的插值,以达到在位置r处插值的波束形成器输出b(t,r)。因此,通过下列步骤计算插值的波束形成器输出:
-通过将各个辅助滤波器应用于所测量的传感器输出信号,来为每个传感器计算与所选择的网格点子组相关联的多个辅助滤波传感器信号,每个辅助滤波器与相应的传感器相关,其中,从相应的传感器中接收所测量的传感器信号,并且每个辅助滤波器取决于至少一个预先计算的滤波器参数;
-为每个网格点组合为各个传感器计算的多个辅助滤波传感器信号,以获得辅助波束形成的输出信号b(t,rn),并且对为各个网格点计算的辅助波束形成的输出信号进行插值,以获得插值的波束形成的输出信号,
可使用已知的任何合适的插值技术来进行插值。简单的插值形式是在角度或线性坐标(1D、2D或3D)中的线性插值。
可通过以下方式,加速波束形成处理:
1)可首先横跨在获得规定的波束形成图中所涉及的所有计算点r上的插值所需要的至少该部分网格点执行波束形成计算。随后,在第二步骤中进行插值。
2)对于FIR滤波器的情况,等式(1)的波束形成器输出是滤波器系数hm(rn)(其现在构成参数向量αm(rn)的主要部分)的线性函数:
b ( t , r n ) = Σ m Σ l = 0 L - 1 p m ( t - ( l + v m ) · T SF ) · h m , l ( r n ) . - - - ( 2 )
其中,TSF是FIR滤波器的采样时间间隔,L是在滤波器中的标号(tab,标签)的数量,vm是在换能器m与所有网格点rn之间共享的整数采样间隔延迟,并且对于给定的n和m,hm,l(rn)构成滤波器系数向量hm(rn)的L个元素。代替对网格点波束形成器输出信号b(t,rn)进行空间插值以获得b(t,r),可插值FIR滤波器系数hm,l(rn)以获得hm,l(r),hm,l(r)随后可被应用在等式(2)中。因此,通过下列步骤计算波束形成的输出信号:
-为每个传感器从与所确定的网格点子组的相应一个相关联的预先计算的滤波器参数hm,l(rn)中,来计算至少一个插值的滤波器参数hm,l(r);
-通过将各个插值的滤波器应用于所测量的传感器输出信号,来为每个传感器计算插值的滤波传感器信号;以及
-组合为各个传感器计算的插值的滤波传感器信号,以获得插值的波束形成的输出信号。
在等式(2)中引入共享延迟量vm,以最小化所需的FIR滤波器长度:假设与DAS波束形成相似,在DAS中需要的延迟接近vm,因此FIR滤波器需要将最大延迟建模为几个样本间隔加上相对由滤波器优化引入的纯延迟的偏离。如果FAS波束形成软件不知道延迟偏移量vm,则将延迟偏移量与FIR滤波器系数的向量hm(rn)一起存储在滤波器库内(即,在αm(rn)中)。
利用与传感器信号相同的采样率计算等式(2)的波束形成器输出信号,使得公式成为对离散卷积的换能器求和:
b ( κ T SP , r n ) = Σ m Σ l = 0 L - 1 p m ( ( κ - l - v m ) · T SF ) · h m , l ( r n ) , - - - ( 3 )
K是对输出的整数采样计数。可使用例如利用FFT的标准重叠相加算法来有效地进行卷积的实际算法。
通过等式(2)的傅里叶变换可获得频域公式:
B ( ω , r n ) = Σ m P m ( ω ) e - jω v m T SF Σ l = 0 L - 1 h m , l ( r n ) e - jωl T SF , - - - ( 4 )
Pm(ω)是换能器声压频谱,并且ω是角频率。引入滤波器的频域表示Hm(ω,rn):
H m ( ω , r n ) ≡ Σ l = 0 L - 1 h m , l ( r n ) e - jωl T SF - - - ( 5 )
允许将频域波束形成表达式(4)写为:
B ( ω , r n ) = Σ m P m ( ω ) H m ( ω , r n ) e - jω v m T SF . - - - ( 6 )
为了简单地进行本描述,仅仅考虑基数2FFT处理,如果对于某个非负整数μ,传感器信号TS的采样时间间隔等于2μTSF,则可使用FFT来计算等式(5)的频域滤波器Hm(ω,rn)。利用传感器信号的记录长度K,在传感器信号的FFT频谱中的频率ωk=2πfK将为fK=k/(KTS),k=0,1,2,…K/2,并且等式(5)变成:
H m ( ω k , r n ) ≡ Σ l = 0 L - 1 h m , l ( r n ) e - j 2 πlk T SF / KT S = Σ l = 0 L - 1 h m , l ( r n ) e - j 2 π lk 2 μ K , - - - ( 7 )
可通过将从长度L到长度2μK的FIR滤波器系数向量进行补零随后进行FFT来为每个组合(n,m)计算该等式。在所测量的信号的采样频率等于或小于在FIR滤波器库内使用的采样频率的大部分情况下,这个可能支持非常有效的频域波束形成。
下面将说明替代性的频域实现方式。尤其地,在下文中,Pm表示由传感器m在给定频率处测量的复杂的传感器信号(例如,声压),并且Hm,n表示应用于来自传感器m的传感器信号的用于在网格点n上聚焦波束形成器的滤波器的复杂频率响应函数(FRF)。然后,在给定频率处的在网格点n上的波束形成的信号Bn
B n = Σ m H m , n P m
而且,在wn表示被应用于在网格点n上的波束形成的结果Bn的插值权重因子以便获得在期望的聚焦点上的波束形成信号B时,在给定频率处的插值的波束形成信号可获得为:
B = Σ n w n B n = Σ n [ Σ n w n H m , n ] P m
因此,为了计算B,首先可计算插值的频率响应,并且随后在波束形成计算中使用插值的频率响应函数。因此,在某些实施方式中,与每个网格点相关联的预先计算的滤波器参数可以是用于与各个传感器相关联的该网格点和各个频率的频率响应。
因为通常仅可获得在一组离散频率处的频率响应,可执行频率的额外插值,以便计算在任意的频率处的波束形成的结果。
针对最小化旁瓣电平的FIR滤波器的优化:
在下文中,将更加详细地描述在FIR滤波器并且在针对最小化旁瓣电平而进行优化的情况下预先计算最佳滤波器系数Hm(rn)的替代性方法。
为了该描述的目的,将描述用于单个聚焦点或方向(例如,网格点rn的其中一个)的一组最佳FIR滤波器系数的计算。将理解,下面描述的方法可用于为所有网格点rn计算滤波器参数。所计算的滤波器参数随后可被存储并用于如本文所描述的方法中。
图3示出了在平面阵列中具有M个换能器的测量设置。l+1个点源从位置ri,i=0,1,2,…l产生入射波,这些位置在此处相对于阵列位于远场区域内。使用FAS波束形成的聚焦能力,在阵列中应提取来自在参考位置处的源i=0的自由场(即,未放置阵列)声压贡献S0(t),同时在最大程度上抑制来自其他“干扰”源的贡献Si(t),i=1,2,…l。基于如从阵列中所看的声波的不同起点(位置ri)处进行提取。
尽管在图3的实例中示出了在自由场中的平面阵列,但是在本部分描述的原理适用于任何阵列几何图形和任何阵列支撑结构,包括嵌入式安装在板内或在刚性球面上。通常,可通过以下方式数学地表示传感器信号:
p m ( t ) = Σ i = 0 I s i ( t ) ⊗ g i , m ( τ ) , - - - ( 8 )
利用已知脉冲响应函数gi,m(τ),从自由场压力贡献Si(t)得到由换能器测量的对压力的实际贡献Pm(t)。对于在自由场中的换能器阵列,其中,可忽略换能器对声场的影响,并且对于仅在远场区域中的源,函数gi,m(τ)恰好表示从参考点到换能器m的声场#i的延迟τi,m
等式(8)的两侧的傅里叶变换导致等效的频域关系:
P m ( ω ) = Σ i = 0 I s i ( ω ) G i , m ( ω ) , - - - ( 9 )
其中,P、S和G分别是p、s和g的傅里叶变换。在gi,m(τ)表示延迟的以上实例中,相应的函数Gi,m(ω)将表示等价的相移(针对在自由场中的阵列和在远场中的源):
G i , m ( ω ) = e - jω τ i , m - - - ( 10 )
为了从源#0中提取贡献S0(t)或S0(ω),我们在该源极的位置r0处聚焦波束形成器。将频域表达式(4)用于波束形成器,由于仅考虑聚焦点r0,所以跳过指数n。用等式(4)中的换能器压力谱替换表达式(9)随后导致:
B ( ω ) = Σ i = 0 I B i ( ω ) , - - - ( 11 )
其中,Bi是来自源编号i的贡献:
B i ( ω ) ≡ Σ m S i ( ω ) G i , m ( ω ) e - jω v m T SF Σ l = 0 L - 1 h m , l e - jωl T SF = S i ( ω ) Σ m Σ l = 0 L - 1 G i , m ( ω ) e - jω ( l + v m ) T SF h m , l = S i ( ω ) u T ( ω , r i ) h , - - - ( 12 )
如下定义列向量u和h:
u ( ω , r i ) ≡ { G i , m ( ω ) e - jω ( l + v m ) T SF } m = 1,2 , . . . , M l = 0,1 , . . . , L - 1 h ≡ { h m , l } m = 1,2 , . . . , M l = 0,1 , . . . , L - 1 , - - - ( 13 )
元素同样地构成并且T表示向量或矩阵的转置。理想地,人们希望找到滤波器向量h,从而横跨感兴趣的频率范围,B0(ω)=S0(ω)并且Bi(ω)=0,i=1,2,…l,根据等式(12),该等式等价于uT(ω,r0)h=1并且uT(ω,ri)h=0,i=1,2,…l。
将描述的用于预先计算最佳滤波器参数的方法的第一实施方式与由S.Yan、C.Hou和X.Ma(同前)描述的方法相似。在该方法中,考虑抑制较小数量的干扰源。因此,将最小化来自相对较少的干扰源的假定位置的贡献。通过最小化波束形成器输出信号的功率,同时要求在输出信号内完全保留来自聚焦点r0的贡献,在平均意义上最小化了来自其他位置(方向)的干扰。
为了最小化输出功率,通过S.Yan、C.Hou和X.Ma(同前)示出该输出功率等于hTRh*,其中,R是传感器信号的协方差矩阵,h是要确定的FIR滤波器系数的向量,并且*表示复共轭。引入R的Cholesky分解R=UHU,其中,H表示共轭(Hermitian)转置,输出功率可表示为其可通过最小化||U*h||2来最小化
参照等式(12),现在,可如下在数学上描述该方法:
满足||U*h||2≤γ
|uTk,r0)h|=1,k=1,2,...,K    (14)
|uTk,rk,i)h≤ε,k=1,2,...,K,i=1,2,...Ik
||h||≤Δ
此处,ε是定义对来自所选择的I个位置的所要求的干扰抑制水平的常数,ωk是施加约束条件的一组频率,rk,i表示应抑制干扰源的一组频率可能变化的点,并且Δ定义FIR系数向量h的范数(通常是2-范数)的最大值,从而限制波束形成器的所谓的白噪声增益(WNG)(见S.Yan、C.Hou和X.Ma(同前))。使用一组依赖频率的点rk,i的原因在于波束宽度将趋于在低频率处很大并且在高频率处很小。问题(14)具有所谓的二阶锥形规划法(SOCP)问题的形式,可使用可用的求解器来有效地解决这个问题,例如见S.Yan、C.Hou和X.Ma(同前)。
与本实施方式相关的一个问题是协方差矩阵R的定义和计算。理想地,它应使用用于波束形成的传感器信号来计算,这意味着必须为每个新测量的每个聚焦点r0解决优化问题(14)。这是相当繁重的任务。例如,对于50个换能器和64个FIR滤波器标号,必须为每个聚焦点确定在向量h中的3.200个变量。为了避免这个问题,可使用在本文中描述的方法的实施方式(其使用与空间插值相结合的预先计算的滤波器的滤波器库)。但是随后,对于预期的应用,协方差矩阵必须是某个假定的平均值。通常,可使用在对角线上具有相同元素的对角线矩阵,这意味着||U*h||2的最小化可被||h||2的最小化代替,这恰好表示WNG的最小化。随后,可去除最后的约束条件-WNG约束条件。
然而,一个主要问题是使用相对较少的干扰源的的假定位置ri。尽管如上所述,通常对在其他位置的干扰源的灵敏度将被最小化,但是在某些频率处,在点ri之间仍可能具有高旁瓣(灵敏度)峰值。相关的问题是选择允许的最大相对旁瓣电平ε。原则上,在存在其他约束条件时,人们不会知道在所选择的频率ωK上可实现哪个旁瓣电平。这需要与给定的阵列和给定的聚焦点相关的某些研究。
通过跳过||h||2的最小化,并且通过对最大的旁瓣电平ε进行最小化,以下替代性实施方式解决了这些问题。尤其地,基于以上讨论,可定义以下优化问题以确定FIR滤波器系数:
满足|uTk,r0)h|=1,k=1,2,...,K    (15)
|uTk,rk,i)h|≤ε,k=1,2,...,K,i=1,2,...Ik
||h||≤Δ
这也是SOCP问题,因此,可使用可利用的求解器来有效地解决该问题。依然可通过最后的约束条件来控制WNG,在位置rk,i之间的旁瓣电平上提供某种程度的“平均控制”。但是在其他方面,建议使用足够密集的点rk,i的网格以将在期望范围内的所有旁瓣控制到某种精确的程度并且从而最小化最高的旁瓣。
在此处的一个主要问题是优化问题(15)的大小,可通过以下实例说明该问题。
实例:
在一个实例中,使用嵌入式安装在半径α=10cm的刚性球面中的M=50个麦克风的阵列。随后,在麦克风之间的平均间距接近5cm,这刚好小于波长的一半到接近2.5kHz。在该频率范围内,可利用50个麦克风在球面上的精选分布来实现低旁瓣电平。因此,在该实例中,将Fmax=2,5kHz被选择为上限频率。
在2.5kHz处,具有ka≡ωa/c≈4.58,k在此处是波数,并且c为声音的传播速度。在球形阵列的方向性图案中,将自然地控制高达大约为该数量N≈ka的球谐度(见M.Park和B.Rafaely的“Sound-field analysis byplane-wave decomposition using spherical microphone array”,J.Acoust.Soc.Am.118(5),2005,11,2005,3094-3103)。然而,允许合理的WNG电平,即,解向量h的合理范数,高达大约为N≈最小值{ka+2,6}的球谐度在优化期间被放大为大约处于相同的电平。随后,波瓣的最小峰值-到-零角半径为θr=π/N,见M.Park和B.Rafaely(同前)。为了不过多地错过旁瓣峰值,已经发现在网格rk,i中的角距Δθ不应超过最小波瓣角半径θr的大约20%。选择具有覆盖整个4π立体角的角距为Δθ=0.2π/Ν的网格需要的点数大约为:
I ≈ 2 1 - cos ( Δθ ) - - - ( 16 )
在2.5kHz处,这导致I≈365个点。但是在这些当中,必须去除位于主瓣区域内的点(达到大约50点),因此,以315个点结束。
如果要建模比延迟和求和的分辨率略微好的低频分辨率,则使用等于16,384个采样/秒的采样频率,根据经验,通常需要FIR滤波器长度L=128。然后,滤波器的时间跨度接近波围绕球体4次所需要的时间。为了确定L个FIR滤波器系数,应使用数量K个(至少等于L/2)频率ωK。由于每个频率对具有实部和虚部的复杂响应有效地提供了约束条件,所以除以2。然而,人们发现某个多元决定(over-determination)有用,因此最后,选择K≈0.7L≈90。在从0Hz到2.5kHz的频率间隔上均等地分布90个频率,并且使用以上方案来计算每个频率处的约束位置rk,i的数量lk,以点rk,i的总数大约为13,600结束。在(15)中的这些大约13,600个不等式约束之外,我们使又一个不等式约束加上K=90个等式约束。在h中可变的FIR滤波器系数的数量(在解决处理中确定)如上所述:ML=50·128=6400。这已经不太适合于在每个波束形成处理期间进行。
但是,实际上通常使用的波束形成高达大约为“波长空间采样的一半的频率”的3倍的因子,在目前的情况下高达7.5kHz。这样做,需要使时间采样频率加倍,即,高达32,768个采样/秒,这就再次需要使FIR滤波器标号的数量L加倍,以便在短时间内保持滤波器长度。而且,在用于选择频数计数K的先前使用的标准之后,该计数也必须加倍,因此以L=256、K=180以及可变FIR系数的数量等于12,800的结束。然而,由于在2.5与7.5kHz之间的添加的频率处需要非常密集的网格rk,i,出现问题大小的最大增加。发现最大球谐度N利用频率进行的以下增加有效:
N &ap; ka + 2 ka &le; 4 6 4 < ka < 6 ka ka &GreaterEqual; 6 . - - - ( 17 )
在0Hz至7.5kHz之间均匀地分布K=180个频率,并且使用先前使用的用于计算在rk,i网格中的密度的方案,以总共大约116,000个不等式约束结束。考虑约束函数的梯度的矩阵,它将具有12,800x116,000=1,484,800,000个元素,在浮点表示中占用6GB的内存。在合理的时间内的波束形成计算中,认为不容易为每个聚焦点处理该问题大小。
从以上实例中得出的结论是,第二实施方式可用于在换能器计数、阵列直径和/或上限频率方面相对较小的问题。对于例如20cm直径的50个元素球形阵列的实际应用,该方法不可行。它可能由强大的计算机来处理以用于将根据本文中公开的方法的实施方式使用的一次性生成的滤波器库。
以下第三实施方式避免对用于预先计算最佳滤波器参数的大计算功率的需要。除此之外,它还解决了在“波长空间采样的一半的频率”之下和之上能够实现截然不同的旁瓣抑制电平的问题。由于在所有的频率处使用相同值ε,所以上述第二替代实施方式可能并不能始终能够适当地处理该问题。人们可尝试定义变化的频率ε,但是将难以了解是否已经定义了“良好的变化”。
根据第三实施方式,在用于源#i对波束形成器输出的贡献的表达式(12)中,我们使用FIR滤波器Hm(ω)的频域表示的定义(5)(再次,为了该描述的目的,跳过指数n,仅考虑单个聚焦点)以获得:
B i ( &omega; ) = S i ( &omega; ) &Sigma; m G i , m ( &omega; ) e - j&omega; v m T SF H m ( &omega; ) . - - - ( 18 )
与在(13)中的向量μ的定义对应,列向量v现在定义为:
v ( &omega; , r i ) &equiv; { G i , m ( &omega; ) e - j&omega; v m T SF } m = 1,2 , . . . , M . - - - ( 19 )
与在第一和第二实施方式中一样,第三实施方式包括最优化在离散的一组频率ωk,k=1,2,…K处的FIR波束形成器的性能。为了便于对本实施方式的矩阵向量进行描述,通过以下方式定义在这些离散频率处的FIR滤波器响应的矩阵H:
并且Hm和Hk分别被用于表示包含在H中的列m和行k的元素的列向量。使用这些定义,在频率为ωk处,单个源对波束形成器输出的贡献的表达式(18)可表示为:
Bik)=Sik)VTk,ri)Hk    (21)
显然,如果应在波束形成器输出中适当地重构源#i,则VTk,ri)Hk的值在所有频率处应等于1。如果应避免其贡献,则VTk,ri)Hk在所有频率处应尽可能较小。
该方法的第三实施方式通常通过以下两个步骤过程来优化FIR滤波器:
1)一次采用一个频率ωk,k=1,2,…K,计算在Hk内的该组复杂换能器权重作为问题的解决方案:
满足vTk,r0)Hk=1    (22)
|vTk,rk,i)Hk|≤εk,i=1,2,...Ik
||Hk-Hk (0)||≤Δ
在完成后,可获得全矩阵H。下面将解释Hk (0)的选择。
2)对于每个换能器m,使具有标号的FIR滤波器
hm≡{hm,l}l=0,1,...,L-1
与频率响应向量Hm相配。这样做,在频率ωk处应用等式(5),再次跳过聚焦点指数n,因为在此处仅考虑单个点(#0):
Hm=Ahm.    (23)
在此处,A是具有K行(频率)和L列(滤波器标号)的以下矩阵:
在标准最小二乘的意义上,求解等式(23)。
要注意的是,在问题(22)中的等式约束包括与在问题(14)中的等式约束相反的相位。这是因为Hmk)的一致“平稳”相位需要跨过频率ωk以在第二步骤中从等式(23)中获得有用的解决方案hm
令人惊讶地,已经证明了,当适当地选择在等式(22)中的Hk (0)时,实际上从上面步骤1)中产生了相位的这种平稳的频率变化。问题(22)是使用迭代优化处理解决的SOCP问题。代替将具有所有零的HK用作起始点,发明人使用满足在(22)中的等式约束的状态良好的分析起始向量Hk (0)。通常,这将是用于“透明”阵列的延迟相加的解决方案以及用于在(例如,所谓的嵌入式安装)刚性球体中的换能器阵列的球谐波束形成解决方案。在一个成功的实现方式中,使用最陡下降算法,以Hk (0)开始,并且当穿过超球面的路径||Hk-Hk (0)||=Δ的情况下,或者在梯度范数变成小于(例如1%)其在Hk (0)处的值时,停止迭代。Δ的典型值为Δ=||Hk (0)||。使用该过程,证明产生可非常容易由相当短的FIR滤波器hm表示的频率响应函数Hm
利用适当选择的频率ωk,关系(23)恰好表示在Hm与hm之间的FFT关系。然而,通过该方式选择频率,不支持多元决定。并且在任何情况下,在步骤2)中的(23)的解决方案是比在步骤1)中的优化问题小得多(消耗少得多的时间)的计算任务。
最大的计算问题是在最高频率处的(22)的解决方案。对于以上结合第二实施方式讨论的50个元素的球形阵列实例,180个频率中的每个涉及具有M=50个复杂变量的(22)的解决方案,并且在7.5kHz处,不等式约束的数量大约为1862。这是大小易于处理的问题。但是,回想必须处理180个频率并且必须为每个聚焦点这样做,它依然是不适合结合每个波束形成计算进行的计算问题。对于与在该实例中使用的阵列类型相结合的实际应用,本文中描述的滤波器库具有相当大的优势。
如上所述,在滤波器库内覆盖的点rn之间的间距应略微小于波束形成器的本地波束宽度。对于以上实例的50个元素的球形阵列(具有大约5°间距,仅加上3个半径距离便从球体超过25cm)适用于在角度(2D)中以及在径向距离中的线性插值。对于径向插值,发现最好将线性变化假设为反距离的函数而不是距离的函数
例如,当分析机器、马达、引擎、车辆(例如,汽车)等的声学特性时,本文中描述的方法和设备可用于识别各种声音/噪声源,诸如,振动物体。
可通过包括几个分立元件的硬件的手段和/或至少部分通过适当编程的微处理器的手段来实施本文中描述的方法的实施方式。在列举几种器件的产品权利要求中,这些器件中的几个可具体化为同一个部件、元件或硬件的条目。在彼此不同的从属权利要求中叙述或者在不同的实施方式中描述某些措施并非表示这些措施的组合这一起码的事实不能被用来获利。
应强调的是,在术语“包括/包含”用于说明书中时,被用来规定存在指定的特征、部件、步骤或元件,但是并不排除一个或多个其他特征、部件、步骤、元件或其组合的存在或增加。

Claims (18)

1.一种确定包括多个传感器的传感器阵列的聚焦输出信号的方法,每个传感器可操作为输出响应于被测量的量的传感器输出信号,所述聚焦输出信号表示所计算的在聚焦点处的量;所述方法包括:
-从每个所述传感器中接收相应测量的传感器输出信号;
-通过执行关于所测量的传感器信号的聚焦计算,来计算所述聚焦输出信号;
其中,所述方法进一步包括确定一组预定网格点的网格点子组,每个网格点具有与其相关的至少一个预先计算的滤波器参数;并且其中,计算所述聚焦输出信号包括执行关于所述网格点子组的插值,以便获得插值的聚焦输出信号。
2.根据权利要求1所述的方法,其中,所述聚焦输出信号表示多个被滤波的传感器信号的组合,每个被滤波的传感器信号表示通过滤波器滤波的所测量的传感器信号中的相应一个,所述滤波器与从中已接收相应的所测量的传感器信号的所述相应的传感器相关联,所述滤波器取决于依赖所述聚焦点的至少一个滤波器参数。
3.根据权利要求1或2所述的方法,其中,计算所述插值的聚焦输出信号包括将一个或多个辅助滤波器应用于相应的所测量的传感器输出信号中的每一个;所述一个或多个辅助滤波器中的每一个与从中已接收相应的所测量的传感器信号的所述相应的传感器相关联,并且所述一个或多个辅助滤波器中的每一个取决于至少一个所述预先计算的滤波器参数。
4.根据前述权利要求中任一项所述的方法,其中,所述插值是线性插值。
5.根据权利要求1至3中任一项所述的方法,其中,所述插值包括在角度坐标中的插值。
6.根据前述权利要求中任一项所述的方法,其中,计算所述聚焦输出信号包括:
-通过将相应的辅助滤波器应用于所测量的传感器输出信号来为每个传感器计算多个辅助滤波传感器信号,每个辅助滤波器与从中已接收所测量的传感器信号的所述相应的传感器相关联,并且每个辅助滤波器取决于至少一个所述预先计算的滤波器参数;
-为每个网格点组合为相应的传感器计算的所述多个辅助滤波传感器信号,以获得辅助聚焦输出信号,并且对为相应的网格点计算的辅助聚焦输出信号进行插值,以获得所述插值的聚焦输出信号。
7.根据权利要求1至5中任一项所述的方法,其中,计算所述聚焦输出信号包括:
-通过将相应的辅助滤波器应用于所测量的传感器输出信号来为每个传感器计算多个辅助滤波传感器信号,每个辅助滤波器与从中已接收所测量的传感器信号的所述相应的传感器相关联,并且每个辅助滤波器取决于至少一个所述预先计算的滤波器参数;
-为每个传感器插值为相应的网格点计算的所述多个辅助滤波传感器信号,以获得插值的滤波传感器信号,并且组合为相应的传感器计算的所述插值的传感器信号,以获得所述插值的聚焦输出信号。
8.根据权利要求6或7所述的方法,包括为多个聚焦点确定相应的网格点子组;并且为所有确定的子组的网格点计算相应的辅助滤波传感器信号。
9.根据权利要求1至5中任一项所述的方法,其中,计算所述聚焦输出信号包括:
-为每个传感器从与所确定的网格点子组中的相应一个相关联的所述预先计算的滤波器参数中计算至少一个插值的滤波器参数;
-通过将相应的插值滤波器应用于所测量的传感器输出信号来为每个传感器计算插值的滤波传感器信号;以及
-组合为相应的传感器计算的所述插值的滤波传感器信号,以获得所述插值的聚焦输出信号。
10.根据前述权利要求中任一项所述的方法,其中,通过最小化最大旁瓣电平来计算与第一网格点相关联的每个预先计算的滤波器参数,所述最大旁瓣电平表示来自不同于所述第一网格点的一组其他位置的并且位于一组预定的频率处的干扰抑制水平。
11.根据权利要求1至9中任一项所述的方法,其中,通过以下步骤计算与第一网格点相关联的每个预先计算的滤波器参数:
-通过最小化最大旁瓣电平来为一组预定频率中的每一个确定一组传感器权重,所述最大旁瓣电平表示来自不同于所述第一网格点的一组其他位置的并且位于所述频率处的干扰抑制水平,每个传感器权重与相应的一个所述传感器相关联;
-通过将由所述预先计算的滤波器参数定义的频率响应与由所述确定的传感器权重的子组构成的频率响应向量相匹配,来为每个传感器确定至少一个预先计算的滤波器参数,所述子组的所述传感器权重与相应的频率相关联。
12.根据前述权利要求中任一项所述的方法,其中,所述网格点被布置为在相邻的网格点之间的距离小于波束形成器的本地波束宽度。
13.根据前述权利要求中任一项所述的方法,其中,所述聚焦计算是滤波和相加计算。
14.根据前述权利要求中任一项所述的方法,其中,确定所述网格点子组包括选择在所述聚焦点附近的网格点。
15.一种用于执行阵列聚焦的处理设备,所述处理设备包括:接口,用于从传感器阵列的相应的传感器中接收一组响应于被测量的量的传感器输出信号;以及处理单元,被配置为执行在前述权利要求的任一项中限定的所述方法的所述步骤,其中,所述处理单元包括用于存储所述一组预先计算的滤波器参数的存储介质。
16.一种声学处理系统,包括在权利要求15中限定的处理设备和一组传感器,所述一组传感器用于测量在一组测量位置处的被测量的量并且在通信连接中可连接至所述设备以将所述被测量的第一声学量转发至所述处理设备。
17.一种计算机程序,包括程序代码方法,所述程序代码方法被适配为:当在数据处理系统上执行所述程序代码方法时,使所述数据处理系统执行根据权利要求1至14中任一项所述的方法的所述步骤。
18.一种计算机可读介质,具有存储在其上的一组预先计算的滤波器参数,所述一组预先计算的滤波器参数的每一个与一组网格点的其中一个相关联,并且当执行根据权利要求1至14中任一项所限定的所述方法的所述步骤时,所述一组预先计算的滤波器参数的每一个被在权利要求15中限定的处理设备使用。
CN201280053952.9A 2011-11-04 2012-11-02 计算上有效的宽带滤波和相加阵列聚焦 Active CN103999151B (zh)

Applications Claiming Priority (5)

Application Number Priority Date Filing Date Title
US201161555597P 2011-11-04 2011-11-04
EP11187871.6 2011-11-04
EP11187871 2011-11-04
US61/555,597 2011-11-04
PCT/EP2012/071717 WO2013064628A1 (en) 2011-11-04 2012-11-02 Computationally efficient broadband filter-and-sum array focusing

Publications (2)

Publication Number Publication Date
CN103999151A true CN103999151A (zh) 2014-08-20
CN103999151B CN103999151B (zh) 2016-10-26

Family

ID=48191409

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201280053952.9A Active CN103999151B (zh) 2011-11-04 2012-11-02 计算上有效的宽带滤波和相加阵列聚焦

Country Status (7)

Country Link
US (1) US9299336B2 (zh)
EP (1) EP2774143B1 (zh)
JP (1) JP2015502524A (zh)
KR (1) KR101961261B1 (zh)
CN (1) CN103999151B (zh)
DK (1) DK2774143T3 (zh)
WO (1) WO2013064628A1 (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111581726A (zh) * 2020-05-11 2020-08-25 中国空气动力研究与发展中心 一种在线的一体化飞行器气动力建模系统
CN113139157A (zh) * 2021-04-22 2021-07-20 中山香山微波科技有限公司 一种dut主能量方向的计算方法和计算机设备

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9299336B2 (en) * 2011-11-04 2016-03-29 Brüel & Kjær Sound & Vibration Measurement A/S Computationally efficient broadband filter-and-sum array focusing
US9615172B2 (en) * 2012-10-04 2017-04-04 Siemens Aktiengesellschaft Broadband sensor location selection using convex optimization in very large scale arrays
FR3000862B1 (fr) * 2013-01-08 2015-01-09 ACB Engineering Dispositifs passifs d'acquisition acoustique large bande et systemes passifs d'imagerie acoustique large bande.
WO2015010850A2 (en) * 2013-07-22 2015-01-29 Brüel & Kjær Sound & Vibration Measurement A/S Wide-band acoustic holography
JP6106571B2 (ja) * 2013-10-16 2017-04-05 日本電信電話株式会社 音源位置推定装置、方法及びプログラム
US9326060B2 (en) * 2014-08-04 2016-04-26 Apple Inc. Beamforming in varying sound pressure level
US20160259050A1 (en) * 2015-03-05 2016-09-08 Navico Holding As Systems and associated methods for updating stored 3d sonar data
US9918321B1 (en) * 2017-03-28 2018-03-13 Sprint Communications Company L.P. Wireless communication system to optimize traffic management on a multi-band wireless base station
CN108735228B (zh) * 2017-04-20 2023-11-07 斯达克实验室公司 语音波束形成方法及系统
DE102017219235A1 (de) 2017-10-26 2019-05-02 Siemens Aktiengesellschaft Verfahren und System zum akustischen Überwachen einer Maschine
CN109507640A (zh) * 2018-12-18 2019-03-22 重庆大学 一种基于实心球阵列的全方位等效源声源识别方法
CN112462415B (zh) * 2020-11-02 2023-07-21 中国电子科技集团公司第三研究所 一种对多振动源进行定位的方法及装置
CN116148770B (zh) * 2023-04-21 2023-07-07 湖南工商大学 基于阵列信号处理的声源定位方法、装置及系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2192061A (en) * 1986-06-27 1987-12-31 Plessey Co Plc A phased array sonar system
CN1866356A (zh) * 2005-08-15 2006-11-22 华为技术有限公司 一种宽带波束形成方法和装置
CN101719368A (zh) * 2009-11-04 2010-06-02 中国科学院声学研究所 高声强定向声波发射的方法及其装置
CN102043145A (zh) * 2010-11-03 2011-05-04 中国科学院声学研究所 基于声矢量传感器均匀直线阵的快速宽带频域波束形成方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS5931466A (ja) * 1982-08-16 1984-02-20 Oki Electric Ind Co Ltd ソ−ナ−聴音装置
JP2886330B2 (ja) * 1990-11-16 1999-04-26 沖電気工業株式会社 信号方位測定装置
AU700274B2 (en) * 1995-06-29 1998-12-24 Teratech Corporation Portable ultrasound imaging system
EP1184676B1 (en) * 2000-09-02 2004-05-06 Nokia Corporation System and method for processing a signal being emitted from a target signal source into a noisy environment
KR100493172B1 (ko) * 2003-03-06 2005-06-02 삼성전자주식회사 마이크로폰 어레이 구조, 이를 이용한 일정한 지향성을갖는 빔 형성방법 및 장치와 음원방향 추정방법 및 장치
US7415117B2 (en) * 2004-03-02 2008-08-19 Microsoft Corporation System and method for beamforming using a microphone array
JP2005265466A (ja) * 2004-03-16 2005-09-29 Nittobo Acoustic Engineering Co Ltd ミニマックス規範に基づくビームフォーミング装置、及びそのビームフォーミング方法
US8345890B2 (en) * 2006-01-05 2013-01-01 Audience, Inc. System and method for utilizing inter-microphone level differences for speech enhancement
US9299336B2 (en) * 2011-11-04 2016-03-29 Brüel & Kjær Sound & Vibration Measurement A/S Computationally efficient broadband filter-and-sum array focusing

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2192061A (en) * 1986-06-27 1987-12-31 Plessey Co Plc A phased array sonar system
CN1866356A (zh) * 2005-08-15 2006-11-22 华为技术有限公司 一种宽带波束形成方法和装置
CN101719368A (zh) * 2009-11-04 2010-06-02 中国科学院声学研究所 高声强定向声波发射的方法及其装置
CN102043145A (zh) * 2010-11-03 2011-05-04 中国科学院声学研究所 基于声矢量传感器均匀直线阵的快速宽带频域波束形成方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
M. KAJALA等: "Filter-and-sum beamformer with adjustable filter characteristics", 《2001 IEEE INTERNATIONAL CONFERENCE ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111581726A (zh) * 2020-05-11 2020-08-25 中国空气动力研究与发展中心 一种在线的一体化飞行器气动力建模系统
CN111581726B (zh) * 2020-05-11 2023-07-28 中国空气动力研究与发展中心 一种在线的一体化飞行器气动力建模系统
CN113139157A (zh) * 2021-04-22 2021-07-20 中山香山微波科技有限公司 一种dut主能量方向的计算方法和计算机设备
CN113139157B (zh) * 2021-04-22 2022-03-08 中山香山微波科技有限公司 一种dut主能量方向的计算方法和计算机设备

Also Published As

Publication number Publication date
CN103999151B (zh) 2016-10-26
KR20140106545A (ko) 2014-09-03
KR101961261B1 (ko) 2019-03-22
EP2774143B1 (en) 2018-06-13
WO2013064628A1 (en) 2013-05-10
US20140313859A1 (en) 2014-10-23
JP2015502524A (ja) 2015-01-22
US9299336B2 (en) 2016-03-29
DK2774143T3 (en) 2018-08-06
EP2774143A1 (en) 2014-09-10

Similar Documents

Publication Publication Date Title
CN103999151A (zh) 计算上有效的宽带滤波和相加阵列聚焦
Zuo et al. Damage identification for plate-like structures using ultrasonic guided wave based on improved MUSIC method
CN102089633B (zh) 用于重建声学场的方法
EP2297557B1 (en) Reconstructing an acoustic field
CN101995574B (zh) 一种近场聚焦波束形成定位法
Engholm et al. Direction of arrival estimation of Lamb waves using circular arrays
Ambroziński et al. Efficient tool for designing 2D phased arrays in lamb waves imaging of isotropic structures
Baravelli et al. Double-channel, frequency-steered acoustic transducer with 2-D imaging capabilities
Engholm et al. Adaptive beamforming for array imaging of plate structures using lamb waves
Lay-Ekuakille et al. Optimizing and post processing of a smart beamformer for obstacle retrieval
CN101860779A (zh) 用于球面阵的时域宽带谐波域波束形成器及波束形成方法
CN105556260A (zh) 宽带声全息
Palmese et al. Three-dimensional acoustic imaging by chirp zeta transform digital beamforming
US8416642B2 (en) Signal processing apparatus and method for removing reflected wave generated by robot platform
Rathsam et al. Analysis of absorption in situ with a spherical microphone array
Movahed et al. Air ultrasonic signal localization with a beamforming microphone array
Touzé et al. Double-Capon and double-MUSICAL for arrival separation and observable estimation in an acoustic waveguide
Li et al. Ghost image suppression based on particle swarm optimization-MVDR in sound field reconstruction
Engholm et al. Using 2-D arrays for sensing multimodal lamb waves
Flynn Remote acoustic sensing of vibrating structures for structural health monitoring
Jung et al. Development of an asymmetric sensor array with beamforming
Chen et al. Method of spatially correlated wideband ambient noise simulation for underwater acoustic array
CN116358692A (zh) 入射声场和散射声场分离方法、系统、装置及其存储介质
Ambrozinski et al. 2D aperture synthesis for Lamb wave imaging using co-arrays
Sohn et al. Development of Noise Source Detection System using Array Microphone in Power Plant Equipment

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant