希尔伯特变换

希尔伯特变换（通用4篇）

希尔伯特变换 篇1

0 引言

实时全息干涉法是全息干涉技术中的主要方法之一。它能够对任意形状、任意材料和粗糙表面的三维物体进行测量,具有全场、实时、非接触、条纹对比度好、测量精度高等特点。通过实时地观察干涉条纹的变化,可以定性地分析出物体表面的形变及其应力等问题[1]。若要自动、定量地进行测量,关键是测量出反映物光场变化的干涉条纹在不同时刻的相位分布,目前对干涉条纹进行自动相位测量的方法有傅里叶变换法、相移法等[2,3]。其中傅里叶变换法适合于单幅条纹的相位分析,只能得到一幅图像中的相对变形量,不适合大尺度物体变形场的测量,而且所提取的相位值为相对相位值(即视场中一点相对于另一点的相位值),不能得出测试过程中该点的相位变化。而相移法需要参考光路中加入相移器,且需要进行精确的控制,一般情况难以实现。这些方法一般都只适合于静态干涉条纹图相位的提取。

希尔伯特变换以其适合于大范围的测量、计算时间短、高分辨力、算法简单等特点在时域信号移相、相位提取等方面得到了广泛的应用[4,5,6]。

本文利用希尔伯特变换法在时间域内对摄像机拍摄的实时全息条纹进行处理,得到视场中每一点在整个测量过程中的相位变化,从而得到全场的相位变化分布情况。

1 实时全息干涉条纹

被测物体在静止状态下完成全息图记录,全息干版H在经过冲洗处理后精确复位。如图1所示,在参考光R与物光O同时照射下,在再现光路中既有由参考光再现产生的包含原物光O′再现像的再现光场UR1,同时又有由物体直接透过全息图形成的实时光场UO0,这两个光场之间具有很好的相干性,形成实时干涉条纹。如果物体发生变化,原参考光再现的原来的物光与变化后的物光叠加的结果,将引起干涉条纹的实时变化,可表示为

式中:φ(x,y,t)=Δφ(x,y,t)-π/2,J为贝塞尔函数,K为与平均曝光量、记录材料等有关系数。当实验条件确定时,系数K、J均为常数[3]。如果物光强O和再现参考光强R近似保持不变,则叠加干涉条纹强度的变化就只取决于物光的相位变化Δφ。如果能够测量出条纹的实时相位φ,就可以根据它与引起相位变化的物理量之间的关系,求出该物理量的变化[1]。

从式(1)可以看出,实时全息干涉条纹是一种动态条纹,在某一点(即x,y确定)的光强将随时间t呈余弦关系变化,因此我们可以在时域内对实时全息干涉条纹进行处理。

2 基于希尔伯特变换的相位提取

希尔伯特变换是完全在时域中进行的一种特殊的正交变换。函数u(t)的希尔伯特变换定义为

其中Hi{u(t)}表示对u(t)作希尔伯特变换,是u(t)的线性函数,它可看作是函数u(t)与(-πt)-1的卷积:

希尔伯特变换等同于一个滤波器,一个信号经希尔伯特变换后,其幅度频谱、功率(能量)谱以及自相关函数和功率(能量)均不变。在频域内,希尔特变换的关系可表示为

式中:V(f)是函数v(t)的傅里叶变换,U(f)是函数u(t)的傅里叶变换,j是虚数符号。由式(4)可知,求一个实函数u(t)的希尔伯特变换,就是将u(t)的所有正频率成份相位移动-π/2,对所有负频率成份相位移动π/2。因此,对f(φ(t))=cos(φ(t))作希尔伯特变换得到sin(φ(t)),利用这一特性,可以得到信号的实时相位:

在Matlab信号处理工具箱中提供了现成的希尔伯特变换函数[7]。

由式(5)直接得到的相位值是随时间变化从-π/2到π/2周期性变化的实时相位,则从i时刻到p时刻的相位变化为

式中:n为条纹变化的级数,即φ(t)从π/2跳变到-π/2的次数;φ(p)和φ(i)分别表示p时刻和i时刻的实时相位。

3 实验与结果分析

本文对模拟动态条纹和实测实时全息干涉条纹进行了相位提取。在模拟实验中,视场大小为200 pixels×200 pixels,初始时刻(t=0 s)物光相位变化为0,物光的相位在41 s内按高斯曲面线性增加到如图2所示的相位分布。初始时刻视场无条纹,如图3(a)所示。在相位增加过程中,视场中间有干涉条纹冒出并不断向外移动,在t=20.48 s、t=40.96 s时刻干涉条纹如图3(b)、(c)所示,点(70,100)的灰度随时间的变化的时域信号如图4(a)所示。

为使背景强度在相位计算中带来的影响最小,本文在进行希尔伯特变换之前,先将信号经过一个高通滤波器以便去除背景强度和噪声的影响。从而得到如图4(b)所示的去除背景后的时域信号。再进行希尔伯特变换之后得到图4(c)所示的波形图,可以看出其波形的信号相位与变换前的相位相差90°。

由式(5)进行相位提取,得到如图4(d)所示相位图。此图为相位值从-π/2到π/2周期性变化的实时相位图。在程序中采用逐级计数的方法来确定条纹变化级数n,即:先设定n=0,再将前一时刻的相位值与后一时刻的相位值进行对比,当发现前一时刻的相位值减后一时刻的相位值大于阈值(实验中取2 rad)时,n加1,从而确定两个时刻间条纹变化的级数。利用式(6)可计算出该点在任意两个时刻间发生的相位变化值。

对整个观察范围内的点进行扫描,求出各点的相位变化值。从t=0到t=40.96 s过程中,各点的相位变化三维图如图5所示。从图中可以看出其相位提取结果与标准值完全相符。

通过大量实验发现,当灰度时域信号噪声增加时,相位提取结果也将出现较大误差,如图6所示为信噪比SNR=20时所得实验结果,从图中可以看出提取结果不仅存在很多随机变化的起伏,而且所得结果一般比标准值要偏大。这是因为当存在较大噪声时,程序在检测式(6)中的n(即φ(t)从π/2跳变到-π/2的次数)时,将可能受噪声影响而大于真实值。

在对实测实时全息干涉条纹提取相位时,作者利用数码摄像机拍摄了铝片受力变形时的实时全息干涉条纹的变化过程,以视频文件保存,利用Matlab软件对实时全息干涉条纹进行处理。因为数码摄像机得到的是彩色图像,在Matlab中每一帧彩色图像用3个二维矩阵分别表示红、绿、蓝三种色彩的图像数据,要处理视频文件就需要建立一个四维矩阵。考滤到实时全息干涉实验中一般采用波长为632.8 nm的He-Ne激光器作为光源,在处理时只取其中的红色图像数据,这样只要一个三维矩阵就可表示出视频文件,大大减小数据处理量。从t=0时刻到t=15.6 s的记录过程中,观察到条纹从上部向下移动,图7为t=0、t=6.4 s和t=15.6 s时刻的实时全息干涉条纹图,其尺寸大小为270×250,像素点(75,25)(以图片左下角顶点作为原坐标原点,以下边缘和左边缘分别作为x轴和y轴)的灰度随时间的变化的时域信号如图8所示。

用以上算法求得从t=0到t=6.4 s、t=11.2 s、t=15.6 s过程中,各点的相位变化三维图如图9所示。从图中可以看出:在没有条纹的部分,得到的相位变化值为0;随着时间的变化,相位变化量增加;在条纹边缘部分从现了一些毛刺现像,这是因为在条纹边缘上条纹的移动不均匀、不连续而引起的误差;在条纹的上部相位变化值比下部的相位变化值要小。这与传统的实时全息干涉条纹人工分析结果是一致的:条纹越密的地方,物光相位变化越大,条纹越疏则相位变化越小。

4 结论

利用希尔伯特变换法能有效提取低噪声情况下动态干涉条纹的相位值及其变化,在时域内对实时全息干涉条纹进行处理,可以得到视场中每一点每一时刻的相位值。进而可求出任意两个时刻间各点发生的相位变化,将为实现全自动、高精度实时全息干涉测量提供有效手段。

参考文献

[1]谢锦平,刘冬梅,王凤鹏,等.用实时全息干涉法分析岩石的侧向受力问题[J].激光杂志,2004,25(3):59-60.XIE Jing-ping,LIU Dong-mei,WANG Feng-peng,et al.Analyzing the problem of rock receive the side directions strength by the way of holographic at any time[J].Laser Journal,2004,25(3):59-60.

[2]Creach K,Schmit J.N-point spatial phase measurement techniques for non-destructive testing[J].Optics and Lasers in Engineering(S0143-8166),1996,24(6):365-379.

[3]吕晓旭,张以谟.实时全息干涉计量中时间序列条纹位相的测量方法[J].中国激光,2004,31(2):219-222.LüXiao-xu,ZHANG Yi-mo.A Novel Phase Measurement Method in Real-Time Holographic Interference Time Sequence Fringes[J].Chinese Journal of Lasers,2004,31(2):219-222.

[4]Violeta D M,Hirofumi K.Dynamic electronic speckle pattern interferometry(DESPI)phase analyses with temporal Hilbert transform[J].Optics Express(S1094-4087),2003,11(6):617-623.

[5]吕捷,王鸣,宦海,等.希尔伯特变换条纹分析法及其在非球面镜测量上的应用[J].光学学报,2005,25(6):781-785.LüJie,WANG Ming,HUAN Hai,et al.Fringe Analysis with Hilbert Transform and Its Application in the Measurement of Aspheric Mirror[J].Acta Optica Sinica,2005,25(6):781-785.

[6]高成勇,周灿林.一种基于希尔伯特变换实施相移的相位解调算法[J].光学技术,2004,30(4):508-512.GAO Cheng-yong,ZHOU Can-lin.A kind of calibration algorithm of phase-shifting based on Hilbert transform[J].Optical Technique,2004,30(4):508-512.

[7]罗军辉.MATLAB7.0在数字信号处理中的应用[M].北京:机械工业出版社,2005.LUO Jun-hui.The use of MATLAB7.0in Digital Signal Processing[M].Beijing:China Machine Press,2005.

希尔伯特变换 篇2

在许多生产机械设备中, 齿轮传动是常见的传动方式, 齿轮箱运行是否正常对整台机器的工作有较大影响, 因此, 对齿轮进行故障诊断具有重要意义。齿轮故障诊断的核心是信号处理技术, 特别是现代数字信号处理技术的发展, 使齿轮诊断故障技术步入了新的阶段。通常, 设计不当、制造不良和维护操作不善是引起振动信号上的调幅及调频现象, 这使齿轮的故障信号难以提取。然而, 希尔伯特变换在很多方面都优于其他的信号处理方法, 非常适用于非线性和非平稳过程, 具有很高的信噪比。本文把希尔伯特变换引入齿轮故障诊断, 为提取齿轮中的故障信号提供一个有效的途径。

二、希尔伯特变换及其解调机理

1. 希尔伯特变换

对于一个连续的信号S (t) 的希尔伯特变换定义为式 (1) , 即H[S (t) ]是连续信号S (t) 与的卷积。

构造解析信号z (t) =S (t) +j H[S (t) ], 它含有实部和虚部, 由希尔伯特变换构成的解析信号, 只含有正频率成分且是原信号正品率成分的2倍。

2. 希尔伯特变换解调机理

设齿轮的振动信号为Sm (t) 。

式中Am——信号振幅

am (t) =mcos (2πfrt) ——信号幅值调制函数

bm (t) =βcos (2πfrt) ———信号相位调制函数

φ0———初始相位

fr———调制频率 (故障齿轮所在轴的旋转频率)

fm——齿轮的啮合频率

将齿轮的振动信号S (t) 进行希尔伯特变换, 把幅值调制信号提出来, 可以得到式 (3) 。

由此可将两者构成解析信号zm (t) 。

Sm (t) 的包络见式 (4) , 即为求得的幅值调制信号。

振动信号的相位函数见式 (5) 。

对上式进行求导可以得到频率调制信号, 见式 (6) 。

对 (4) 和 (6) 式进行FFT (快速傅里叶变换) , 就会在频谱上很容易识别出幅值调制信号和频率调制信号的频率, 从而提取出重要的故障信息。

三、齿轮故障现场诊断实例

某石化公司芳烃厂氧化反应搅拌器G2-301结构简图及其测点布置如图1所示, 机组输入轴转速为1495r/min, 减速机各齿轮副齿数见表1。经现场测试, 发现振动值波动, 各测点的加速度值最大为3.5g, 速度值最大为5.9mm/s。

测点3处的振动值最大, 振动值超标, 说明齿轮副存在故障, 由于齿轮结构和工作原理的一些特点, 其振动信号较为复杂, 齿轮工作过程中的故障信号频率基本上表现为两个部分:一部分为齿轮啮合频率及其谐波构成的载波信号;另一部分为低频成分的幅值和相位变化所构成的调制信号。以频谱分析技术、希尔伯特变换解调技术为手段, 对其进行故障诊断分析。

在测得的功率谱中出现一级啮合频率747.5Hz, 二级啮合频率220.5Hz及高次谐波440Hz、660Hz, 三级啮合频率相对较小, 各级啮合频率见表2。从频谱图中看到二级啮合频率及其高次谐波成分在频谱中占主要成分, 幅值较大, 明显高于其他频率成分, 有丰富的边频带成分, 且分布较宽, 这些频带是由所在齿轮轴调制而成, 由此可以判定二级啮合齿轮存在磨损故障。

对故障进一步分析, 220.5Hz两边分布的边频带进行分析, 细化后的边频带结构, 寻找故障的特征信息。利用希尔伯特变换解调原理, 从测点3测得的解调谱中 (图2) 可以看到:主要调制频率为Ⅲ轴的转频 (3.5Hz) , 二次谐波2 (7Hz) , 三次谐波3 (10.5Hz) , 幅值相对较大, 由此可以判断出测点3的振动信号调制来源于中间轴Ⅲ, 即该轴的大齿轮发生故障。

四、结论

讨论希尔伯特变换在齿轮故障中的算法并加以验证, 在信号分析中, 它不会造成信号能量的扩散和泄露, 是自适应的信号处理方法。频谱分析是齿轮故障信息的最基本的研究方法, 但仅靠频谱分析不能有效判断齿轮故障的程度和部位。本文在解决齿轮故障问题上利用了希尔伯特变换解调分析技术, 该算法优于其他信号处理方法, 可准确地提取故障信息, 是一种行之有效的齿轮故障诊断方法。

摘要：信号处理技术是诊断齿轮故障的核心, 通过对一台搅拌器齿轮故障的成功诊断实例, 将希尔伯特变换引入其中, 阐述希尔伯特变换的定义和解调原理。该方法能够及时、准确地提取故障有效信息, 是一种可靠的故障诊断方法。

关键词：希尔伯特变换,齿轮故障,解调

参考文献

[1] 杨国安.机械设备故障诊断实用技术[M].北京:中国石化工业出版社, 2007

[2] 廖伯瑜.机械故障诊断基础[M].北京:冶金工业出版社, 2007

[3] 钟佑明, 秦树人, 汤包平.一种振动信号新变换法的研究[J].振动工程学报, 2002

希尔伯特变换 篇3

稳态量的各种选线方法在工程实际上并不实用,接地故障时暂态电流比稳态电流大很多,基于暂态分量的选线方法选线精度高、效果好且不受消弧线圈的影响[1,2]。但是由于电流系统的结构和运行状态不尽相同[3,4],而导致系统发生单相接地故障后,暂态信号不明显或不稳定[5,6],如果只使用单一特征分量的选线方法不可能彻底解决小电流接地系统的选线难题[7],也不可能做到100%正确选线。所以,必须通过对各种选线方法进行有机综合,才能够建立出完善的故障选线方案。

本文分析故障暂态零序电流特征,运用希尔伯特变换求取瞬时频率,然后进行小波分解提取特征频带计算频带能比因子,根据计算所得的频带能比因子k判断接地电流的成分,对于高频分量占主要成分的情况,采用故障线路暂态能量最大原理进行故障选线;对于低频分量占主要成分的情况,则采用直流分量法,因故障线路零序电流中含有较强的直流分量,而非故障线路中几乎不含直流分量,因此可以根据这个差异性进行选线。该综合选线方法可以自适应地选定小波分解层数提取信号特征,减少了算法时间,为用于实际创造了条件。

1 单相接地故障暂态零序电流特征分析

发生单相接地故障时小电流接地系统的结构图[8]如图1 所示,Q处于断开状态时为中性点不接地系统,Q处于闭合状态时为中性点经消弧线圈接地系统。图中,L0i、R0i、C0i分别为第i条线路的电感、电阻和分布电容;u0为母线零序电压;L为消弧线圈零序电感;R为消弧线圈串联电阻;Rf0为零序过渡电阻;Uf0为故障接地点位置等效电压;n为馈出线数。

令u0Uφmsin(t ),则接地电流的暂态分量为

式中:w0为回路的自振荡角频率;L0为零序回路中的等效电感;为零序电压的初相角;;tL、tC分别为电感和电容回路的时间常数; d 为自由振荡分量的衰减系数,

由式(1)可知,在相电压初相角为不同值时, 接地电流中的电容电流和电感电流差别很大,暂态电感电流由直流分量和工频分量构成[9]。暂态容性电流频率为暂态电流中的高频段分量,而暂态感性电流频率为低频段分量[10]。因此相电压初相角不同时,暂态电流信号中的低频段分量与高频带分量的比例差别较大,为反映这种差别,引入频带能比因子,它是暂态电流信号中的低频段分量与高频段分量的能量之比,即

其中:εL为暂态电流中低频成分的能量;εH为暂态电流中高频成分的能量;xL为小波包分解低频段的小波包分解系数;xH为小波包分解高频段的小波包分解系数。

为保证选线的正确性,本文根据频带能比因子是否大于1 来判断零序电流中的主要成分,然后采用不同的方法进行选线。当零序电流的主要成分是高频分量时,采用暂态方法进行选线;当零序电流的主要成分是低频分量时,采用暂态电流分量法, 利用小波重构提取衰减直流分量,比较衰减直流分量的能量与基频信号的能量比进行选线。

2 选线原理

2.1 希尔伯特变换

一个实值函数的希尔波特变换是将信号s ( t)与1 /(πt) 做卷积,以得到s'(t),其表达式为

其中,h(t) =1/ (πt)。

则其反变换为

易知,s (t) 与s'(t)是一对共轭复数,则假设

则由相位公式可知,瞬时频率为

2.2 希尔伯特算法流程

将小波包分解作为希尔伯特变换以及特征频带提取的预处理,通过逐层对零序电流进行小波分解,分解成高、低频信号,运用希尔伯特变换求取高、低频信号的瞬时频率,根据求得的频率判断是否符合我们所需的特征频带,如果符合,终止小波分解;如果不符合,进行下一层的小波分解,直到满足我们所需的要求。我们设定300~1500 Hz为我们所需的高频频带,0~50 Hz为我们所需的低频频带,其算法流程如图2 所示。

假设零序电流i (t) 的第i层第k个子频带分量di ,k(t ) 的瞬时频率为fi ,k(t ) ,其中t =1,···, N 。设( fA, fB) 为特征频带的频率,则提取条件为

式中是瞬时频率的平均值。

2.3 选线原理

采用db10 小波对零序电流进行小波分解,采样频率为4 000 Hz,运用希尔伯特变换求取瞬时频率, 当满足

时,停止小波高频分解,当满足

时,停止小波低频分解。然后求取频带能比因子, 如果频带能比因子最小值小于1,则认为零序电流主要成分是高频分量时,采用小波能量极大值法进行选线,比较小波高频能量,找出能量最大的那条出线,再将最大能量与其他线路能量进行比较;当最大能量大于等于其他能量之和时,则认定能量最大线路j发生单相接地故障。反之,如果所有线路的暂态高频能量相差不大,则认为是母线故障。如果每条线路的频带能比因子都大于1,则暂态电流中的主要成分是低频分量,对零序电流进行小波分解到一定层数提取相应的直流分量和工频分量,对这两部分信号进行单枝重构,然后比较重构信号的能量,即直流分量的能量与工频信号的能量比,令,则故障线路m1,非故障线路m1,如果所有线路的m都小于1,则为母线故障。

实现步骤的流程如图3 所示。

3 仿真结果及分析

假设系统为一个有4条线路的10 kV配电网络。 线路仿真参数设置为

线路1~线路4 的长度分别为6 km、9 km、12km、20 km,低压侧中性点开关断开时系统不接地, 开关闭合时系统经消弧线圈接地, L计算公式为

其中:l为系统线路总长;fN为工频50 Hz;c0为馈线零序电容;v为失谐度,v (ICIL) IC。

在上述仿真条件下,改变相电压初相角以及故障时刻,则频带能比因子将是不同的,不同初相角和故障时刻下各出线的频带能比因子如表1 所示。

由表1 可知,频带能比因子与电压初相角有关。 在Φ=0o附近时,故障线路和非故障线路的频带能比因子都大于1,说明低频能量比高频能量多,此时低频能量占主要成分,采用暂态直流分量法选线。 在Φ=90o附近时,故障线路和非故障线路的频带能比因子都小于1,说明高频能量比低频能量多, 此时高频能量占主要成分,则采用高频能量极大值方法选线。在Φ=0~90之间时,故障线路和非故障线路的频带能比因子不是都大于1 或者不是都小于1,则最小的频带能比因子小于1 时,认为暂态电流中高频能量多,也采用高频能量极大值方法。

对图1 所示系统在Matlab中搭建模型,在不同线路长度、补偿度、过渡电阻、故障点进行大量仿真,结果如表2 所示。

由表2 可知,对下述几种情况采用高频能量极大值法,如表3 所示,都满足, 因此可以正确选线。

表2 剩余情况采用暂态电流分量法,如表4 所示,满足故障线路m1,非故障线路m1,如果都小于1,则认为是母线故障,因此也可以正确选线。

m同时还受过渡电阻的影响,假设故障发生在中性点经消弧线圈线路4 的10 km处,采用过补偿8%,通过改变过渡电阻的值,可以得到如图4 的变化曲线。由图可见,随着过渡电阻阻值的增加,故障线路的m减小,当过渡电阻为10 k 时,m还是大于1,同非故障线路的m相比还是要大很多, 这表明用该方法适用性强。

综上分析,基于希尔伯特变换的综合选线方法不受接地方式、过渡电阻、相电压初相角的影响, 可靠性高,能够准确选出故障线路和母线。

4 结论

希尔伯特变换 篇4

通信系统中,经常需要对一个信号进行正交分解,即分解为同相分量和正交分量。由于希尔伯特变换可以提供90°的相位变化而不影响频谱分量的幅度,即对信号进行希尔伯特变换就相当于对该信号进行正交移相,使它成为自身的正交对[1]。因此,希尔伯特变换在通信领域获得了广泛应用。

在传统的设计中,希尔伯特变换器可由一个FIR 滤波器和一个时延模块实现,也可由一组滤波器对实现,而实现FIR型希尔伯特变换器的一个简单方法就是对原型低通滤波器作正弦/ 余弦变换[2]。但是,无论哪种方法都需要通过计算对低通滤波器的系数进行转换,其计算繁琐且存在一定的误差。Matlab作为滤波器设计的基础软件,不仅可以快速有效地实现希尔伯特变换器的设计、分析仿真和最优化,而且可以直接计算出希尔伯特变换器的系数,加之Matlab具有强大的接口功能,为后续的设计提供了方便。

1 希尔伯特变换器的基本原理

连续时间信号x(t)的希尔伯特变换(t)定义为[3]:

$\begin{array}{l}\widehat{x}(t)=\frac{1}{\text{π}}{\displaystyle {\int }_{-\infty }^{+\infty }\frac{x(\tau )}{t-\tau }}\text{d}\tau =\frac{1}{\text{π}}{\displaystyle {\int }_{-\infty }^{+\infty }\frac{x(t-\tau )}{\tau }}\text{d}\tau \\ =x(t)\otimes \frac{1}{\text{π}t}(1)\end{array}$

由式(1)可得单位冲击响应h(t)=1/(πt),由于jh(t)=j/(πt)的傅里叶变换是符号函数sgn(ω),所以希尔伯特变换器的频率特性为:

${\rm H}(\text{e}{}^{\text{j}\omega })=-\text{j}\text{s}\text{g}\text{n}(\omega )=\{\begin{array}{l}-\text{j},\omega >0\\ \text{j},\omega <0\end{array}(2)$

记${\rm H}(\text{j}\omega )=|{\rm H}(\text{j}\omega )|\text{e}{}^{\text{j}\phi (\omega )}$,当$|{\rm H}(\text{j}\omega )|=1$时:

$\phi (\omega )=\{\begin{array}{l}-\text{π}/2,\omega >0\\ \text{π}/2,\omega <0\end{array}(3)$

信号x(t)的希尔伯特变换可以看成是信号x(t)通过一个幅度为1的全通滤波器输出,信号通过希尔伯特变换器后,其负频率成分作+90°的相移,而正频率成分作-90°的相移。

这类滤波器要求滤波器的零频响应为0,若滤波器阶数为偶数,则还要求Nyquist频率(归一化频率为1)处的响应为0。即如果滤波器的阶数为偶数,那么增益在频率为0 Hz和fs/2处必须降为零,希尔伯特滤波器必须是一个带通滤波器。如果滤波器的阶数为奇数,那么增益在频率为0 Hz处必须降为零,希尔伯特滤波器必须是一个高通滤波器[4]。

2 希尔伯特变换器的Matlab设计

2.1 直接程序法

Matlab信号处理工具箱提供了firls函数和remez函数,它们的调用格式语法规则相同,只是优化算法不同,函数firls利用最小二乘法使期望的频率响应和实际的频率响应间的误差最小;函数remez实现Park-McClellan算法,这种算法利用remez交换算法和Chebyshev近似理论设计滤波器,使实际频率响应拟合期望频率响应达到最优[5]。

函数调用格式为b=remez(n,f,m,′h′)或b=firls(n,f,m,′h′)[6],其中,n为滤波器的阶数;f为滤波器期望频率特性的频率向量标准化频率,取值0～1,是递增向量,允许定义重复频点;m为滤波器期望频率特性的幅值向量,向量m和f必须同长度且为偶数;b为函数返回的滤波器系数,长度为n+1,本文将采用remez函数法。

下面设计一个希尔伯特变换器,要求采样频率为2 000 Hz,通频带为50～950 Hz,滤波器阶数为60阶。实现程序如下:

clf;n=60;

f=[0.05 0.95];m=[1]; % 理想滤波器的幅频特性

fs=2000; % 采样频率

b=remez(n,f,m,′h′); % 采用remez设计Hilbert变换器

[h,w]=freqz(b,1,512,fs); % 计算Hilbert变换器的脉冲响应

figure(1);

plot(w,20*log10(abs(h)));grid; %计算幅频特性

axis([0 1000 -40 10]); %绘制正半轴频谱

xlabel(′频率/Hz′);

ylabel(′幅值/dB′);

title(′Hilbert变换器的幅频特性′);

figure(2);

plot(w,angle(h)/pi*180); % 计算相频特性

grid;

xlabel(′频率/Hz′);

ylabel(′相角/′);

title(′Hilbert变换器的相频特性′)

设计的希尔伯特变换器的特性如图1,图2所示。

从仿真结果可以观察到增益在0 Hz和1 000 Hz处降为零,即为带通滤波器;同时具有严格的线性相位特性,符合设计要求。在设计中如果特性不满足要求,原有的参数必须作相应的调整,在程序中只需对参数进行重新设定,就可以得到所需要的希尔伯特变换器。

2.2 利用FDATool工具设计法

FDATool是Matlab信号处理工具箱专用的滤波器设计分析工具,操作简单、灵活,可以采用多种方法设计不同的滤波器,同时可以实现滤波器的最小阶数设计[7]。在Matlab命令窗口输入FDATool后回车就会弹出FDATool界面。

根据2.1中的设计实例,首先在Filter Type栏中选择Hilbert Transformer,在Design Method栏中选择Equiripple法,在filter order中选择60,在Frequency and Magnitude Specifications中设置F=[50 950];M=[1];Fs=2 000,最后点击Design Filter,通过菜单选项Analysis可以在特性显示区看到滤波器的各种特性,如图3～图5所示。

在幅频特性和相频特性满足要求的同时,由图3可知单位脉冲响应为奇对称,即h(n)=-h(N-n-1),也符合希尔伯特变换器的特性。若设计不满足要求,则可以直接在FDATool界面中改变参数,在设计满足要求后,还可以把希尔伯特变换器的系数导出为Matlab变量,文本文件或C语言头文件等[8],这为后续的设计提供方便。

3 希尔伯特变换器的效果验证

3.1 直接程序验证法

对于所设计的希尔伯特变换器,频率为10 Hz的振动作为输入信号,采样频率为100 Hz,验证对应的数据点是否满足相位相差90°的特点,主要实现程序如下:

clf;n=60;

f=[0 0.05 0.95 1];a=[0 1 1 0]; % 理想滤波器的幅频特性

b=remez(n,f,a,′h′); % 采用remez 设计Hilbert变换器

[h,w]=freqz(b,1,512); % 计算滤波器的脉冲响应

t=0:1/100:3;x=sin(2*pi*10*t); % 时间序列和输入信号

figure(1)

subplot(2,1,1),plot(t(1:100),x(1:100))

title(′输入信号′)

y=filter(b,1,x); % 运用remez设计的Hilbert变换器进行滤波

subplot(2,1,2),

plot(t([1:100]+20/2),y([1:100]+20/2),′r-′) % 绘制与输入信号对应的输出信号,考虑了延迟效应[9]

xlabel(′时间/s′)

title(′输出信号′)

由图6可知,该希尔伯特变换器对阻带和通带波纹进行了控制,输出信号的相位比输入信号的相位前移了90°,符合希尔伯特变换的性质。如若把输入信号的频率变为-10 Hz,则输出信号较之输入信号会后移90°,如图7所示。通过验证可知,该希尔伯特变换器实现了π/2移相。

3.2 Simulink仿真法

利用Matlab提供的Simulink工具对该希尔伯特变换器进行仿真,输入信号以3.1节中的要求为例,仿真系统如图8所示,在建立仿真系统时,注意对所选择的模块进行相应的参数设计。仿真结果如图9,图10所示。

4 结 语

在Matlab平台上,采用直接程序法和FDATool工具法均可以快捷有效地完成希尔伯特变换器的设计,可以随时对比设计要求和希尔伯特变换器的特性,以使设计达到最优化。由于Matlab具有强大的接口功能,设计的结果可以很方便地移植到DSP,FPGA等器件中[10]。在实际使用中,只需按要求修改参数,即可实现不同的希尔伯特变换器,实用性较强。

摘要：在通信系统中,希尔伯特变换是被广泛应用的重要变换。为了实现数字解调,通常需要借助希尔伯特变换器对信号进行分解,利用Matlab设计希尔伯特变换器是一种最为快捷、有效的方法。通过具体的设计、仿真及对原始信号和经过希尔伯特变换器输出延迟信号的比较,说明Matlab是一个在滤波器设计方面很有力的工具。

关键词：Matlab,Hilbert变换器,remez,FDATool工具

参考文献

[1]汪璇,曹万强.Hilbert变换及其基本性质分析[J].湖北大学学报:自然科学版,2008,30(1):53-55.

[2]张兵,汤清华,吴国安,等.基于DSP的数字正交解调器的设计与实现[J].微电子学与计算机,2009,26(7):32-35.

[3]王薇,王晓茹,黄晓青,等.无功功率测量的Hilbert数字滤波器研究[J].电测与仪表,2007,44(3):9-12.

[4]陈啸晴.基于Hilbert数字滤波的无功功率表设计研究[D].长沙:中南大学,2009.

[5]刘波,文忠,曾涯.Matlab信号处理[M].北京:电子工业出版社,2006.

[6]张照明,刘郑波,刘斌.应用Matlab实现信号分析和处理[M].北京:科学出版社,2006.

[7]飞思科技产品研发中心.Matlab 7基础与提高[M].北京:电子工业出版社,2007.

[8]任淑萍,王欣峰.基于Matlab的数字滤波器设计与仿真[J].机械管理开发,2009,24(3):181-183.

[9]万永革.数字信号处理的Matlab实现[M].北京:科学出版社,2007.