采样示波器（精选5篇）
采样示波器 篇1
1 系统硬件设计
数字示波器总体结构如图1所示。
1.1 FPGA硬件电路设计
FPGA完成时基的产生和采集时序控制功能。由于所设计的示波器输入频率范围较宽,系统采用等效和实时两种采样方式。通过FPGA对已处理过的输入信号进行测频,若输入频率小于1 MHz,选用实时采样;反之选用等效采样。
1.1.1 A/D驱动模块
LTC1196是一款串行A/D,在时钟的控制下,每12个节拍即可进行一次A/D转换。当CS跳变为低电平时,即开始转换过程,转换结果由dout脚串行输出。可通过编写相应的模块来控制LTC1196的工作。此模块采用状态机实现。模块通过AdEn启动转换,当rdy为“1”时转换完成,选中dout读取采样值。AD空闲时idle为“1”。具体见图2。
1.1.2 采样时序控制模块
采样时序控制模块是本设计的核心,由锁相环产生120MHz的时钟驱动,是等效采样中高精度延时的发生器,并能根据触发信号和AD的状态精确地控制信号的保持和AD转换的启动。具体的工作过程如图3所示。
系统进入等效采样模式后,当第一次触发信号到来时,模块内部对120M时钟计数“1”,计满后将当前的信号保持,然后启动AD进行数模转换;当第二次触发信号到来时对时钟计数“2”,然后保持,启动AD,以此类推。这样就能用等效采样的原理采集一个波形,当AD转换的数据写满FIFO后由ARM给出复位信号,模块重新从“1”开始计数。具体框图见图4。
1.2 信号调理电路设计
信号的调理是将输入的信号经放大、量程转换、保护、滤波、线性化等操作,使其变成一个能够被A/D正常转换的信号。本设计中,信号调理电路分电压跟随、电压放大、加法电路及比较选择电路四部分,电路图分别如图5、图6所示。示波器的垂直分辨率分为三档:1V/div,0.1V/div,10mV/div。放大倍数分别设计为1.1倍、10倍和100倍,根据输入信号来选用不同的放大倍数。由于放大电路采用了正负双电源,所以应在A/D前将放大后的电压抬高2.5V,使得输出时信号在A/D转换所要求的0~5V范围之内。此外,在输入端加入电压跟随,可以提高输入电阻,提高系统性能。
1.3 采样保持电路设计
采样保持器跟踪输出来自输入端的模拟信号,在给出保持命令的时刻,将输入端的信号无延时、准确、长时间地保持在其输出端,直到保持命令撤消为止。由于常见的集成采样保持器(如LF398)的性能不符合本系统的要求,所以需要用分立元件按照采样保持器结构图搭建一个采样保持电路,如图7所示。
1.4 A/D转换电路设计
LTC1196是一款采样率可达1MHz的8位串行A/D,能工作在3V或5V的电压下,典型环境下消耗功率仅为50mW。LTC1196通过A/D通过三个I/O口与FPGA相连,分别输出片选、时钟、数据信号,由FPGA控制其转换的起停与数据的读取。其电路如图8所示。
2 系统软件设计
数字示波器采用以C8051F020为控制核心,完成人机交互、数据提取、数据处理和波形显示等功能,相关软件采用C语言在Keil集成开发环境中编写、调试和下载。
主程序负责数据的读取、处理、更新,以及人机交互方面的处理,其流程如图9所示。
定时中断负责波形的刷新。波形的显示是用D/A在模拟示波器的显示屏上画点实现的,为使画出的点均匀就需要用定时器来控制D/A的时间,定时器中断流程如图10所示。图中Da_flag为D/A显示完成标志符,n、m是用于计取数据的序号变量。在D/A显示好一屏数据之后,将D/A值置0,在屏幕的最下方画直线直至新一屏数据准备完毕。
3 测试结果及分析
3.1 系统带宽的测量
测量方法:用信号发生器产生特定频率的正弦信号输入到被测示波器中,用另外的双踪示波器分别对被测示波器输入端与输出端的波形进行测量,重点读取其输入电压峰峰值Vpp和输出电压峰峰值V'pp。
测量数据见表1。
通过表1可以看出,选择不同的垂直分辨率档位系统的带宽是不同的。其中垂直分辨率为1V/div时,带宽约为2.3MHz、垂直分辨率为0.1V/div时,带宽约为1.6MHz,垂直分辨率为10mV/div时,带宽约为1.3MHz。系统具备测量10Hz~1MHz波形的能力。
3.2 垂直灵敏度的测量
测量方法:向被测示波器输入一个正弦信号,从被测示波器读取该正弦信号峰峰值,同时用模拟示波器测取该值作为标准电压。通过比较两个数值得到测量误差。
测量数据见表2。
通过表2可以看出,1V/div档和0.1V/div档的相对误差均小于5%;0.1V/div档时,相对误差较大。
3.3 扫描速度的测量
示波器设计的不同扫描速度如表3所示。在被测信号为500Hz、5kHz、50kHz和500kHz时,波形刚好为各档位的一个水平扫描周期,周期测量误差可以用标准信号进行调整,完全满足波形周期测量误差≤5%的要求。
4 结束语
同时采用等效采样和实时采样两种方式,等效采样率可达100MSa/s,有效地扩展了示波器的工作带宽,能复现频率大大超过奈奎斯特极限频率的信号波形。
参考文献
[1]朱英时.数字示波器原理与测试分析[J].[出版地不详]:测试技术,2007(Z1).
[2]贾春霞,张洪艳.数字存储示波器现状初探[J].天津:仪器仪表用户,2005(1).
[3]蒋焕文.电子测量[M].北京:中国计量出版社,1988,5:29-35.
[4]陈尚松,等.电子测量与仪器[M].北京:电子工业出版社,20-05-01:56-67.
[5]顾德英,等.计算机控制技术(第2版).北京:北京邮电大学出版社,2007-04:48-50,156-161.
[6]韩峰.单片机等效采样示波器的设计[J].鸡西:鸡西大学学报,2006(12).
[7]刘国林.电子测量[M].北京:机械工业出版社,2003:103-152.
[8]顾德英,等.计算机控制技术(第2版).北京:北京邮电大学出版社,2007-04:48-50,156-161.
采样示波器 篇2
关键词：非理想带通采样,正交变换,低通滤波法,优化设计,FPGA
在信号处理领域,正交相干检波保留了信号复包络的所有信息,并在雷达、通信等电子系统中得到了广泛的应用。图1所示为传统I/Q正交化实现方法,传统I/Q正交化通过模拟移相器移相90°实现,但模拟器件性能的限制使得I/Q两路并不能完全正交并存在I/Q幅度的不一致性。图2中数字中频接收机可以较好的解决这种不一致性。随着数字信号处理技术的发展及高速A/D、DSP芯片的不断出现,使得直接中频采样成为可能,即直接将模拟的中频信号用A/D变成数字信号,在数字域将其分成I和Q信号。与模拟方法相比,直接中频采样具有非常高的精度和稳定性。
对于某些特定信号或由于工艺原因,中频信号会产生频偏,给采样及正交检波带来不利影响,文中主要针对中频信号的频偏对采样及正交检波的影响进行了分析,并提出解决方法。
数字中频采样主要有以下几种方法:低通滤波法、数字乘积检波法、希尔伯特变换以及贝赛尔插值法。这几种方法本质上相同,其中以低通滤波法在工程中最为常用。
1 带通信号采样定理
1.1 中频采样定理的局限性
信号进行采样的必要条件是采样后的频谱不发生混频。设频率带限信号x(t),频率限制在(fL,fH)内,采样率满足fs≥2B,且
式中M取能满足fs≥2(fH一fL)的最大整数。采用式(1)可以简化数字正交检波器的设计。但在某些特定情况下,产生的信号会有频偏,如高斯最小键控频移信号(GMSK),由于其信号的固有特性,使其产生的频偏不能被消除,这时载频加上频偏就得到一个非整数的中频信号f0,如仍按上式选取采样频率会存在一定问题,如当f0=30.333MHz,B=6 MHz,当M=1,2,3,4,5时代入式(1)得到的fs值。
式(1)中并没有考虑混频器的参考频率,而实际上,为了避免组合频率的干扰,系统对接收机中各混频器有严格要求,同时,非相参振荡源对系统也有不良影响[1]。理论上说取采样频率可以通过本振信号分频或者倍频来得到。但假设本振信号频率为100 MHz,要通过分频产生表1中的17.333 MHz采样频率fs是很困难。在工程中产生这样频率的采样信号是很不方便,因此采用式(1)来得到采样频率具有局限性[2]。
1.2 带通信号采样定理
带通信号采样定理在工程中应用身份广泛。设带通信号中心频率为f0,带宽为B,则其上下截止频率分别为fH=f0+B/s和fL=f0-B/2,由采样值不失真地重建信号的充要条件可得采样率fs满足
其中m=1,…,mmax,mmax=[fH/B],[x]为不大于x的最大整数。
式(2)表明带通信号的采样频率取值范围由mmax个互不重合的区间Sm=[2fH/m,2fL/(m-1)]组成,即。S1,S2…,Smax对应不失真采样频率范围,最低不失真采样频率为fsmin=2fH/mmax。S1=[2fH,+∞)对应低通信号采样定理所确定的采样频率范围,若将低通信号看作频谱分布下界为零的带通信号,则带通信号采样定理包含了低通采样定理[3]。
文中所讨论雷达接收机接信号中频f0为30.416 667 MHz,带宽B为6.667 MHz,则接收到的带通信号上截止频率fH为33.750 167 MHz,下截止频率fL为27.083 167 MHz。按照上述定理可得mmax=5,则m∈{1,2,3,4,5},表2为对应不同m值的取值范围。
带通信号采样定理确定了采样频率的取值范围,在实际工程中,非整数的采样率较难实现,因此可取fs为40 MHz。但采用这个频率有不利影响:
(1) I/Q正交中的系数不再是{1,0,-1,0,…}和{0,-1,0,1,…}序列。本文中载频和采样频率比值等于fs/f0=40/30.416 667,虽然采样后正弦或余弦函数仍为周期序列,但此时周期为N=4×107,给实际工程中应用带来了不便;
(2)对低通滤波器的设计要求高,信号负频谱不处于两相邻的正频谱中间,导致两边的过渡带宽度不等,需要选取窄的过渡带宽作为滤波器的过渡带,这就要求对滤波器进行优化设计。图3为信号中频采样后左右过渡带宽不等的示意图,其中虚线是设计低通滤波器时需要满足的频域特性。
2 中频采样滤波器的优化设计
2.1 滤波器参数的确定
输入信号基带和倍频分量的频率范围,基带与倍频间距以及带通采样定理是设计低通滤波器的依据。过渡带越宽,滤波器越容易实现。由文献[4]计算可以得到带通信号的参数。
基带频率范围对应为滤波器的通带为:[0,0.166 7π];倍频分量频率范围对应为滤波器的阻带为:[0.791 7π,π];基带与倍频间距对应为滤波器的过渡带宽为:min{0.625π,0.708 4π}。基带与其两边倍频间距不一样,过渡带宽只能取基带与倍频间距的最小值才能满足设计要求。
2.2 滤波器设计方法
滤波器参数确定以后,就可根据实际需要选用滤波器设计方法。滤波器分为FIR和IIR两种。IIR数字滤波器利用模拟滤波器成熟的理论及设计图进行设计,保留了一些典型模拟滤波器优良的幅度特性,但为了得到线性相位,需要另外增加相位校正网络,使滤波器设计变得复杂,成本也高。与IIR滤波器相比,FIR的实现是非递归的,稳定性好,精度高;更重要的是FIR滤波器在保证幅度特性满足要求的同时,很容易做到严格的线性相位[5,6],因此,FIR滤波器在高保真的信号处理中得到广泛应用,文中采用FIR滤波器进行设计。
FIR滤波器的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题,设计方法主要有窗函数法,频率采样法和等波纹逼近法等。频域采样法适合于设计具有任意幅度特性的滤波器,但边界频率不易控制,要获得较好的边界频率只有通过增加采样点数,但采样点数的增加则必须使用更多的硬件资源,不利于工程实现,因此这里不宜采用频域采样法。同时考虑到所采用的FPGA器件的乘法器资源,希望在满足性能指标的前提下滤波器的阶数尽量少,故本文中滤波器系数取8阶,下面主要对相同阶数下3种不同滤波器的性能进行分析。
2.3 滤波器设计方法及性能比较
文中采用Matlab的FDATool(Filter Design&A-nalysis Tool)来设计滤波器,FDATool是Matlab信号处理工具箱中专门用于滤波器设计与分析的工具,他可以设计几乎所有的常规滤波器,并且使用十分简便。这里主要针对于工程,因此,滤波器阶数不能很大,最好为确定值,将其定为8阶,而在FDATool中滤波器阶数并不能手动调节,需要多次试探后才能得到所需要的滤波器阶数,可以将设计好后的滤波器进行性能分析,看是否满足要求,如不满足则要求重新设计。
以同一指标(N=8,通带频率wp=0.166 7π,采样频率fs=40 MHz,截止频率wc=0.479 2π),分别用窗函数法中的Kaiser窗和Hamming窗以及等波纹逼近法中的切比雪夫逼近法设计低通滤波器,将生成的滤波器系数导出到Matlab中并进行性能分析。
图4是3种低通滤波器幅值特性结果对比,从图中可以看出用Kaiser窗设计的主瓣最宽,旁瓣起伏最小;用Hamming窗设计的滤波器主瓣宽度比Kaiser窗的窄,但旁瓣起伏要大,用切比雪夫逼近法所设计滤波器主瓣宽度在Hamming窗与Kaiser窗之间,旁瓣起伏最小,且最小阻带衰减最高,而实际应用中的滤波器要求在保证主瓣宽度达到一定要求的前提下,适当牺牲主瓣宽度来换取旁瓣波动的减少,因此从幅频响应图中可以看出切比雪夫逼近法设计的滤波器性能较窗函数法更加优秀。图5为分别采用这3种方法设计的滤波器作为数字正交检波器的低通滤波器,各滤波器的镜频抑制比(IR),从图中可以看出,切比雪夫逼近法抑制镜频效果最好,平均镜频抑制比达到了-75.542 dB。
从对图4以及图5的分析中可以看出,采用相同的阶数,切比雪夫逼近法得到的滤波器具有更好的性能。
文中采用fs=40 MHz造成采样后信号左右过渡带宽不等,因此,滤波器设计时需考虑过渡带宽的选择问题。针对文中过渡带为0.625 0π和0.708 4π而言,图6为以切比雪夫等波纹法设计阶数N=8,wp=0.166 7π的两组滤波器镜频抑制能力比较。从图中可以看出当采样后信号基带与其两边倍频间距不等时,过渡带宽取基带与倍频间距的较大值时,其镜频抑制比采用最小值时要差,因此,针对非理想带通采样,在设计滤波器时过渡带宽应采用较小值才能更好的满足设计要求。
左右过渡带设计的滤波器镜频抑制能力对比
3 工程实现
FPGA器件是由大量逻辑宏单元构成的,通过配置这些宏单元可以形成各种硬件结构,进而构成不同的电子系统。用FPGA构成的电路可以运行于并行和串行结构,利用FPGA诸多优势可以直接完成信号处理功能[7]。文中讨论的雷达接收机采用14 bit的A/D器件,即进入数字正交检波器的信号为14 bit。总体考虑实现正交变换所包含的乘加运算以及给后续雷达信号处理传输数据所需的中间存储量之后,正交检波器实现采用Altera公司的CycloneⅡ系列EP2C35芯片。该芯片有33 216个逻辑单元,105个M4K RAM块,嵌入式18×18乘法器35个,4个PLL。
文中信号中频f0为30.416 667,当采样频率为40 MHz时,正交插值的I/Q两路系数周期很大,存储一个周期数据需要极大存储空间。经仿真发现,可以将系数中连续96个点周期扩展进行正交变换,并与采用原系数进行正交变换后的结果对比,相对误差为10-4,满足工程需要,因此可以只存储96个系数,而不必存储所有系数,从而极大地减少了系数的存储容量。
图7为I/Q正交在FPGA中实现顶层框图,I/Q两路系数分别放入只读存储器(ROM)中,即图中的coefficient＿I和coefficient＿Q模块。数据经过A/D采样后进入FPGA,在FPGA中与I/Q系数相乘,I/Q两路系数以96个点为周期,因此图中ROM容量为128个16位数据。乘法器输出的两路数据即为正交变换的I/Q两路。
正交变换之后的I/Q两路需要通过低通滤波器,在频域上滤除掉等价于高频的镜频分量,而在时域上是将正交变换后的实虚部通过FIR滤波器进行滤波。工程中滤波器可采用直接型,级联型和转置型来实现,文中所设计滤波器阶数少且系数对称,可采用FIR标准型如图8,实现。图8中滤波器的实现充分考虑了滤波器系数对称的特点,将原本需要的8个乘法器减少到了4个,加法器也有所减少,从而优化了资源配置。
滤波后数据带宽为2 B,因此回波经过正交检波之后一般还要对数据进行抽取。抽取的目的是为了降低频谱间的间隔,提高频带利用率,最后得到的就是所需信号的复包络。
4 结束语
文中以带通采样定理为基础介绍了数字正交检波器的设计。特定信号以及工艺等原因往往会造成信号中频不理想,文中对非理想的中频信号对采样以及正交变换带来的影响进行了分析,随后介绍了正交变换中滤波器的设计,对非理想中频信号滤波器系数的选择问题进行了分析,并给出选择标准,最后对工程中数字正交变换进行了设计。虽然带通采样定理具有很大的灵活性,但针对非理想中频信号,采样率的选择以及采样后对正交变换的影响是一个不可忽视的问题,应根据实际情况对其进行方针分析,确定合适的采样率、正交变换方法以及滤波器系数。文中的设计方法已成功应用于国内某新型雷达信号处理机中,并对该雷达中非理想中频信号成功地进行信号处理效果令人满意。
参考文献
[1]牛亚莉，苏涛．正交插值与脉冲压缩算法及其滤波器设计[D]．西安：西安电子科技大学，2006．
[2]徐仰彬．雷达中频采样理论的改进[J]．现代电子技术，2006，22(20)：124-126．
[3]苏涛，何学辉，吕林夏．实时信号处理系统设计[M]．西安：西安电子科技大学出版社，2006．
[4]王飞雪，王新春，雍少为，等．带通信号采样定理和全数字式正交检波器的设计[J]．电子科学学刊，1999，21(3)：207-310．
[5]丁玉美，高西全．数字信号处理[M]．西安：西安电子科技大学出版社，2001．
[6]曾涛，龙腾，毛二可．一种新的中频采样滤波器的设计与实现[J]．电子学报，2000，28(6)：51-54．
采样示波器 篇3
常用的窄带干扰抑制方法可分为时域处理和频域处理两类。时域处理方法主要是基于线性或非线性的预测滤波技术[2,3]，当干扰快速变化时，这种技术就不再适用。频域处理方法主要是利用了在变换域映射后，扩频信号和背景噪声的频谱十分平坦，而窄带干扰的频谱呈脉冲状，集中在某些频点上，可以进行快速干扰识别及滤除。但在实际应用中，如果对信号直接分块进行DFT运算，会产生严重的频谱泄露[4]，它使干扰能量在频域中扩散，影响后续的干扰抑制处理。为了减轻频谱泄露问题，大多数频域干扰抑制技术都是对分块数据先进行时域加窗再进行DFT变换。窗函数的引入，扭曲了数据块边缘处的有用信号，带来了信噪比损失的问题，且使得处理后恢复出的信号失真较大。为此，MITRE公司的CAPOZZA等人提出50%重叠加窗DFT处理算法[5]改善了信噪比插入损耗问题。后来JONES等人提出用临界采样DFT滤波器组进行频域处理抑制窄带干扰[6]。参考文献[7]对上述两种算法进行了详细的比较，认为在采样相同DFT长度的情况下，50%重叠加窗DFT算法具有更小的插入损耗，临界采样DFT滤波器组算法则具有更小的频谱泄露。为了兼顾频谱泄露抑制能力与插入损耗，本文提出了基于过采样DFT滤波器组的窄带干扰抑制方法。与CAPOZZA的加窗DFT处理方法相比，DFT滤波器组具有频谱泄露小的优点，具备更强的干扰抑制能力;而与JONES的临界采样DFT滤波器组相比，过采样DFT滤波器组能够在减轻频谱泄露的同时，实现信号的近似完全重构(NPR)，插入损耗更小。
1 DFT滤波器组
滤波器组的基本思想是将输入的全带信号，经过一组分析滤波器分解成若干个子带信号，然后针对各子带信号的特点分别进行处理，最后再通过一组综合滤波器重构出原信号[8]。合理地分析/综合滤波器组的设计，可以使得信号完全准确重构。图1为M通道临界采样滤波器组结构图。
分析滤波器组主要负责频谱分解。有一种分析滤波器hm(n)设计思路是由一个原型低通滤波器h(n)，经过线性调制，在频谱上依次移位衍生而成。即:
其中，k=0,1,…,M-1。设原型滤波器h(n)的长度为L,L为子带个数M的I倍。则经过分析滤波器组后生成的各子带信号为：
再经过M倍降采样器得到的信号为:
利用多相结构的理论继续分析，将输入信号x(n)分成M个子序列，令l=Mi+m,i=0，…，I-1,m=0,1,…,M-1,则有：
记作：
再记tm(n)=il=-Σ01xm(n-i)pm(i)，得:
这即为M点的DFT形式，类似的分析也可将组合滤波器组变换成逆DFT的形式。于是，得到M通道临界采样DFT分析滤波器组的多相结构如图2所示，其中降采样率K=M。
型低通滤波器h(n)经M倍抽取而得到。一种可行的近似完全重构的滤波器组设计方法是，分析/综合滤波器组基于同一个原型低通滤波器衍生而来，利用窗函数法设计出此原型低通滤波器[9]，方法如下：
其中w(n)为某长度为L的窗函数。综合滤波器组原型滤波器设计为g(n)=h(-n)。图2中的临界采样滤波器组，在子带混叠小与信号完全重构上难以同时兼顾。如果图2中的降采样率K
2 基于DFT滤波器组的干扰抑制方法
2.1 基于DFT滤波器组的干扰抑制模型
接收到的GNSS基带信号经过DFT分析滤波器组之后，各子带的输出即为其频谱输出，在频域进行干扰检测时，对存在干扰的谱线进行适当处理，达到抑制窄带干扰的目的，然后经过综合滤波器组，重构出滤除掉干扰后的原GNSS信号。基于DFT滤波器组的窄带干扰抑制模型如图3所示。
2.2 频域干扰抑制处理算法
在接收到的导航信号中，有用信号淹没在背景噪声下，因此在带内没有干扰的情况下，主要是高斯噪声在该频带内的频率分量。由于DFT对信号进行的是线性变换,所以DFT后的带内每一个频率分量仍是高斯分布的。设原时域高斯噪声均值为0，方差为σn2，则可算得M点DFT得到的各频率分量满足均值为0，方差为Mσn2的高斯分布，其幅度值满足Rayleigh分布，此Rayleigh分布概率密度函数与概率分布函数分别为：
计算幅度的均值为：
代入上述概率分布函数表达式可得：
则可计算出无干扰信号时的谱线幅度分布概率情况如表1所示。
由此可知，当幅度门限为3μ时，幅度小于门限的频域分量占所有频域分量的99.91%。对于慢衰落的卫星信道，由以前的噪声功率可以大致估计出当前的噪声功率。自适应估计噪声门限Th(n)的过程借鉴参考文献[10]中方法：
其中，vk(n)为DFT分析滤波器组输出的各谱线复数值，Th(n)表示估算的第n帧的背景噪声门限，Ak(n)为此帧中所有低于当前门限的谱线幅度，Sk(n)对这部分谱线进行计数。当阈值估计趋向稳定时，便得到背景噪声包络，即为门限。对幅度值超过门限的谱线，认为是干扰信号，直接置零滤除干扰。
3 性能分析与仿真比较
3.1 性能分析与对比
当临界采样DFT分析滤波器组的原型低通滤波器长度L=M且h(n)=w(n)时，此分析滤波器组就退化为加窗后的DFT;如果L=M，且h(n)≡1，则此分析滤波器组就进一步退化为单纯的DFT。其实单纯的DFT变换本身就可以看作是幅频响应为|sin(ωM/2)/sin(ω/2)|的低通滤波器衍生出的一个均匀窄带滤波器组[8],由于sinc函数旁瓣电平高,且只在圆频率为2π/M的整数倍处为0,所以滤波器组内各子带之间存在较多混叠,这也就是直接进行DFT会造成较多的频谱泄漏的原因。当干扰信号的圆频率不是2π/M的整数倍时,各子带滤波器在此处旁瓣电平都不为0,则干扰频率就会泄漏到所有的子带上。
加窗时的泄漏分布取决于所采用窗函数的频域特性，不加窗相当于使用矩形窗。矩形窗主瓣宽带窄但旁瓣电平高(最大-13.56 dB)，干扰信号很强时，其泄露出的旁瓣会在很宽的频带范围内都高于背景噪声。而对于非矩形窗，由于其旁瓣电平低，减轻了干扰信号的频谱泄露。以Blackman-Harris(4-term)窗为例，其主瓣宽度为12π/M[4]，即无论DFT长度M为多少，其主瓣都会占用约7个频点，而其最大旁瓣电平却低至-57 dB，当应对强于背景噪声50 dB的单音干扰时，Blackman-Harris(4-term)窗处理后的频谱中仅有约7个频点的电平高于背景噪声。但是窗函数的使用会带来相当的信噪比损失，损失程度可表示为式(16)[12]。仍以256点的此窗为例，信噪比损失达3.04 dB。
而对第二节中分析的DFT滤波器组方法而言，可以将其分析滤波器组看作是一个用L长的时间窗函数截取信号后做M点DFT，例如L=3M时，式(6)等效于在做了L点DFT之后，再在频域进行3:1抽取，显然这比传统加窗DFT的频谱泄露的点数会更少。而过采样滤波器组能在临界采样滤波器组基础上更进一步，既满足频谱泄露少，又能信号近似完全重构，使得信噪比处理损耗很小。
3.2 仿真结果对比
下面以2.046 MHz带宽的GPS C/A码信号进行比较实验，设定输入信噪比为-20 dB，再加入两个干信比都为70 dB的单音干扰，分别采用未加窗DFT、加窗DFT、DFT滤波器组三种方法进行窄带干扰抑制实验。加窗DFT的窗函数采用参考文献[5]推荐的Blackman-Harris窗;滤波器组的原型滤波器长度L=3 M，所用窗函数与上述相同。
图4中实线为采用相同的频域自适应估计算法分别得到的干扰门限。所加入的两个单音干扰，第一个的圆频率恰好为2π/M整数倍，后一个不是;而在现实场景中由于信号多普勒频移的不确定性以及干扰的不可知性，窄带干扰一般都不在2π/M整数倍频点上。
通过图4可以看到，直接DFT得到的频谱在第二个干扰源附近有很严重的频谱泄露，如果采用频域滤除的方法会滤除掉相当大一部分有用信号，所以这种方法极不实用。若采用加窗后的DFT，每个强单音干扰要滤除约7个频点。而采用(过采样或者临界采样)DFT滤波器组方法，每个强干扰只需滤除约3个频点。一般认为，低于25%的GNSS信号带宽被滤除后的信号导航性能都仍然能够接受[7]。因此，如果采用M=256点的DFT，则加窗DFT方法能够承受约25%×256/7≈9个这样的单音干扰，而DFT滤波器组方法能够承受约25%×256/3≈21个这样的单音干扰。由此可知，在相同DFT长度情况下，过采样和临界采样的DFT滤波器组，都比加窗DFT有更强的干扰抑制能力。
下面比较过采样和临界采样两种DFT滤波器组的处理方法，以及50%重叠加窗对信噪比的影响。把不含干扰的GNSS信号分别经历这三种不同处理流程，三种方法的DFT长度都为M,其中过采样DFT滤波器组的降/升采样率K取M/2,M=256。对通过几种常用窗函数设计出的滤波器组，分别进行对比仿真实验，得到如表2所示的插入损耗结果。可见，在采用相同窗函数设计出的原型滤波器的情况下，过采样DFT滤波器组的插入损耗比临界采样DFT滤波器组明显要小，但略高于50%重叠加窗方法。
通过仿真实验与对比分析表明，这种方法在抑制频谱泄露方面的能力，与临界采样滤波器组相同，都优于CAPOZZA提出的重叠加窗DFT处理方法[5];而在插入损耗方面，与重叠加窗DFT相当，都优于JONES等人提出的基于临界采样滤波器组的处理方法[6]。
综上所述，基于过采样DFT滤波器组的窄带干扰抑制方法兼具了上述二者的优点，兼具优越的频谱泄露抑制性能和微小的插入损耗,在GNSS抗干扰领域具有广泛的应用前景。
摘要：对于全球导航卫星系统(GNSS)接收机而言,窄带干扰十分常见并且危害较大。对此提出一种基于过采样离散傅里叶变换(DFT)滤波器组的频域窄带干扰抑制技术。与传统的基于加窗DFT处理的方法相比,这种方法能够更好地减小干扰信号频谱泄露问题;而与基于临界采样DFT滤波器组方法相比,这种方法能够更加有效地降低导航信号的畸变问题,特别适合在卫星导航接收机中应用。理论分析和仿真结果表明,基于过采样DFT滤波器组的新方法具有更强的窄带干扰抑制能力和更小的插入损耗。
关键词：全球导航卫星系统,窄带干扰,干扰抑制,DFT滤波器组,过采样
参考文献
[1]KAPLAN.GPS原理与应用[M].寇艳红,译北京:电子工业出版社,2007:181-208.
[2]POOR H V,RUSCH L A.Narrowband interference supp-ression in spread spectrum CDMA[J].IEEE Pres.ThirdQuarter,1994(8):14-27.
[3]VIJAYAN R,POOR H V.Nonlinear techniques for inter-ference suppression in spread spectrum systems.IEEETrans.Commun.,1990,38(7):1060-1065.
[4]应启珩.离散时间信号分析和处理[M].北京:清华大学出版社,2001:103-105,228-234.
[5]CAPOZZA P T,HOLLAND B J.A single-chip narrow-band frequency domain excisor for a global positioning sys-tem(GPS)receiver.Custom Integrated Circuits Conference,May 16-19,San Diego,California.
[6]JONES W W,JONES K R.Narrowband interference supp-pression using filter-bank analysis/synthesis techniques[C].IEEE MILCOM Conference,San Diego,California,Paper1992,38(1):1.
[7]RIFKIN R,VACCARO J J.Comparison of narrowbandadaptive filter technologies for GPS[R].Position Locationand Navigation Symposium,IEEE,March 2000.
[8]胡广书.现代信号处理教程[M].北京:清华大学出版社,2004:219-225.
[9]YIU K FC,GRBIC N,NORDHOLM S.Multicriteria designof oversampled uniform DFT filter banks[J].Signal Proce-ssing,IEEE,2004,11(6):541-544.
[10]TANAKA T.A direct design of oversampled perfect recon-struction FIR filter banks.Signal Processing,IEEE,2006,54(8):3011-3022.
[11]梁继业,刘会杰.卫星扩频通信中一种简化的窄带干扰抑制方法[J].电路与系统学报,2004,9(6):86-89.
采样示波器 篇4
关键词：编队导航,采样点,状态估计,正交变换
0 引言
无人机编队飞行需要编队导航系统提供飞行器之间的编队位置、编队速度和编队姿态信息。目前通常采用GPS与惯性测量装置的组合来获得位置和姿态信息。由于GPS信号容易受干扰，因此对不依赖GPS信号的编队导航系统成为近年来的研究热点。基于视觉的导航系统可以取代GPS，提供编队导航以及绝对导航信息。本文以视觉编队导航为研究背景，对确定采样型滤波器在编队导航中的应用展开研究。
编队导航系统有如下几个特点，对滤波器提出了不同的要求：
(1）系统方程中的状态包括编队姿态、编队位置、编队速度和惯性器件误差，系统状态维数为16，会出现权值为负或非局部效应。例如，采用UKF则会出现很大的负权值ω0=-13 3；而采用CKF，虽然可以避免这一情况的发生，但是由于采样点离采样中心过远，可能会存在非局部效应的问题；
(2）状态向量中，各分量之间的数值大小差异较大，在采样时应当避免对数值较小分量的忽视；
(3）采用四元数描述姿态，在采样时会遇到问题，由于四元数要求q=q02+q12+q22+q32=1，而得到的采样点中，这一条件无法满足。因此需要采用无约束的状态变量来替代四元数进行采样。
1993年，Shuster提出了一种由3个变量组成的姿态误差向量来表示姿态四元数的误差向量[1]，这3个变量可有多种表示方法，包括Gibbs向量（180°出现奇点）和修正Rodrigues参数（360°出现奇点），虽然存在奇点，但是由于是用来表示误差的，所以不会发生奇点出现的情况。本文采用无约束的Rodrigues参数来描述编队姿态[2,3]。
1 确定采样滤波器对编队导航参数的估计
1.1 编队导航运动方程
1.1.1 编队姿态运动方程
对于姿态的描述有多种描述方式，如旋转矢量、方向余弦矩阵、四元数和修正罗德里格参数。通常采用方向余弦矩阵和四元数来表示坐标系之间的转动关系，由于方向余弦矩阵有9个未知参数，而四元数只有4个，因此在姿态更新算法中采用四元数的计算量更小、精度更高[4]。当向量在不同坐标系下变换时，习惯上采用方向余弦矩阵，而方向余弦矩阵和姿态四元数之间可以相互转换，因此文中采用四元数法进行姿态更新，在坐标变换时仍然采用方向余弦矩阵来表示。
解算坐标系（b系）是指进行编队位置解算的坐标系。编队飞行中，通常选用导航坐标系或体坐标系。解算坐标系相对惯性空间的转动引起的姿态变化通过姿态运动方程来描述，可以用方向余弦矩阵或四元数来表示。
四元数的基本运算：
其中：q0为标量部分；[q1q2q3]为向量部分。
四元数具有如下性质：
(2）共轭四元数q-1=q0-iq1-jq2-kq3。
姿态四元数的微分方程为：
式中：
定义长机体坐标系l到僚机体坐标系f的方向余弦矩阵为Clf，通过Clf可以将坐标系l中的向量ul投影到坐标系f中，得到uf=Clful。
陀螺仪测得长机的惯性角速度在长机体坐标系下为ωlil,ql为惯性坐标系到长机体坐标系的姿态四元数：
陀螺仪测得僚机的惯性角速度在长机体坐标系下为ωfif,qf为惯性坐标系到僚机体坐标系的姿态四元数：
长机体坐标系到僚机体坐标系的相对四元数为：
对上式两边求导：
式（5）右侧第一项：
由姿态四元数定义：
对上式求导：
将式（2）代入上式得到：
结合式（4）和式（7）得到：
将式（6）和式（8）代入式（5）得到：
式中：
1.1.2 编队质心运动方程
编队飞行示意图如图1所示。
在僚机体坐标系下的编队位置向量定义为：
对上式两边求导得到编队速度方程：
编队加速度方程：
根据牛顿法则，惯性加速度为：
在地球坐标系中，重力加速度为：
编队加速度方程可改写为：
1.2 视觉量测方程
视觉量测系统由安装在长机上的光点和僚机上的成像装置组成，通过成像装置对光点的感测来获得量测信息。基于对效率和精度的要求，采用视线角信息的视觉导航系统，至少需要三个点。
僚机上的成像装置对长机上的光点Pi进行感测，在投影平面上得到投影点P′i。在成像空间坐标系下，Pi的位置坐标为(xic,yic,zic),P′i的位置坐标为(αi,βi,-f)。如图2所示，以P1为例，有如下几何关系成立：
视觉量测系统的输出量为[αi,βi]，根据上式得到：
在解算坐标系下，假设相机安装位置为成像空间坐标系的原点P0，坐标为(x0b,y0b,z0b)，光点Pi的位置坐标为(xib,yib,zib)，已知解算坐标系到成像空间坐标系的方向余弦矩阵为Cbc，通过坐标变换可以得到(xic,yic,zic):
光点在长机体坐标系下的安装位置(xil,yil,zil)是已知的，长机在解算坐标系下的位置坐标为(xlb,ylb,zlb)，长机体坐标系到解算坐标系的方向余弦矩阵为Clb，则(xib,yib,zib)为：
通常情况下，相机在飞行器上的安装位置是固定的，成像坐标系与飞行器体坐标系之间也是相对固定的，因此，可以假设成像坐标系与僚机体坐标系是重合的，即(x0b,y0b,z0b)=(xfb,yfb,zfb),Cbc=Cbf。由式（12）和式（13）得到：
式（14）中Cbf由僚机导航系统提供，Clf由编队姿态确定，rb=(xb,yb,zb)为解算坐标系下的编队位置。
解算坐标系为僚机体坐标系，则Cbf=Cff=I,rb=rffl=(xf,yf,zf)，由式（11）和（14）得到量测方程为：
由于式（15）的非线性很强，因此，通过单位化降低非线性强度：
1.3 确定采样滤波器设计
1.3.1 状态方程
以编队导航的编队姿态、位置和速度组成状态向量：
根据编队姿态运动方程（9）得到：
编队速度方程为：
在编队飞行中，长机和僚机的位置距离所引起的重力加速度的差异可以忽略，因此gfe=gle，由编队加速度方程（10）可得：
式中：
理想惯性器件输出根据实际输出和器件误差得到：
惯性器件误差bfa和bfg作为噪声的主要来源，在系统方程中并不是加性噪声，但是可以将其作为被估计状态变量加入到系统状态向量中：
器件误差可以看作是常值误差和随机误差的组合：
其中，ηfgu和ηfau为高斯白噪声。
1.3.2 量测方程
量测向量y由视觉量测装置的输出得到，假设有N个光标可见，那么量测向量为2N维：
假设焦距f=1，由视觉量测方程式（15）和式（16）得：
式中Pj为第j个光标的安装位置。
1.3.3 编队导航参数估计
定义误差四元数为δq=[δρTδq4]T，对应的Rodrigues参数为：
参数a的取值在0～1之间，f为比例因子。当a=0,f=1时，式（17）为Gibbs向量；当a=f=1时，式（17）为修正Rodrigues参数。对于协方差中姿态部分的小误差与任意旋转顺序下的姿态估计误差相近时，可以采用简化因子。例如，Gibbs向量线性化成半角，Rodrigues参数线性化为四元角。选择f=2(a+1)，对于小误差，δp等于旋转角。δp到δq的变换为：
通过δp完成姿态四元数的传递与更新，重新定义滤波状态向量：
k时刻，采样中心点为χ0=x̂k，对应的姿态四元数为q̂fl(0)=q̂fl,k，第i个采样点为：
对应的姿态四元数为：
δq̂fl(i)由σiδpfl根据式（18）得到：
将得到的姿态四元数代入姿态运动方程得到预测更新后的姿态四元数q̂fl,k+1 k(i)，误差姿态四元数为：
通过式（17）得到每个采样点对应的误差Rodrigues参数，将采样点中的编队位置、编队速度和惯性器件误差部分代入系统状态方程从而得到预测更新后的采样点χi。
在整个滤波过程中，姿态四元数用于编队导航状态的更新，而Rodrigues参数用于获得采样点和计算均值及协方差。
2 应用举例
选取一个固定的地理坐标系n，给出初始经度λ=38°，纬度ϕ=-77°，在该地理坐标系下长机和僚机的飞行轨迹为[5]:
长机位置轨迹：
僚机位置轨迹：
式中：僚机的旋转角速度为ωfn=[0 0ωr]T。陀螺仪的初始误差为bfg=[0.8-0.75 0.6]Tdeg/h；加速度计的初始误差为bfa=[-0.002 0.037 5-0.004]Tm/s2。光标位置见表1。
采样周期和滤波步长为0.1 s，仿真时长为300 s。考虑到计算量的问题，采用3阶采样策略，首先分别采用UKF和CKF对编队导航状态进行估计，均出现滤波发散的情况，这与前面分析的中心采样点权重过大以及非局部效应对滤波器的稳定性产生影响。因此，采用正交变换的方法对采样点进行优化，优化后的滤波器实现了收敛，仿真结果如图3～图5所示。
3 结论
通过构建僚机体坐标系下基于视觉量测的无人机编队导航系统模型，采用确定采样型滤波器对编队导航状态进行估计。采用姿态四元数来完成编队姿态的更新，由于四元数有归一化要求，而在滤波过程中，获得的采样点中的姿态部分并不能满足这一要求，因此需要采用无约束的参数来表示姿态。文中在获取采样点时以广义罗德里格参数来代替四元数，通过仿真对不同的采样策略下的编队导航参数的估计效果进行了分析。仿真结果表明，本文提出的改进方法能够满足编队导航状态估计的要求。
参考文献
[1]SHUSTER M D.Survey of attitude representations[J].Journal of Astronautical Sciences,1993,41:439-517.
[2]CHEON Y J,KIM J H.Unscented filtering in a unit quaternion space for spacecraft attitude estimation[C]//IEEE International Symposium on Industrial Electronics.Vigo,Spain:IEEE,2007:66-71.
[3]MARKLEY F L,CRASSIDIS J,CHENG Y.Nonlinear attitude filtering methods,AIAA 2005-5927[R].San Francisco,California,USA:AIAA,2005.
[4]MORTENSEN R E.Strapdown guidance error analysis[J].IEEE Transactions on Aerospace and Electroic Systems,1974,10(4):451-457.
[5]FOSBURY A M.Control and Kalman filtering for relative dynamics of a formation of uninhabited autonomous vehicles[D].Buffalo,USA:State University of New York,2006.
采样示波器 篇5
1/3倍频程谱值的计算多采用频域算法和多采样率数字滤波器法。频域算法首先对噪声采样信号进行快速傅里叶变换，计算出功率谱，然后计算每个中心频率的带宽内功率谱总和，得到1/3倍频程频谱值[1]。为对噪声信号进行准确分析，通常采用很大的FFT长度(大于或等于8 192点)来提高分析的准确性[2]。而多采样率数字滤波器算法对信号进行分析前会先进行信号抽取，减少需要处理的数据和滤波运算量，然后输入带通滤波器进行计算求出1/3倍频程频谱[3,4]。1/3倍频程分析中，有30个频带谱值，因此频域算法在硬件上比数字滤波器算法实现更简单，计算量相同的情况下运算速度比数字滤波器要快[5]。但是1/3倍频程分析中心频率为20 Hz的带宽只有4.6 Hz，在低频率段，频域算法对噪声信号的频率分辨率无法满足测量要求。只有提高采样点数才能做到，而通常噪声测量中每次测量的采样点数不会发生变化[6]。采用多采样率数字滤波器算法能够有效避免低频段声音分辨率不足的问题，减少计算量，得到较为理想的低频处的频谱分析[3]。
针对频域算法和多采样率数字滤波器算法的优缺点，采用LabVIEW设计多采样率数字滤波器组实现1/3倍频程分析，降低多采样率数字滤波器算法的设计难度，同时提高声音信号低频段1/3倍频程频分析的准确度。
1 多采样率数字滤波器原理
多采样率数字滤波器法计算噪声1/3倍频程谱值时，将20 Hz到20 kHz的声频范围划分为30个频带[7]，其原理框图如图1所示。图1中，h(n)为低通滤波器即抗混叠滤波器，↓2表示对信号进行抽取，IIR为数字带通滤波器。由图1可知，中心频率为25 Hz和31.5 Hz两频带数据经过四次抗混叠滤波和信号抽取，40～80 Hz四频带数据经过三次抗混叠滤波和信号抽取，100～160 Hz三频带数据经过两次抗混叠滤波和信号抽取，200 Hz频带的数据经过一次抗混叠滤波和信号抽取，250 Hz～20kHz二十频带数据均不经过抗混叠滤波和信号抽取，最后将各频带的数据送入数字带通滤波器后分别求均方根值(root mean square,RMS)，得到1/3倍频程谱值。
1.1 信号抽取
噪声信号分析一般采用48 kHz的采样率，要实现实时分析1/3倍频程谱值，在约为20μs的采样间隔内，需要通过30组带通滤波器，计算量太大[3]。因此通过信号抽取来减少处理的数据量，降低运算的复杂度。
若令其采样周期为Ts，将x(n)中每2点中抽取一个，依次组成一个新的序列z(n)为
经推导后得到抽取信号的频谱为
已知信号x(n)最高频率为fc，根据抽样定理可知，将x(t)抽样变为x(n)时，若fs≥2fc，则抽样结果不会发生频谱混叠。对x(n)做2倍抽取得到z(n)，若保证由z(n)重建出x(t)，则要求fs/2≥2fc，即fs≥4fc。当fs≥4fc的条件得不到满足时，信号抽取后的频谱Z(ejω)发生混叠[8]。Z(ejω)高频部分的频谱能量混叠至低频部分，使低频部分的频谱能量增加，将导致1/3倍频程谱值计算出现偏差。
1.2 抗混叠滤波
为避免信号频谱混叠，必须在信号抽取前进行抗混叠滤波即对信号进行低通滤波。设h(n)为低通滤波器，其频率响应为
设x(n)通过h(n)滤波后的输出为v(n)，其表达式为
对v(n)进行2倍抽取后的序列y(n)变为
经推导后，抽取后信号的频谱变为
经过抗混叠滤波后再进行信号抽取的频谱变化如图2所示。
如图2所示，x(n)经过抗混叠滤波器h(n)后变为v(n)，再对v(n)进行2倍抽取后变为新序列y(n)。结合图2和式(2)、式(6)可知，加上频带为(-π/2，π/2)的低通滤波器后，可以避免抽取后频谱的混迭。
1.3 数字带通滤波
由图1可知，通过数字带通滤波器后对各频带数据求均方根值的计算可得到各频带声压级。理想情况下，计算中心频率为fm(n)的1/3倍频程频带声压级要求该频带信号能量全部留在中心频率为fm(n)的频带范围内，即数字带通滤波器对该频带信号实现最大能量衰减。实际情况无法保证1/3倍频程带通滤波器能够实现对信号能量的最大衰减。由文献[7,9]可知，1/3倍频程滤波器实际衰减值如图3所示。
国标GB/T 3241—2010对1/3倍频程滤波器相对衰减限值做了明确规定，要求数字带通滤波器对信号的实际衰减限值不能超过图3所示的最小衰减限值。8阶巴特沃斯无限脉冲响应滤波器(infinite impulse response,IIR)在通频带内的频率响应曲线最大限度平坦，而在阻频带则逐渐下降为零[10,11,12]，与GBT 3241—2010要求的带通滤波器特性非常相近，适合设计1/3倍频程数字带通滤波器。
2 基于LabVIEW的滤波器设计
为实现多采样率数字滤波器法计算1/3倍频程，应用LabVIEW软件开发平台对其进行设计。
2.1 抗混叠滤波器设计
采用有限长单位冲激响应半带滤波器(finite impulse response,FIR)设计抗混叠滤波器，其归一化截止频率为0.5，通带与阻带相对于1/2奈奎斯特频率对称，因此近一半滤波器系数为零[5,13]。当抽取因子为2时，信号通过抗混叠滤波器的计算量减少。经过推导得到滤波器系数集L(n)为如下形式[14]:
采用LabVIEW的数字滤波器设计工具包设计抗混叠滤波器，抗混叠滤波器设计的主要参数有采样频率fs，通带和阻带纹波系数和阻带衰减δs[15]。根据GB/T 3241—2010对于1级滤波器频谱混叠的能量衰减规定，要求通带和阻带的纹波系数小于0.1，阻带衰减δs=80 dB。以fs=48 kHz的抗混叠滤波器设计为例，采用最优纹波法设计抗混叠滤波器，其程序框图如图4所示。
利用图4所示设计程序得到抗混叠滤波器系数集L(n)的大小为103，即该混叠滤波器为102阶。该抗混叠滤波器满足半带滤波器的设计要求，滤波器系数除L(51)=0.5，其余奇数项均为零。将系数集L(n)输入到FIR滤波器中，得到抗混叠滤波器。
由图4抗混叠滤波器设计程序得到采样频率为48 kHz的抗混叠滤波器幅频响应曲线，如图5所示。
由图5可知，该抗混叠滤波器满足通带和阻带的纹波系数小于0.1，阻带衰减δs=80 dB的设计要求。
中心频率25～200 Hz的10个频带数据分别经过的抗混叠滤波器和信号抽取次数不同，每进行一次信号抽取，抽取后采样频率变为原来的1/2。因此设计抗混叠滤波器系数时，将前一级信号抽取后的采样频率作为抗混叠滤波器设计程序的输入，而通带和阻带的纹波系数及阻带衰减保持不变，得到不同采样频率情况下的抗混叠滤波器系数，最后分别将系数输入FIR滤波器中得到抗混叠滤波器。
2.2 带通滤波器设计
2.2.1 滤波器参数计算
信号通过抗混叠滤波器和信号抽取后，采样频率发生变化，需根据不同的采样频率设计数字带通滤波器。直接利用MATLAB的FDATOOL滤波器设计工具包[16]完成对滤波器参数计算。数字带通滤波器的设计选用二阶巴特沃斯IIR滤波器，通过级联来实现。以中心频率50 Hz的1/3倍频程数字带通滤波器设计为例，根据GB/T 3241—2010的规定，中心频率为50 Hz时，其带通滤波器的上、下限截止频率分别为56.2 Hz和44.7 Hz。由图1和抽样理论可知，中心频率为50 Hz的频带数据经过3次信号抽取后，该频带数据的采样频率由48 kHz变为6 k Hz。在FDATOOL设计界面输入信号抽取后的采样频率和带通滤波器上、下限截止频率，计算得到该滤波器为4个二阶巴特沃斯IIR滤波器级联及各级IIR滤波器的前向系数和反向系数，见表1。
2.2.2 数字带通滤波器设计
选用LabVIEW中带初始条件的二阶IIR滤波器进行4级级联得到8阶巴特沃斯IIR滤波器，将表1中1/3倍频程滤波器各级IIR前向系数和反向系数按顺序输入各级带初始条件的二阶IIR滤波器中，得到中心频率为50 Hz的1/3倍频程数字带通滤波器，程序框图如图6所示。
其他中心频率点的1/3倍频程滤波器也可根据以上方法进行设计。根据国标GB/T 3241—2010对于1/3倍频程滤波器1级滤波器衰减限值的要求，所设计的数字带通滤波器的衰减值必须满足国标规定。
对25 Hz～20 kHz的1/3倍频程数字带通滤波器做相对衰减值实验验证。选取8 kHz、1 kHz、250Hz、200 Hz和25 Hz五个中心频率点进行对比。对五个中心频率点的数字滤波器做相同的4 096点的相对衰减值计算，得到各频带数据信号抽取次数M及其对应的数字滤波器衰减值，见表2。其中，“NaN”表示超出测量范围的归一化频率点，无需测量。
由表2可知，中心频率大于等于250 Hz的频带数据均不经过信号抽取，其数字带通滤波器相对衰减值在所有归一化频率点均满足国标GB/T 3241—2010对于1/3倍频程滤波器1级滤波器衰减限值的要求。中心频率小于等于200 Hz的十频带数据，若分别对其进行信号抽取的次数不足时，将导致数字带通滤波器衰减值不满足国标规定。例如，当中心频率为200 Hz的频带数据不进行信号抽取即M=0和中心频率为25 Hz的频带数据进行3次信号抽取即M=3，其数字带通滤波器的相对衰减值在归一化频率点5.392和0.185处均不满足GB/T 3241—2010对于1/3倍频程滤波器相对衰减限值的规定。分别对其提高信号抽取次数后，其相对衰减值在所有的归一化频率点满足要求。
3 系统应用
在实现噪声1/3倍频程分析的基础上利用测量麦克风，外置声卡和计算机搭建了噪声分析系统。系统组成框图如图7所示。
该噪声分析系统，不仅能够对噪声信号进行1/3倍频程分析和计算，还实现对噪声常用评价量、倍频程和功率谱的测量，实时显示和保存测量结果。根据实际需要可在该系统原有功能模块的基础上对其进行扩展。
利用噪声分析系统进行1/3倍频程分析时，首先利用测量麦克风和外置声卡采集外界声压信号，信号经过声卡的放大/衰减后送入计算机噪声分析系统，先降噪然后送入1/3倍频程分析模块进行计算得到1/3倍频程谱。最后将1/3倍频程谱值以柱状图的形式显示。
为验证噪声分析系统的可行性和稳定性，对噪声分析系统进行测试。被测信号为1 kHz纯音信号，分别在加白噪声和不加白噪声两种条件下进行测量。噪声分析系统主界面及1/3倍频程分析结果分别如图8、图9所示。
由图8和图9可知，在中心频率1 kHz频带处信号能量始终为最高。添加白噪声后，其他频带的能量均增加，这是白噪声能量和原始信号波形叠加的结果。通过1/3倍频程分析能够更加直观的了解声音信号的能量分布，得到最大声级和频带总声级。
4 结论