CN114333753A - 一种基于麦克风阵列的风机管道anc系统参考信号构造方法 - Google Patents
一种基于麦克风阵列的风机管道anc系统参考信号构造方法 Download PDFInfo
- Publication number
- CN114333753A CN114333753A CN202111610076.2A CN202111610076A CN114333753A CN 114333753 A CN114333753 A CN 114333753A CN 202111610076 A CN202111610076 A CN 202111610076A CN 114333753 A CN114333753 A CN 114333753A
- Authority
- CN
- China
- Prior art keywords
- array
- data
- noise
- frequency
- matrix
- 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.)
- Pending
Links
- 238000010276 construction Methods 0.000 title claims description 22
- 238000000034 method Methods 0.000 claims abstract description 40
- 239000011159 matrix material Substances 0.000 claims abstract description 39
- 238000011068 loading method Methods 0.000 claims abstract description 12
- 238000009499 grossing Methods 0.000 claims abstract description 11
- 230000009466 transformation Effects 0.000 claims abstract description 11
- 238000007781 pre-processing Methods 0.000 claims abstract description 8
- 239000013598 vector Substances 0.000 claims description 33
- 230000008569 process Effects 0.000 claims description 11
- 238000005070 sampling Methods 0.000 claims description 9
- 238000001914 filtration Methods 0.000 claims description 7
- 230000005059 dormancy Effects 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 5
- 230000009471 action Effects 0.000 claims description 4
- 238000013461 design Methods 0.000 claims description 4
- 230000006870 function Effects 0.000 claims description 4
- 230000004044 response Effects 0.000 claims description 4
- 238000012935 Averaging Methods 0.000 claims description 3
- 238000003491 array Methods 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims description 2
- 238000003379 elimination reaction Methods 0.000 claims description 2
- 206010057239 Post laminectomy syndrome Diseases 0.000 claims 3
- 230000001629 suppression Effects 0.000 abstract description 8
- 238000000605 extraction Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 8
- 230000000694 effects Effects 0.000 description 7
- 238000001228 spectrum Methods 0.000 description 7
- 230000009467 reduction Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 230000003044 adaptive effect Effects 0.000 description 3
- 230000002411 adverse Effects 0.000 description 3
- 230000001427 coherent effect Effects 0.000 description 3
- 230000000737 periodic effect Effects 0.000 description 3
- 238000013459 approach Methods 0.000 description 2
- 238000000265 homogenisation Methods 0.000 description 2
- 238000009434 installation Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 230000036961 partial effect Effects 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 238000004378 air conditioning Methods 0.000 description 1
- 230000004888 barrier function Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000033228 biological regulation Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000005111 flow chemistry technique Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000012887 quadratic function Methods 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000005236 sound signal Effects 0.000 description 1
Images
Landscapes
- Circuit For Audible Band Transducer (AREA)
Abstract
一种基于麦克风阵列的风机管道ANC系统参考信号构造方法。以麦克风阵列设备作为前馈式ANC系统的参考信号传感器，算法基于FFT和二次插值法估计管道风机噪声中的旋转噪声频率，在完成旋转噪声分量提取等阵列数据预处理后，对于均匀环形阵列，首先进行圆谐波变换，对于均匀线性阵列或变换所得虚拟线性阵列，进行子阵划分并执行前后向空间平滑，对平滑后的阵列数据协方差矩阵进行对角加载，进而执行MPDR或LCMP波束成形，最终基于最小二乘法对波束成形所得的数据序列进行振幅和相位的估计，配合定时器构造出高质量风机旋转噪声参考信号。本发明的方法提高了FxLMS等前馈式ANC算法的稳定性和风机旋转噪声的抑制效果。
Description
技术领域
本发明属于主动噪声控制(Active Noise Control,ANC)领域，涉及一种基于麦克风阵列的风机管道ANC系统参考信号构造方法。
技术背景
在噪声控制技术中，主要有“被动噪声控制”和“主动噪声控制”两种手段。其中被动噪声控制主要采用特殊吸声材料包裹或填充噪声传播路径来完成噪声的控制，此方法主要用于高频噪声的抑制。对于低频噪声，由于声波波长较长，被动噪声控制对于材料性能和尺寸的要求较高，往往难以实现，且在很多应用场景下，被动噪声控制设备的安装将受到限制。而主动噪声控制与之相比，具有设备简单、安装便捷及对于低频噪声有良好抑制效果的特点，正逐渐成为噪声控制的主流技术。
随着社会经济的进步，人们愈发追求学习、工作环境中良好的听觉体验。不难想象，在写字楼、教学楼和图书馆等环境中，若中央空调等风机设备在满足出风量的同时，发出的噪声能得到有效控制，则此时的环境将更有利于人们专注地工作和学习。风机管道噪声控制这一经典的课题便是在这样的大背景下越来越受到重视。
一些研究风机噪声的文献，如武汉理工大学的伍先俊等人于2001年发表的期刊论文《风机叶片噪声机理及降噪》及郑冬明于1998年发表的《风机噪声的控制》等，共同论述了风机噪声产生的机理和噪声特点，并指出风机噪声主要由中低频为主的旋转噪声(周期性噪声)和中高频为主的湍流噪声所组成，前者在风机高转速低负荷的情况下构成了风机噪声的主要元素。对于湍流噪声，可用被动噪声控制的手段进行抑制，而对于旋转噪声，则需要引入ANC技术来进行有效抑制。
在ANC技术中，常见的控制系统分为前馈型和反馈型结构，二者的主要区别在于反馈型ANC系统仅有误差信号参与自适应滤波器系数的调节，而前馈型ANC系统除了误差反馈信号外，还有初级噪声源的前馈参考信号共同参与自适应滤波过程，因此前馈型系统通常具有更优的系统稳定性和降噪性能。前馈式ANC主要采用FxLMS、FxNLMS及FxRLS等经典算法，其中ANC控制器为兼顾算法逐点自适应这一严苛的实时性要求，通常采用专用于实时音频信号处理的DSP或FPGA设备，一般需要两路音频输入接口和一路音频输出接口。
在前馈型ANC系统中，参考信号主要存在次级扬声器到参考麦克风的声反馈、声波多径传播的声反射两种高度相干的干扰及其他方向上的如人类语音等非相干干扰，将影响系统的稳定性与降噪性能。因此如何获取高质量的参考信号成为前馈型ANC系统的研究重点。
一类规避参考信号干扰方法如中国发明专利申请，公开号CN113257214A所述，采用转速计这一非声学传感器来替代参考麦克风，获取风机旋转噪声的基波频率用以构造参考信号。
对比此类系统，本发明原始参考信号的采集利用麦克风这一声学传感器完成，结合相应算法，使得构造的信号可以保留旋转噪声的振幅信息和相位信息，此二者事实上构成FxLMS等ANC算法的主要调节对象。
在参考传感器采用声学传感器时，常见的相关硬件部分如哈尔滨工业大学的于梦娇于2013年发表的硕士学位论文《管道低频噪声主动控制系统的DSP实现》所述，参考麦克风往往仅由单只麦克风构成，且初级噪声往往由ANC控制器内部生成，参考信号的频率等参量被视为已知量。
对比此类系统，本发明以参考麦克风阵列设备替代单只麦克风，同时结合相应估计算法，无需获取信号频率和相位的先验信息。
对于参考信号的声反馈干扰问题，中原大学Kuan-Chun Chen等人在2017年发表的论文《Active noise control in a duct to cancel broadband noise》中给出了较直接的声反馈抑制方案，此方案需要对声反馈通道进行建模，声反馈抑制的效果与建模精度密切相关。
对比此类系统，本发明无需引入声反馈通道的建模及声反馈抑制算法，且将参考麦克风阵列设备和参考信号构造算法独立于ANC控制器设备，参考麦克风阵列设备可用低成本的Linux设备实现，降低了ANC控制器的运算压力。
还有一类规避声反馈的方案是采用心形麦克风，以其物理特性带来的空间指向性定向提取初级噪声源方向的信号，同时抑制次级扬声器方向的信号，然而此类麦克风一般成本不菲，多见于专业的录音室和舞台，且其幅频响应特性一般不如全向麦克风平坦。
对比此类方案，本发明的参考麦克风采用普通的全向型麦克风，结合阵列信号处理技术中的波束成形算法，使得阵列获取类似的空间指向性。
将阵列引入ANC系统的技术较早见于Daniel J.Sebald于1996发表的会议文章《Application of MVDR beamforming to reject turbulence noise in a duct》，此文采用阵列波束成形技术的主要目的在于降低风机管道内湍流噪声对参考信号的影响。
对比此类系统，本发明则充分考虑了风机管道系统这一应用场景中污染风机旋转噪声的各种主要相干及非相干干扰因素，同时考虑了参考信号的延时补偿问题以提高所构造的参考信号的质量。
发明内容
本发明要解决的技术问题是为风机管道前馈式ANC系统提供高质量的参考信号，本发明的参考信号构造算法能够运用阵列的空间指向性定向提取风机产生的旋转噪声，同时抑制次级扬声器的声反馈及其他方向的干扰信号，并充分保留旋转噪声的振幅与相位信息，提高ANC系统的稳定性与降噪性能。
本发明的技术方案如下：
一种基于麦克风阵列的风机管道系统ANC参考信号构造方法，在风机管道内布置好前馈式ANC系统，ANC控制器完成两路音频输入振幅校准以及次级通道离线辨识后，依次打开麦克风阵列设备和ANC控制器；具体流程如下：
步骤1：对阵列阵元位置矢量、对角加载量和FFT窗口等参数及用于存放算法输出序列等数据的容器进行初始化。
步骤2：设定一个周期的最大用时Tps，并启动计时器开始计时。
步骤3：若为首次启动，即处理第一帧数据，则求各通道基于FFT的旋转噪声频率估计值，并以各通道频率估计值的平均值作为下一帧数据的旋转噪声先验频率信息；若非处理第一帧数据，则接步骤4。在基于FFT的旋转噪声频率估计流程中，数据序列进行FFT变换之前需要先进行加窗操作，进而用二次插值法，结合FFT模值序列的最大值点或极大值点估计出风机旋转噪声(周期性噪声)的频率。首帧数据完成频率估计后，计时器计算得出此时算法流程的流逝时间telapse，流程休眠(Tps-telapse)后返回等待下次调用。由于首帧估计未进行实质的参考信号构造，该步骤输出的构造信号为步骤1所述容器的初始化值(一般初始化为全0序列)。
步骤4：下次调用(一帧数据采集完毕)，由步骤2开始进入流程，进行阵列数据预处理。此步骤以前一帧得到的频率估计值为中心频率设计巴特沃斯陷波器参数及结合阵列坐标系得出初级噪声源方向的阵列流形矢量等阵列参数。并对原始数据序列进行前后向滤波操作得出底噪声，用阵列原始数据减去所得底噪声获取风机噪声中的旋转噪声分量。接着对各通道旋转噪声分量进行均匀化以消除振幅差异，并对均匀化后的各通道数据求希尔伯特变换以使数据带有瞬时相位信息。
步骤5：对步骤4所得的阵列数据进行基于FBSS的MPDR或LCMP波束成形。若阵列配置为均匀环形阵列，则首先对阵列进行圆谐波变换得到虚拟线性阵列。接着对均匀线性阵列或虚拟线性阵列，进行子阵划分并执行前后向空间平滑(FBSS)操作，对平滑后的数据协方差矩阵进行对角加载以减小阵列失配的不良影响。进而执行MPDR或LCMP波束成形，取序列的实部作为波束成形器的输出。
步骤6：对步骤5波束成形器的输出进行基于FFT的旋转噪声频率估计。
步骤7：依据步骤6旋转噪声频率估计值、多频点信号表达式和步骤5波束成形器输出构造矩阵方程。用最小二乘法得出该矩阵方程的解向量，进而解出该解向量封装的信号的振幅与相位信息。至此流程已经获取了旋转噪声的频率、振幅、相位和一帧数据固有时长四个信息，依据这些信息，配合计时器与算法时长设定值，可构造出ANC控制器所需的风机旋转噪声参考信号。往后当前帧所得频率估计值和振幅估计值都将作为信号的先验信息为下一帧所用。该步骤描述的振幅估计值和步骤3描述的基于FFT的旋转噪声频率估计值均可添加滑动平均操作，使参数更加稳定。
步骤8：计时器计算得出此时算法流程的流逝时间telapse，流程休眠(Tps-telapse)后输出步骤7所得风机旋转噪声参考信号并返回，等待下次调用。下次调用由步骤2开始进入流程。
本发明的有益效果：对于风机管道ANC系统而言，本发明的参考信号构造算法能使得旋转噪声避免被次级扬声器的声反馈等因素污染，且保留了参考信号的振幅与相位信息，极大地提高了ANC系统的稳定性及降噪效果。参考麦克风阵列子系统可用低成本的Linux设备实现，独立于ANC控制器，由于所构造的参考信号质量较高，可使ANC控制器从参考信号预处理的运算压力中解放出来，执行FxRLS等更复杂的前馈式ANC算法以及进行更精确的次级通道建模。
附图说明
图1是整体系统结构框图。
图2(a)是阵列坐标系示意图。
图2(b)是ULA阵元布置示意图。
图2(c)是UCA阵元布置示意图。
图3是本发明参考信号构造算法的主流程图。
图4是本发明基于FFT的旋转噪声频率估计流程图。
图5(a)是本发明阵列数据预处理流程图；
图5(b)是本发明基于FBSS的MPDR或LCMP波束成形流程图；
图5(c)是本发明基于LS的振幅、相位估计与信号构造流程图；
图6是采样自真实风机的单频旋转噪声示例图；(a)是时程图；(b)是时程局部放大图；(c)是功率谱图；
图7是采样自真实风机的含谐波旋转噪声示例图；(a)是时程图；(b)是时程局部放大图；(c)是功率谱图；
图8是一个小尺寸阵列原始数据通道振幅差异现象示例图；(a)是时程差异对比图；(b)是功率谱差异对比图。
具体实施方式
以下结合发明内容和说明书附图详细说明本发明的具体实施方式。
图1为整体系统结构框图，指示了系统在风机管道内的位置示意与工作原理。图中参考麦克风阵列设备、ANC控制器(包括三个音频接口和接口以下部分)、次级扬声器和误差麦克风共同构成前馈式ANC系统。其中参考麦克风阵列为Linux设备，ANC控制器为DSP或FPGA设备。参考麦克风阵列阵元由全向麦克风构成，误差麦克风为一只全向型驻极体麦克风。在FxLMS等自适应ANC算法的作用下，ANC控制器通过次级扬声器产生与初级噪声(来自于管道风机)振幅相同但相位相反的次级噪声，二者相干叠加，达到风机管道噪声抑制的目的。ANC控制器需要预先完成两路音频输入的振幅校准(即同一信号相同声压激励下，两音频输入通道需得到相同采样值)和次级通道的离线辨识。
进行阵列布置时，可将参考麦克风阵列布置为均匀线性阵列(Uniform LinearArray,ULA)或均匀环形阵列(Uniform Circular Array,UCA)，阵列不可过于靠近初级噪声源以保证声波到达阵列时，波阵面为平面波模型。同时在构造算法具体实现之前，首先为算法明确阵列坐标系。如图2所示，无论是ULA还是UCA，通常以阵列的中心(圆心)作为阵列坐标系原点，且将第0号阵元布置在三维直角坐标系x轴正向的第一个位置上。参考麦克风阵列实物安装时，将x轴正向对准初级噪声源方向(即图1所示管道风机方向)，x轴负向对准次级扬声器方向，即在ULA情形下，阵元均匀分布在x轴上，初级噪声源方向信号与声反馈干扰信号形成“端射”。对于窄带平面波，选定阵列参考系后，阵元接收到的频域数据表现为如下形式：
其中F(ω)为平面波信号波形的频域形式，
称为阵列流形矢量。
为阵元在三维直角坐标系中的位置矢量，
为波动方程定义的波数矢量，
为标准球坐标系相对于三维直角坐标系的方向余弦，c为声速。
在窄带波束成形技术中，波束成形器往往表现为一系列复权值
由线性系统时域卷积对应频域乘积可知，波束成形器输出的频域表现形式为：
从(2)式可见阵列流形矢量
封装了阵列的空-时-频信息，此参量将直接运用于复权值
的计算。同时，将ULA情形带入，则可以发现这一参量具有Vandermonde形式，这是执行前后向空间平滑去相干的先决条件，UCA则可以通过变换得到虚拟线性阵列。
图3为本发明参考信号构造算法的主流程图。算法首先对阵列阵元位置矢量、对角加载量和FFT窗口等参数及用于存放算法输出序列等数据的容器进行初始化。在处理数据帧时，算法首先启用计时器开始计时，同时设定一个算法流程的最大用时Tps。此设定值应当小于一帧音频数据的最大采样延时，设一帧音频数据的长度为Nsample，采样频率为fs，则理想情况下一帧数据第一点延时为0，最后一点对应最大延时
例如，以8KHz为采样率进行采样，设定一帧音频数据为4096点，则对应的Tmaxdelay＝512ms，此时Tps可取400ms。在完成一帧数据处理后，再次启用计时器并计算得出此刻算法的耗时telapse，算法流程休眠Tps-telapse后返回，等待下次调用。其中一次算法流程得到的频率估计值和振幅估计值将作为下一次算法流程的旋转噪声信号先验信息。首次调用参考信号构造算法时(首帧估计时)，算法输出所得构造信号为相应容器的初始化序列，有效的输出将在后续的调用中得到。
由图3可见，本发明参考信号构造算法主要由“基于FFT的旋转噪声频率估计”、“阵列数据预处理”、“基于FBSS的MPDR或LCMP波束成形”和“基于LS的振幅、相位估计与信号构造”4个部分组成，各部分的流程图见图4和图5。
具体如下：
1.基于FFT的旋转噪声频率估计，流程图如4所示，其目的是获取图6(a)-(b)所示旋转噪声谱峰对应的频率，旋转噪声为单频噪声或多频点噪声时均适用：
为避免频谱泄漏及栅栏效应带来的不良影响，本发明基于二次函数顶点公式(二次插值)进行高效估计。此方法需要FFT振幅谱最值点或选取的极值点及其前后各一个数据点完成方程的构造，即不在一条直线上的三点确定一条抛物线。其中最大值点子方法确定的是旋转噪声的基频，其整数次谐波可直接乘以谐波序号得出估计；极大值点子方法可由差分法实现，可估计旋转噪声的基频和直接估计谐波频率，在谐波频率因为管道声学特性发生偏移时仍然适用，可避免最大值法对谐波估计的误差，但极大值点的选取必须满足高度阈值条件(所选极值点的幅度应大于一定高度)和距离阈值条件(两个极大值点的间隔应当大于一定频率)。由于实序列的FFT序列具有共轭对称的特性，故寻找最大值或极大值时，仅用振幅谱的正频率部分即可。在FFT变换之前，应当对原始数据序列进行加窗以提高频率估计精度，设寻得某个最大值点或极大值点数据组分别是(xf,vf)、(xc,vc)和(xb,vb)，其中xc是最大值点或极大值点(索引根据频率分辨率换算为频率)，vc表示该点幅值，则可构造线性方程组：
解出二次抛物线参数a、b及c后可根据抛物线顶点的横坐标公式求得频率估计值，为避免此范德蒙矩阵求逆，也可直接解析抛物线顶点横坐标，具体公式为：
当所加的窗为高斯窗时，式(4)中vf、vc和vb应取自然对数。对于一次算法流程的估计值，无论是频率还是后续提及的振幅，均可以用滑动平均窗口进行平滑，使得在稳态环境中(认为旋转噪声的频率是不变的或仅仅缓慢地在一个小范围内变化)，估计参数的变动减小，算法输出数据帧与帧之间的衔接更加光滑。
2.阵列数据预处理，流程图如图5相应部分所示，主要完成风机噪声中旋转噪声分量的提取、阵列数据通道均匀化和希尔伯特变换：
步骤1：以前一帧数据估计出的频率fest为中心频率设计巴特沃斯陷波器参数、声波波数及波长意义下的阵元位置矢量等MPDR/LCMP波束成形算法需要使用的参量。设阵列接收的原始数据矩阵为
利用巴特沃斯滤波器对Rx的各行(相当于阵列的各通道)进行滤波得到底噪声数据矩阵Rn，利用Rs＝Rx-Rn得出旋转噪声数据矩阵。巴特沃斯陷波器的阶数不可过高，同时阻带宽度不可过窄，以一帧数据有4096点为例，滤波器的阶数可取6阶，截止频率可设置为fest±2.5Hz。滤波应当采用前后向滤波的形式进行，目的在于依据线性系统理论抵消数字滤波器带来的信号延时，同时为避免滤波后数据首尾畸变严重，可保留部分的前一帧原始数据(如前一帧原始数据的后2048个点)共同参与滤波。
步骤2：对Rs各通道进行均匀化，即只考虑不同阵元间的相位差异，不考虑振幅差异。图7展示了小尺寸阵列原始数据通道振幅差异现象的示例，这相当于为阵元加入了一个不确定的滤波窗口，影响波束成形器的性能。本设计采用将Rs的FFT序列各点模值归一到平均值来实现均匀化，这样做可充分保留旋转噪声的幅值信息。
步骤3：将Rs各通道进行希尔伯特变换，得到信号的解析形式数据矩阵Ra，使得数据序列带有瞬时相位信息，使波束成形器得到的复权值矢量可对其进行增益-移相操作。
3.基于FBSS的MPDR或LCMP波束成形，流程图如图5相应部分所示，目的是使阵列仅提取初级噪声源(即管道风机)方向上的旋转噪声，而抑制其他方向上(如次级扬声器方向)的干扰：
步骤1：若阵列设备配置为均匀线性阵列，则进行基于前后向空间平滑(Forward-Backward Spatial Smoothing,FBSS)的MPDR或LCMP波束成形；配置为均匀环形阵列，则须先进行圆谐波变换化为虚拟线阵，使得其阵列流形矢量具备范德蒙形式。按照系统确定的阵列坐标系，定义圆谐波变换矩阵T为：
其中n＝-N,-N+1,…,0,…N-1,N为圆谐波序号，N为选取的最大圆谐波阶数；m＝0,1,…,M-1,M为阵元数量，
为UCA第m号阵元在阵列坐标系中的方位角。定义圆模态响应消去阵J为：
其中Jn(kr)为第一类贝塞尔函数，k为波数矢量
的模值，r为均匀圆环阵列的半径，则对于UCA，阵元数据域向圆谐波数据域转换的过程为：
步骤2：为避免相干干扰对波束成形器性能带来的不良影响，需要使用前后向空间平滑处理。首先需要对均匀线性阵列或经过上述圆谐波变换形成的虚拟线性阵列进行子阵划分，以阵元数M＝5的均匀线性阵列为例，若其阵元依次表示为ARRAY＝{e0,e1,e2,e3,e4}，则可划分为3个具有三个阵元的子阵为ARRAYsub0＝{e0,e1,e2}、ARRAYsub1＝{e1,e2,e3}及ARRAYsub2＝{e2,e3,e4}，分别对应经过前述步骤处理得到的解析形式数据子矩阵Ra0、Ra1及Ra2。其中第i号子阵的采样数据协方差矩阵为：
空间平滑处理为：
其中K+1是子阵数量，引入前后向平均得到最终的前后向空间平滑采样数据协方差矩阵为：
其中Js为与子阵阵元数同阶的副对角线为1，其余位置为0的交换矩阵，Cx,ss *表示对Cx,ss元素取共轭，JsCx,ss *Js又称为后向采样数据协方差矩阵。
步骤3：对Cx,fbss进行对角加载(Diagonal Loading,DL)以降低阵元位置失配或波达方向失配等问题发生时，对波束成形器性能带来的消极影响。其本质是对采样数据协方差矩阵人为地注入一功率相对较低的白噪声，即加载方法为：
其中
称为对角加载量，I为与Cx,fbss同阶的单位矩阵。最优对角加载量在不同情形下无统一的解析表达形式，通常可依据原始信号的信噪比或经过前述步骤处理得到的数据序列功率进行选择，灵活性较高。
步骤4：进行波束成形器权值矢量设计，其中MPDR最优权值矢量为：
其中
表示期望方向的阵列流形矢量(若依据步骤1配置阵列坐标系，则其对应x轴正向，即初级噪声源的方向)，L表示子阵阵元数量。
LCMP最优权值矢量为：
其中
称为约束值矢量，而F称为约束矩阵。约束设置较为灵活，常见的有无畸变约束和零点约束，在本发明的应用场景中，约束可设置为：
其中
分别为关于阵列坐标系原点对称的子阵的期望方向和零陷方向的阵列流形矢量。如前所述，在本发明的应用场景中，
对应初级噪声源方向，
对应次级扬声器方向。
本发明推荐使用LCMP最优权值矢量进行波束成形，经过上述步骤的处理，设关于阵列坐标系原点对称的子阵的数据矩阵为Ra,sym，则某频率点n对应的最终的波束成形器输出为：
4.基于LS的振幅、相位估计与信号构造，流程图如图5相应部分所示，目的是提取旋转噪声的振幅和相位，并结合一帧数据的音频延时及设定的算法流程的最大用时Tps完成延时补偿，使得参考麦克风阵列输出给ANC控制器的信号由已经流逝的信号转变为将要发生的信号：
具有N次谐波的周期性信号或N+1点频率的多频信号可表示为：
由(16)式可见，各分量的频率、振幅、相位和作用时间后，即可构造最终满足上述需求的送入ANC控制器参与运算的参考信号。作用时间构造为时刻向量：
基于最小二乘法，由(16)式构造参数矩阵：
解矩阵方程
此处
为式(15)所有频率点合成的最终信号，其中
则显然振幅an为un和vn的2-范数，根据
及
的反函数符号及数值可解出无歧义相位
至此算法流程已获取旋转噪声的频率、振幅、相位和作用时间四个构造要素。其中频率和振幅可加入滑动平均操作，使得参数估计更为稳定。
最终的参考信号构造为：
其中“+Tmaxdelay”一项表示对正弦序列进行左移以补偿音频采样带来的延时，“+Tps”一项表示对序列进行左移以补偿参考信号构造算法的耗时。
Claims (4)
1.一种基于麦克风阵列的风机管道ANC系统参考信号构造方法，其特征在于，步骤如下：
步骤1：初始化参数和数据容器；
步骤2：设定一个算法周期最大用时Tps，并启动计时器开始计时；
步骤3：首次启动，处理第一帧数据，求各通道基于FFT的旋转噪声频率估计值，并以各通道噪声频率估计值的平均值作为下一帧数据的旋转噪声先验频率信息；
非首次启动时，直接进入步骤4；
步骤4：预处理阵列数据；以前一帧得到的频率估计值为中心频率设计巴特沃斯陷波器参数，结合阵列坐标系得出初级噪声源方向的阵列流形矢量；对原始数据序列进行前后向滤波操作得出底噪声，用原始数据序列减去所得底噪声，获取风机噪声中的旋转噪声分量；接着对各通道旋转噪声分量进行均匀化以消除振幅差异，并对均匀化后的各通道数据求希尔伯特变换以使数据带有瞬时相位信息；
步骤5：对步骤4预处理后的阵列数据进行基于前后向空间平滑操作FBSS的MPDR或LCMP波束成形；阵列设备配置分为均匀环形阵列和均匀线性阵列；阵列设备配置为均匀环形阵列时，首先对阵列进行圆谐波变换得到虚拟线性阵列；接着对配置为均匀线性阵列或虚拟线性阵列的阵列设备进行子阵划分并执行前后向空间平滑操作，对平滑后的数据协方差矩阵进行对角加载；进而执行MPDR或LCMP波束成形，取序列的实部作为波束成形器的输出；
其中
为波束成形权值矢量，Ra,sym为关于阵列坐标系原点对称的子阵的数据矩阵，real(·)表示取实部操作，(·)T表示转置操作，Nsample表示一帧数据的点数；
步骤6：对步骤5波束成形器的输出进行基于FFT的旋转噪声频率估计；
步骤7：依据步骤6旋转噪声频率估计值、多频点信号表达式和步骤5波束成形器，输出构造矩阵方程；
多频点信号表达式：
其中an与
分别是频率为fn的信号分量所对应的振幅与相位，风机旋转噪声频率fn由步骤6估计得出，式中
作用时间构造为时刻向量：
基于最小二乘法，构建参数矩阵方程：
解矩阵方程
此处
为式(1)所有频率点合成的最终信号，其中
振幅an为un和vn的2-范数，根据
及
的反函数符号及数值解出无歧义相位
基于旋转噪声的频率、振幅、相位和时刻向量四个信息，配合计时器与算法时长设定值，构造风机旋转噪声参考信号；
其中“+Tmaxdelay”表示对正弦序列进行左移以补偿音频采样带来的延时，“+Tps”一项表示对序列进行左移以补偿参考信号构造算法的耗时；
步骤8：计时器得出此时算法流程的流逝时间telapse，流程休眠(Tps-telapse)后输出步骤3各通道噪声频率估计值的平均值作为步骤1所述数据容器的初始化值等待下次调用；输出步骤7所得风机旋转噪声参考信号并返回，等待下次调用，下次调用由步骤2开始进入。
2.根据权利要求1所述的基于麦克风阵列的风机管道ANC系统参考信号构造方法，其特征在于，所述步骤5中基于FBSS的MPDR或LCMP波束成形，具体步骤如下：
步骤5.1：阵列设备配置为均匀线性阵列，进行基于前后向空间平滑FBSS的MPDR或LCMP波束成形；配置为均匀环形阵列，先进行圆谐波变换化为虚拟线性线阵；定义圆谐波变换矩阵T为：
其中n＝-N,-N+1,…,0,…N-1,N为圆谐波序号，N为选取的最大圆谐波阶数；m＝0,1,…,M-1,M为阵元数量；
为均匀环形阵列第m号阵元在阵列坐标系中的方位角；定义圆模态响应消去阵J为：
其中Jn(kr)为第一类贝塞尔函数，k为波数矢量
的模值，r为均匀环形阵列的半径，j为虚数单位；
阵元数据域向圆谐波数据域转换的过程为：
其中(·)H表示共轭转置操作，
为阵元数据域下的阵列流形矢量，
为经过阵列数据预处理，各通道数据序列以解析形式存在的数据矩阵，
及RCH分别为均匀环形阵列圆谐波域下的阵列流形矢量和数据矩阵；
步骤5.2：前后向空间平滑处理；首先对均匀线性阵列或经过上述圆谐波变换形成的虚拟线性阵列进行子阵划分，其中第i号子阵的采样数据协方差矩阵为：
其中Rai为第i号子阵的数据矩阵；
空间平滑处理为：
其中K+1是子阵数量，引入前后向平均得到最终的前后向空间平滑采样数据协方差矩阵为：
其中Js为与子阵阵元数同阶的交换矩阵，其副对角线为1，其余位置为0；Cx,ss *表示对Cx,ss元素取共轭，JsCx,ss *Js又称为后向采样数据协方差矩阵；
步骤5.3：对Cx,fbss进行对角加载；加载方法为：
其中
称为对角加载量，I为与Cx,fbss同阶的单位矩阵；
步骤5.4：进行波束成形器权值矢量设计，其中MPDR最优权值矢量为：
其中
表示期望方向的阵列流形矢量，L表示子阵阵元数量；
LCMP最优权值矢量为：
其中
称为约束值矢量，F称为约束矩阵。
3.根据权利要求1或2所述的基于麦克风阵列的风机管道ANC系统参考信号构造方法，其特征在于，所述步骤6中基于FFT的旋转噪声频率估计值流程，原始数据序列加窗进行FFT变换；取FFT序列前半段的模值并将序列索引依据FFT频率分辨率换算为实际频率；取模值序列的最大值点或同时满足高度阈值条件与距离阈值条件的极大值点及其前后各一点用于二次插值；当FFT变换所加窗为高斯窗时，先对所取点的模值取自然对数；再利用二次插值法，结合所取点计算出风机旋转噪声频率估计值；风机旋转噪声频率估计值为二次插值求出的频率估计值直接输出或经滑动平均化后输出；首帧数据完成频率估计后，计时器得出此时算法流程的流逝时间telapse，流程休眠(Tps-telapse)后返回等待下次调用。
4.根据权利要求1所述的基于麦克风阵列的风机管道ANC系统参考信号构造方法，其特征在于，所述最大用时Tps设定值小于一帧音频数据的最大采样延时。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111610076.2A CN114333753A (zh) | 2021-12-27 | 2021-12-27 | 一种基于麦克风阵列的风机管道anc系统参考信号构造方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111610076.2A CN114333753A (zh) | 2021-12-27 | 2021-12-27 | 一种基于麦克风阵列的风机管道anc系统参考信号构造方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114333753A true CN114333753A (zh) | 2022-04-12 |
Family
ID=81012629
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111610076.2A Pending CN114333753A (zh) | 2021-12-27 | 2021-12-27 | 一种基于麦克风阵列的风机管道anc系统参考信号构造方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114333753A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114842821A (zh) * | 2022-05-06 | 2022-08-02 | 齐鲁工业大学 | 基于改进FxLMS算法的噪声主动控制方法及系统 |
| CN116624793A (zh) * | 2023-07-25 | 2023-08-22 | 上海电机学院 | 一种双指向性超低压气体管道泄漏声波信号检测方法 |
-
2021
- 2021-12-27 CN CN202111610076.2A patent/CN114333753A/zh active Pending
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114842821A (zh) * | 2022-05-06 | 2022-08-02 | 齐鲁工业大学 | 基于改进FxLMS算法的噪声主动控制方法及系统 |
| CN116624793A (zh) * | 2023-07-25 | 2023-08-22 | 上海电机学院 | 一种双指向性超低压气体管道泄漏声波信号检测方法 |
| CN116624793B (zh) * | 2023-07-25 | 2023-10-31 | 上海电机学院 | 一种双指向性超低压气体管道泄漏声波信号检测方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105590631B (zh) | 信号处理的方法及装置 | |
| CN108172235B (zh) | 基于维纳后置滤波的ls波束形成混响抑制方法 | |
| Benesty et al. | Conventional beamforming techniques | |
| Fisher et al. | Near-field spherical microphone array processing with radial filtering | |
| Yan et al. | Optimal modal beamforming for spherical microphone arrays | |
| CN108600894B (zh) | 一种耳机自适应有源噪声控制系统及方法 | |
| CN114333753A (zh) | 一种基于麦克风阵列的风机管道anc系统参考信号构造方法 | |
| US8712075B2 (en) | Spatially pre-processed target-to-jammer ratio weighted filter and method thereof | |
| KR100856246B1 (ko) | 실제 잡음 환경의 특성을 반영한 빔포밍 장치 및 방법 | |
| KR102573843B1 (ko) | 음성 제어를 갖는 낮은 복잡도의 다중 채널 스마트 라우드스피커 | |
| CN102440002A (zh) | 用于传感器阵列的优化模态波束成型器 | |
| JP2001510001A (ja) | 複数ソースを伴うオーディオ処理装置 | |
| Saito et al. | Influence of modeling error on noise reduction performance of active noise control systems using filtered-x LMS algorithm | |
| CN105981404A (zh) | 使用麦克风阵列的混响声的提取 | |
| JPWO2017002525A1 (ja) | 信号処理装置、信号処理方法、および信号処理プログラム | |
| Daniel et al. | Time domain velocity vector for retracing the multipath propagation | |
| CN112331226B (zh) | 一种针对主动降噪系统的语音增强系统及方法 | |
| JPH10207490A (ja) | 信号処理装置 | |
| CN110111802B (zh) | 基于卡尔曼滤波的自适应去混响方法 | |
| CN113889136A (zh) | 一种基于麦克风阵列的拾音方法、拾音装置及存储介质 | |
| CN111312269A (zh) | 一种智能音箱中的快速回声消除方法 | |
| Thomas et al. | Eigenvalue equalization filtered-x algorithm for the multichannel active noise control of stationary and nonstationary signals | |
| CN111263291B (zh) | 一种基于高阶麦克风阵列的声场重构方法 | |
| Wang et al. | Kronecker product adaptive beamforming for microphone arrays | |
| CN113948101A (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 |