信号周期范文

信号周期范文（精选6篇）

信号周期 第1篇

频谱分析在数字信号处理中用途广泛:如滤波、检测等方面, 这些都需要DFT (Discrete Fourier Transform) 运算[1,2,3]。信号的Fourier变换建立了信号的时域与频域之间的一一对应关系, 可以借助DFT来分析。有限长序列的DFT可以由数字方法直接计算, 且DFT存在快速算法, 便于用计算机处理[4,5,6]。本文介绍具体的连续非周期信号结合Matlab软件来分析其频谱。

1 连续非周期信号谱DFT分析

在已知连续信号数学解析式的情况下, 非周期信号的频谱可以根据Fourier变换的定义进行解析计算。实际应用中的多数信号不存在数学解析式, 信号的频谱无法利用傅里叶分析公式方法直接计算, 一般需采用数值方法进行近似计算分析频谱, 在进行数字计算时, 需对计算的连续变量进行离散化。由于连续非周期信号x (t) 的频谱函数X (jω) 是连续函数, 因此, 需要对其进行离散化处理得到x[n]以近似分析相应的频谱。通过建立序列x[n]的离散Fourier变换X[m]与连续非周期信号x (t) 的Fourier变换X (jω) 之间的关系, 可以利用DFT对连续非周期信号频谱进行近似分析, 此近似分析过程中一般将会出现三种现象:混叠现象、泄漏现象和栅栏现象[7,8]。这些现象与应用中信号和DFT的参数选择有关。下面分别讨论近似过程中可能出现的问题及其解决方法。

2 DFT分析过程中出现的若干问题

2.1 混叠现象

由DFT计算出的频谱是信号x (t) 的频谱X (jω) 周期化的抽样值, 如果连续信号不是带限信号, 或者抽样频率不满足抽样定理, 在连续信号离散化时, 就会出现信号频谱的混叠。解决连续信号离散化过程中的频谱混叠主要有两种方法:对于带限连续信号, 只要提高抽样频率使之满足时域抽样定理;对于非带限连续信号, 可根据实际信号对其进行低通滤波, 使之成为带限信号。工程实际中的连续信号一般都不是带限信号, 连续信号在抽样前通常都经过一个模拟低通滤波器 (称为抗混叠滤波器) 进行低通滤波, 以减少混叠误差, 提高频谱分析精度。

2.2 泄漏现象

对连续非周期信号的采样序列x[n]进行DFT运算时, 时间长度总是取有限值, 在将信号截短即时域加窗处理的过程中, 出现了分散的扩展谱线的现象, 称为频谱泄漏。对离散序列的加窗实际上是将离散序列与窗函数相乘, 加窗后信号的频谱是加窗前信号的频谱与窗函数频谱的卷积, 造成截短后信号的频谱与截短前信号的频谱不同, 所得的频谱在原来没有频谱的区间出现了频谱。原来比较尖锐的谱峰变得比较平缓, 当两个不同频率的谱峰靠得比较近时, 可能显现不出两个明显的峰值。特别是强信号谱的旁瓣可能淹没弱信号的主谱或误认为是另一假信号的主谱线。矩形窗的旁瓣幅度大, 谱间干扰严重。频谱泄漏使频谱变模糊, 分辨率 (事实上通常规定DFT的频率分辨率为fsN, fs为采样频率, N是指信号x[n]的有效长度) 变差, 泄漏程度与窗函数幅度谱主瓣宽度有关。窗型一定, 窗口越长, 主瓣越窄, 频谱泄漏越小;窗口长度一定, 矩形窗主瓣最窄, 频谱泄漏最小, 但其旁瓣的幅度最大。因此为了尽量减少泄漏现象, 应选用旁瓣幅度小、主瓣窄, 即“泄漏”小的窗函数。相对而言, 布莱克曼窗的旁瓣幅度比矩形窗小, 谱间干扰小, 但其主瓣过渡带宽, 分辨率差。采样频率或采样周期是在满足混叠误差前提下选取的, 当采样频率或采样周期确定后, 适当增加窗口长度有利于减小泄漏误差。

2.3 栅栏现象

DFT得到的频谱X[m]只能是连续非周期信号频谱X (jω) 上的有限离散频点采样, 由于X[m]是离散序列, 因而无法反映抽样点之间的细节, 就如同隔着百叶窗观察窗外的景色, 这种现象称为栅栏现象。栅栏现象是利用DFT分析连续非周期信号频谱过程中无法克服的现象, 有时频谱中的某些重要信息恰好就在抽样点之间, 将被错过, 而检测不出。为了改善栅栏现象, 把被“栅栏”挡住的频谱分量检测出来, 可在原记录序列后面补零, 增加DFT的长度, 即增加频域X (jω) 上的采样点数N, 改变离散谱线的分布, 就可能检测出原来看不到的频谱分量。

3 利用DFT进行谱分析的参数选择

在利用DFT分析连续时间信号的频谱时, 涉及频谱混叠、频率泄漏及栅栏现象。频率混叠与连续信号的时域抽样间隔有关, 频率泄漏与信号的时域加窗截短的长度及窗型有关, 栅栏现象与DFT的点数有关。在大多数情况下, 一般已知待分析连续信号的最高频率, 以及希望的DFT分析的频率分辨率。下面根据信号的Fourier变换的理论, 讨论利用DFT进行谱分析的参数 (抽样频率、持续时间、样点数等) 选择的原则。

首先确定信号抽样频率fs, fs应满足时域抽样定理, 即fs2fm, 其中fm为待分析的连续信号的最高频率, 抽样间隔T应满足:

然后确定抽样信号的长度N, N应满足频率分辨率Δf的要求, 即:

矩形窗时取c=1, Hamming窗时取c=2。

根据谱线间隔Δfd确定DFT的点数L, 即:

L一般取满足式 (3) 的2的整数幂次。

4 结合实例分析[9,10]

实例1:已知一连续信号为x (t) =cos (2πf0t) +cos (2πf1t) , 其中f0=100 Hz, f1=130 Hz。现以频率fs=600 Hz对该信号进行抽样, 试利用DFT分析其频谱。

由于抽样频率fs大于信号x (t) 的最高频率f1的2倍, 故抽样过程没有造成混叠, 抽样后的序列

由于x[n]为无限长序列, 可采用矩形窗对其进行加窗截短处理。为了能够分辨这两个间隔为Δf=f1-f0=30 Hz的相邻谱峰, 由式 (2) 可得矩形窗的长度N应满足N≥20。现应用Matlab软件进行频谱分析。分别取N=10和N=20时, 由DFT计算出的频谱做对比, 如图1所示。

程序如下:

由于N=10不满足式 (2) , 所以由由图1 (a) 几乎分辨不出信号中两个频率分量。而由图1 (b) 可见, 当信号的长度满足式 (2) 时, 可清晰地分辨出信号中的两个不同频率分量。

实例2:已知一连续信号为y (t) =cos (2πf0t) +0.15 cos (2πf1t) , 其中f0=50Hz, f1=100 Hz。现以频率fs=400 Hz对该信号进行抽样, 试利用DFT分析其频谱。采样频率满足抽样定理, 所以采样后的频谱不会产生混叠现象, 但由于信号y (t) 中存在一个较弱的频率分量f1, 若采用矩形窗函数加窗, 则由于其旁瓣泄漏较大, 很难检测出信号y (t) 幅度较小的频率分量f1, 因而采用Hamming窗函数, 如图2所示。

其实现程序如下:

图2 (c) , 图2 (d) 与图2 (a) , 图2 (b) 相比, Hamming窗以增加主瓣宽度来降低旁瓣能量, 检测出较弱的频率分量f1。在图2 (c) , 图2 (d) 中当采用同一种窗时, 截取信号长度N=30个样点时, 由于所取的数据点较少, 频率分辨率低, 由DFT计算出的频谱仍不能显示幅度较小的频率分量。当截取信号N=50个样点时, 能够清楚地显示出幅度较小的频率分量。

栅栏现象是利用DFT分析连续频谱过程中无法克服的现象, 有时频谱中的某些重要信息恰好就在抽样点之间。为了改善栅栏现象, 以观察到频谱中更多的细节, 常用的方法是在截短后的序列后补0, 构成一个L>N的序列, 序列补0后对应的频谱函数不变, 但对序列进行L点的DFT时, 两个相邻抽样点之间的谱线间隔Δfd就变为fsL, 因为抽样频率fs不变, L>N, 所以图1和图2是将N点截短信号补0后L=512点可减小频谱分析时的谱线间隔, 从而所计算出频谱将会显示更多的细节。

图3分别给出N=20, L=32, 64, 128, 256时, 利用矩形窗计算出的实例1中的信号的频谱。

从图3中可以看出, 随着补0的增加, 可以更多地显示出频谱的信息, 从而减少了栅栏现象。

需要说明一点通过补0只能提高信号频率的显示分辨率。换句话说, 如果连续信号在离散化或时域加窗过程中, 由于混叠或泄漏等过程已经造成信号频谱中信息的失真, 则无论怎么补0也无法再恢复已损失的信息。

5 结语

本文利用Matlab实现连续信号DFT频谱的计算。可见利用Matlab软件分析频谱方便、快捷, 减少了计算量。对连续非周期信号先要离散化, 然后截短处理, 最后做DFT分析。在这些近似分析过程中可能会出现三种现象, 避免及减少这些问题, 只有合理选择DFT的一些参数, 可以完全使分析结果在工程误差允许范围内。

参考文献

[1]李莉.数字信号处理原理和算法实现[M].北京:清华大学出版社, 2010.

[2]汉泽西, 姚英彪.用DFT分析正弦信号频谱时应注意的几个问题[J].西安石油学院学报, 2003, 18 (2) :67-70.

[3]井敏英, 郭佳林, 李大伟, 等.连续时间信号频谱分析研究及Matlab实现[J].大众科技, 2011 (4) :42-43.

[4]陈怀琛.数字信号处理教程:Matlab释义与实现[M].北京:电子工业出版社, 2004.

[5]林爱英, 滕红丽, 袁超, 等.DFT在信号谱分析中的应用[J].安徽工业大学学报, 20011, 28 (2) :192-196.

[6]吴镇扬.数字信号处理[M].2版.北京:高等教育出版社, 2010.

[7]陈后金.数字信号处理[M].2版.北京:高等教育出版社, 2008.

[8]黄寒华.关于滑动DFT算法中的频谱泄漏问题的探讨[J].现代电子技术, 2007, 30 (24) :202-204.

[9]徐明远, 刘增力.Matlab仿真在信号处理中的应用[M].西安:西安电子科技大学出版社, 2007.

信号周期 第2篇

数字信号处理, 是从20世纪50年代起, 随着信息学科和计算机学科的高速发展而迅速发展起来的一门学科, 它的重要性, 在各个领域的发展和应用中日益表现出来, 数字信号处理技术已成为每一个电子信息科学工作者必须掌握的重要知识。数字信号处理课程是国内外大学信息、通信、航天、地震、天文、声学、生物医学工程、语音处理、机械工程等专业本科生与研究生的必修专业课程。数字信号处理是把信号用数字或符号表示成序列, 通过计算机或通用 ( 专用) 信号处理设备, 用数值计算方法进行处理, 提取有用信息, 或将信号按照需要进行改善, 以便于应用。例如, 滤波、检测、变换、增强、估计、识别、调制/解调、编码/解码、参数提取、频谱分析等, 被广泛地应用于国民经济和社会生活的各个方面。

但是由于《数字信号处理》是一门用数学工具来分析解决典型物理问题的课程, 所以要求学生掌握高等数学和高等代数等相关的数学知识, 对学生的数学功底要求很高。数字信号处理的理论、算法比较抽象, 存在公式多、推导复杂、抽象难懂等特点, 学生理解和接受起来有一定的困难。而MATLAB具有算法编程功能强大和图形绘制功能丰富的特点, 所以把数字信号处理的理论教学和MATLAB的应用结合起来, 改革课堂教学模式, 利用MATLAB进行计算机仿真辅助教学, 使学生加深对基本概念和基本理论的理解, 增进学习的深度和广度, 培养学生独立思考的习惯, 锻炼学生的动手能力, 使学生在这门课程中系统地理解和掌握数字信号处理的知识。

二、MATLAB 语言的特点

MATLAB是矩阵实验室 ( Matrix Laboratory) 的简称, 是美国Math Works公司出品的商业数学软件, 主要包括MATLAB和Simulink两大部分, 它将计算、可视化和编程功能集成在非常便于使用的环境中, 是一个交互式的、以矩阵计算为基础的科学和工程计算软件。

MATLAB是一种应用编程软件, 编程语言简洁明了, 易于掌握, 通过运用高级语言编程思想解决了信息中的问题, 然后借助丰富的图形绘制功能做出仿真结果, 实现所见即所得的功能, 为信号的分析与处理提供了强有力的辅助工具。

三、MATLAB 在数字信号课程中周期方波信号的频率分解的辅助教学

在数字信号处理中, 一个周期为2m的函数, 其傅里叶级数展开公式为:

在公式 ( 1) 中:

在区间 ( 0, 2m) 里等间隔地取N个取值点, 采样时间间隔设置为Δt, 其中m = NΔt/2, 则其离散傅里叶变换形式为:

在公式 ( 2) 中:

由上述可知, 在数字信号处理的教学过程中, 涉及到大量的利用离散傅里叶变换、快速傅里叶变换等进行的数学运算和推导, 用于频率分析和滤波。如果仅仅用公式计算来表达复合信号的分离过程, 会显得枯燥无味, 很难引起学生的兴趣, 导致教学效果不佳, 学习效率低下。

MATLAB具有强大的编程性能, 可以解决上述计算的编程需要, 利用MATLAB自带的函数, 编写实现离散傅里叶变换的算法, 再通过MATLAB强大的图形功能绘出二维图像, 完成对频谱分析问题的模拟, 能生动形象地解决数学问题, 增强教学效果。

周期方波信号是一种非常重要的周期信号, 它实际上是由不同频率和不同振幅的周期正弦信号叠加而成, 利用离散傅里叶变换能从方波信号中分离出原来的正弦信号。

将振幅为1. 5的方波进行Fourier分析, 研究能否从中分析出含有不同频率不同振幅的信号。

方波信号的函数为:, 如图1所示。

利用公式 ( 2) , 在MATLAB中编程, 通过求解Fourier系数a ( k) 和b ( k) , 实现离散傅里叶变换, 得到的振幅频谱图如图2所示。

从图2中可以看到: 方波信号的离散傅里叶变换后, 得到的系数几乎是平均分布在整个频率轴上, 也就是说方波信号实际上是包含了多种频率的正弦信号的叠加。其中频率为0. 25Hz的正弦信号振幅最大, 几乎是下一个频率为0. 5Hz的正弦信号振幅的3倍。而且, 随着频率的成倍提高, 振幅也不断地衰减。

从方波信号的频谱图中可以得出这样的结论: 方波信号实际上是由多个频率的正弦信号叠加而成的, 取的频率分量越多, 即k值越多, 与原始的方波信号越接近。但是学生要理解这个结论并不容易, 不能得到直观的印象。

修改调整MATLAB中的谐波次数, 分别把参数k设为1、3、9和33, 对比得到的结果。

当把参数k设置为1时, 只是把方波分解为一个正弦振动, 其波形与方波信号相差很大, 如图3所示。当把参数k设置为3时, 在原来一个正弦振动的基础上, 增加了一次谐波, 得到的振动波形比较接近方波信号, 如图4所示。

当把参数k设置为9时, 把方波分解为五个正弦振动, 其波形与方波信号接近了不少, 如图5所示。当把参数k设置为33时, 谐波振动增加到17个, 波峰和波谷基本被拉平, 得到的振动波形已经近似于方波信号了, 如图6所示。

通过这个例子, 学生对方波信号实际上是由无数个频率不同、振幅不同的正弦信号叠加而成的含义理解得更加透彻, 对离散傅里叶变换的概念和算法掌握得更加牢固, 而且通过编写MATLAB程序, 锻炼了分析和解决实际问题的能力。

四、结束语

通过具体实例说明了在数字信号处理课程中周期方波信号的频率分解的MATLAB辅助教学过程, 通过图形展示了不同频率不同振幅的正弦信号叠加为周期方波信号的过程, 取得了非常好的教学展示效果。MATLAB非常适合作为教学辅助工具, 它不仅提供强大的编程和运算功能, 而且能够以丰富的图形系统将结果形象地展示出来。因此有必要将数字信号处理的理论教学和MATLAB的应用结合起来, 改革教学方法, 使学生在学习基本理论的同时学会应用MATLAB进行信号与系统的分析, 把枯燥的公式化为生动的图形, 并且通过调整参数来优化显示结果, 这将极大地提高学生学习钻研的兴趣, 加深对基本概念和基本理论的理解, 增进学习的深度和广度。

摘要：周期方波信号的频率分解是数字信号处理课程中合成频率信号分解部分的重要案例, 有利于学生加深对周期信号和离散傅里叶变换的理解和运用。在传统课堂教学中, 板书讲解数学公式存在形式单一、内容枯燥、学生兴趣不高、不易理解的特点。提出了以MATLAB为手段辅助数字信号处理课程教学方法, 提高教学效率和学生兴趣, 通过对周期方波信号频率分解问题的具体分析, 运用MATLAB强大的编程功能实现信号的时间域和频率域仿真, 把离散傅里叶变换得到的结果绘制成图形图像, 调整参数逼近信号波形, 形象生动地展现了MATLAB在数字信号处理课程辅助教学中的重要作用, 使数字信号处理课程真正成为数学基本理论与实际应用相联系的纽带与桥梁。

关键词：MATLAB,数字信号处理,辅助教学

参考文献

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

[2]曾伟梁, 刘颖.论MATLAB三维仿真在数学模型课程教学实践中的作用[J].哈尔滨师范大学自然科学报, 2006, (3) .

非均匀时域序列中多周期信号探测 第3篇

为了解决普通傅立叶变换在含噪非均匀时域序列应用中的局限性, 迄今为止, 人们做了一些相关研究。1999年余建航, 张曾锠等人采用时域平均处理的方法从非均匀的时域序列中提取出弱周期信号[3], 时域平均处理是从混有噪声的复杂周期信号中提取简谐分量的有效方法之一, 但是研究结果发现周期测量误差对简谐分量的幅值和相位有较大影响。2001年贾焕玉等人用周期折叠法从非均匀时域序列中寻找周期信号, 并给出其相位和振幅, 但是该方法只能给出序列中所存在周期信号的周期、振幅、相位, 不能对周期信号的显著性进行准确估计, 对于多目标信号情形面临信号分离的困难。由lomb[4]创立, 经Scargle[5]发展的lomb-Scargle Periodogram算法目前被广泛应用于天文、地球物理和生物医学等学科领域的非均匀实验观测数据的频谱分析, 但是该方法存在探测效率低, 精度差的缺点[6,7]。本文采用加权叠代理论对该法进行优化处理, 得到一种更适合与非均匀时域信号中多周期信号探测的有效方法Weithted lomb-Scargle Periodogram, 最后用模拟的方法检验了该方法与其他方法的优越性。

1 Weighted Lomb-Scargle Periodogram

对于某给定的时域序列X (t i) , i=, 1, 23L, n, 传统Lomb-Scargle Periodogram方法是以最小二乘法为理论基础, 定义其拟合方程:

这里X (t i) 是离散实验数据;ti是离散实验数据的取样时间, n为实验数据统计量。离散的频率变量f∈[fmin, fmax], 其中fmin为序列的极限频率, fmax必须不大于序列的尼奎斯特频率fN, 频率取样部长∆f=fmin。

下面定义序列中各频率分量的幅度为如下变量:

其中K=fmax/∆f。采用最小二乘拟合估计:

则各频率分量的估计值:

其中

则Lomb-Scargle Periodogram功率谱定义为:

Lomb-Scargle Periodogram方法在处理非均匀时域序列时, 会在真实信号两侧产生虚假谱峰。另外, 多信号间的相互调制、背景噪声和取样起始时间的平移对谱结构有重大影响。

首先, 引入时间平移不变量τ (f) , 建立如下新的理论模型:

其中τ (f) 为时间平移不变常量定义如下:

考虑到其他频率分量和噪声对该频率信号的影响, 定义加权因子:

其中σi为各取样时刻的误差。

改进后的Weighted Lomb-Scargle Periodogram谱分析法可以有效地减小各频率信号间的相互调制和噪声对功率谱结果的影响, 大大地提高信号参数的估计结果。引入时间平移不变量后消除了取样起始时间变化对谱分析结果的影响。

2模拟分析与讨论

为了分析迭代加权处理后的Weighted Lomb-Scargle Periodogram算法在处理非均匀时域序列中多目标信号探测的优越性, 我们采用Mont Carlo模拟的方法, 产生具有三个周期调制信号 (频率分别为0.1、0.4和0.41Hz, 信噪比依次为6.0、12、14分贝) 的100个含噪非均匀序列, 然后分别使用Lomb-Scargle Periodogram和Weighted Lomb-Scargle Periodogram算法对其进行频谱分析, 结果如下图所示, 圆圈为模拟信号的信噪比, 散点为其傅立叶变换结果。

从图1, 我们可以看到上述两种方法都能准确地探测到模拟序列中的三个周期信号, 其频率估计值与模拟信号频率一样, 但是在信号强度的估计上, 加权处理后的Lomb-Scargle Periodogram的处理结果与模拟序列的真实结果几乎一致, 而Lomb-Scargle Periodogram算法的信噪比估计值与真实值存在较大的差异。图a显示, 模拟信号的信噪比越低Lomb-Scargle Periodogram算法的探测精度也越差。但是Weighted Lomb-Scargle Periodogram算法的探测精度与信号的信噪比几乎没用影响。

Weighted Lomb-Scargle Periodogram算法不仅可以较准确地从非均匀序列中探测到各周期信号, 该方法还能给出各周期信号在一个时间后期的相位分布, 然后再采用曲线拟合, 可以估计各周期信号的振幅和初相位。图2给出了模拟序列中的一个频率为1.61Hz的周期信号的相位分布。

模拟序列中, 人为加入的正弦周期信号的幅度为0.8, 除相位为1.83, 均值为0。从上图可以看到, Weighted Lomb-Scargle Periodogram算法得到的估计参数与真实参数非常接近, 足以说明其探测精度很高。其正弦曲线拟合结果也显示, 拟合曲线全部落在95%的置信区间。

3结语

MontCarlo模拟序列的频谱分析和参数估计结果表明, 改进后的Weighted Lomb-Scargle Periodogram傅立叶变换可以有效地减小时域序列的非均匀性在傅立叶变换中光谱泄露和各目标信号间的相互调制的影响, 提高信号探测的效率和参数估计的准确性。引入时间平移不变量后消除了序列的取样起始时间变化对谱分析和相位分布结果的影响。因此该方法是非均匀序列中多目标信号探测和参数估计的有效方法之一, 它解决了天文、地球物理和生物医学等学科领域中长期普遍存在的非均匀实验序列中多目标信号探测的问题。

摘要：天文、地球物理、生物医学等学科研究中, 面临着诸多的时间取样不均匀实验观测序列, 其中携带有丰富的有规律的周期信号, 从非均匀的实验观测序列中提取出周期信号是非常必要的。本文采用迭代加权理论对广泛应用于非均匀序列的周期信号探测的Lomb-Scargle Periodogram算法加以优化处理, 得到一种更优化的多目标信号探测方法。最后, 用MontCarlo模拟法对该方法与其前身的优越性予以比较, 并给出模拟序列中周期信号的相位分布。

关键词：非均匀时域序列,周期信号,探测

参考文献

[1]S.Bourguignon, H.Carfantan, and T.B hm.Astronomy&Astrophysics, 2007, 46:379～387.

[2]Yi Liu, JunHui Fan, HongGuang Wang.Chinese Journal of Astronomy and Astrophysics, 2006, 6:353～356.

[3]余建航, 张曾锠.周期测量误差对周期信号时域平均处理结果的影响[J].振动工程学报, 1999, 12 (2) :338.

[4]Lomb, N.R., Least-squares frequency analysis of unequallyspaced data.Atrophy space Sci, 1976, 39:447～462.

[5]J.D.Scargle.Studies astronomicaltime series analysis.Ⅱ.TheAstrophysical Journal, 1982, 263:835～853.

[6]简念川, 王广利, 李金岭, 等.基于Lomb算法的非等间距VLBI观测资料频谱分析研究[J].天文学报, 2006, 47 (3) :336～347.

[7]K.Hocke and N.Kampfer.Gap filling and noise reduction ofunevenly sampled data by means of the Lomb-Scargleperiodogram.Atmospheric Chemistry and Physics Discussions, 2008, 8:4063～4623.

[8]M.Schimmel, Emphasizing difficulties in the detection of rhythmswith Lomb-Scargle periodogram.Biol.Rhythm Res., 2001, 32:341～345.

信号周期 第4篇

关键词：车速波动信号,断纸预测,短时功率谱,周期统计平均值,有色噪声

1 引 言

造纸车速波动信号是一类重要的非平稳周期随机信号, 它可以用一个周期和幅度都呈正态分布的周期函数来表示。非平稳周期信号的周期或频率是一个随机变量, 而周期或频率的统计平均值 (以下简称平均周期和平均频率) 在工程实践中十分重要, 造纸机的车速波动信号平均周期$\overline{{\rm T}}$可以作为故障诊断的主要依据, 能够判断断纸的发生[1]。

由于非平稳周期信号的非平稳性, 其周期的统计平均值在不同时刻可能不同, 计算$\overline{{\rm T}}$并不容易。考虑到非平稳周期信号的平均周期可能具有时变性, 本文采用短时功率谱分析这样以时间和频率为变量的二维车速波动信号, 以获得当前时刻的时间或频率的统计平均值[2]。同时本文从时频分辨率、时频聚集性、抗干扰能力和动态性能几个方面对短时功率谱方法进行研究, 得到明确的结论。

2 非平稳周期信号及其周期的统计均值

2.1 非平稳周期信号

有一类随机信号可以表示成如下形式:

s (t) =As (t+nT) +v (t) (1)

式 (1) 中周期T是一个随机变量, 它围绕着平均周期$\overline{{\rm T}}$呈正态分布;幅值A也是随机变量, 它围绕均值EA (值为1) 呈正态或均匀分布;v (t) 是个0均值噪声。对公式 (1) 取数学期望, 由于A与T互不相关, 因此有:

$Es(t)=E[As(t+n{\rm T})+v(t)]=Es(t+n\overline{{\rm T}})(2)$

因此s (t) 是一个统计平均值具有周期性的非平稳信号, 我们称之为非平稳周期信号[3]。

非平稳周期信号在实践中广泛存在。例如, 对具有一定速度的气体-固体颗粒的二相混合物进行超声波测量, 获得的超声波回波信号就是一个非平稳周期信号。造纸机车速信号中的波动成分, 也是非平稳周期信号。

2.2 非平稳周期信号的周期 (频率) 统计平均值

对于非平稳周期信号, 其周期T的统计平均值$\overline{{\rm T}}$ (=ET) 十分重要。如上所述的造纸机的车速波动信号s (t) 的$\overline{{\rm T}}$能够判断断纸的发生[4]。

由于非平稳周期信号的非平稳性, 其周期的统计平均值在不同时刻可能不同, 计算$\overline{{\rm T}}$并不容易。笔者在研究中发现, 对s (t) 进行时频变换, 得到二维谱S (t, f) 的峰值所对应的频率fmax就是t时刻平均频率$\overline{f}$的良好估计值, 因此, 对非平稳随机信号的时频分析是获得平均频率$\overline{f}$或平均周期$\overline{{\rm T}}$的有效途径。

3 车速波动信号s (t) 的时频分析

以车速波动信号s (t) 为研究对象, 图1是部分s (t) 信号。显然从时域中获得其$\overline{f}$或$\overline{{\rm T}}$都非常困难。

与短时傅立叶变换相比, 短时功率谱结合了短时傅立叶变换和功率谱估计的特点, 先对非平稳信号进行加窗平稳化, 再进行功率谱估计, 并通过时间窗的移动, 对每个时刻都进行这种运算:

Px (t, f;$h)=\left|\frac{1}{\sqrt{2\text{π}}}{\displaystyle {\int }_{-\infty }^{\infty }s}(u)h{}^{*}(u-t)e{}^{-\text{j}2\text{π}fu}\text{d}u\right|{}^2(3)$

这里h (t) 是以t=0为中心的分析窗, 它可以有效抑制时刻u=t的邻域外的信号, 所以短时功率谱获得的谱图是信号在时刻t的邻域内的局部频谱。短时功率谱忽略了相位信息, 具有很强的统计特性。虽然和短时傅立叶变换一样, 短时功率谱也会对频谱造成泄漏, 但只要它不引起峰值频率的位置偏移和峰值的大幅度下降, 就不会对非平稳信号平均频率的计算造成影响。图2是对车速波动信号s (t) 的短时功率谱图[5]。

由图2可见, 图2 (a) 的短时功率谱图中峰值频率的谱线 (白色箭头指出) 清晰, 反映出其良好的时频分辨率, 由此可以获得每一时刻t精确的$\overline{f}$值。但是在t=9～11 s之间谱线出现断裂, 在t=13.5 s附近的谱线出现分叉和偏移。在图2 (b) 中经过惯性滤波处理后, 断裂、分叉的谱线被很好地修复, 但是谱线变宽, 说明分辨率略有下降, 同时惯性滤波不可避免地会对运算过程的动态性能略有影响。由图2 (b) 测出的频率约57.8～58.6 kHz, 与车速波动信号s (t) 在这个时段实际频率58.4 kHz的相对误差仅为1%, 精度较高。

4 噪声对非平稳随机信号时频分析的影响

非平稳周期信号s (t) 由于受到造纸车间环境的影响不可避免地会引入噪声, 式 (1) 中v (t) 可能是白噪声或有色噪声[6], 下面分两种情况分别进行讨论。

4.1 白色噪声

在非平稳周期信号s (t) 中加入幅度与它相同的高斯白噪声, 成为白噪声信号s′ (t) , 分析短时功率谱的抗干扰能力。图3分别是白噪声信号s′ (t) 的短时功率谱和经过较强惯性滤波的短时功率谱。

由图3 (a) 可知短时功率谱受到白噪声的影响, 失去对信号s′ (t) 的分析能力。图3 (b) 是经过较强的惯性滤波的短时功率谱图, 它的时频分辨率好, 谱线清晰, 能够将信号与噪声完全分离, 并能够对噪声成分略有抑制作用, 其抑制作用能够保证谱图中噪声的幅度远低于非平稳周期信号峰值的幅度。

总之, 白噪声对短时功率谱有较大的影响, 其中加强惯性滤波后, 短时功率谱可以克服白噪声的影响。

4.2 有色噪声

在非平稳周期信号s (t) 中加入幅度相当于0.8|s (t) |max的有色噪声, 成为有色噪声非平稳周期信号, 分析短时功率谱时频方法的抗干扰能力。

图4是有色噪声非平稳周期信号s″ (t) 经过强惯性滤波的短时功率谱。由图4可见, 经过惯性滤波的短时功率谱图, 时频分辨率好, 谱线 (白色箭头所指) 较为清晰, 谱图中信号与噪声完全分离, 但对噪声谱 (黑色箭头所指) 的抑制作用不足, 噪声部分的峰值甚至强于谱线, 为此还需要对噪声谱进一步处理。

总之, 有色噪声对具有惯性滤波短时功率谱有较强的影响, 短时功率谱可以通过加强滤波强度和采用进一步的噪声处理来减弱影响的强度, 但需要以牺牲时频分析过程的动态性能为代价。

5 结 论

非平稳周期信号的平均频率或周期具有十分重要的物理意义。在为获取平均频率或周期进行的时频分析方法中, 短时功率谱估计的时频分析能力和抗干扰能力都较差, 但具有惯性滤波作用的短时功率谱, 可以根据噪声的强度调整滤波程度, 在牺牲部分动态特性后, 能够获得比较好的时频分辨率和抗干扰能力。因此, 具有惯性滤波的短时功率谱时频分析能够很好地应用于要求高精度和快速动态性能、现场干扰强的断纸预测控制中。

参考文献

[1]周强, 韩九强.基于功率谱估计的造纸生产线断纸预测[J].化工自动化及仪表, 2007, 34 (2) :13-16.

[2]李羽中, 程虹霞, 林建辉.一种有效实现弱平稳随机信号完整功率谱估计方法[J].中国测试技术, 2006, 32 (1) :114-116.

[3]雷良育, 张卫平, 王江涛.非平稳随机路面模拟与分析[J].传感技术学报, 2007, 20 (4) :877-880.

[4]贾明平.周期平稳故障信号分析[J].东南大学学报, 2006, 36 (3) :346-349.

[5]胡昌华, 周涛, 夏启兵, 等.基于MATLAB的系统分析与设计——时频分析[M].西安:西安电子科技大学出版社, 2002:120-140.

信号周期 第5篇

目前, 针对DSSS 信号的扩频波形和信息序列的估计都是分开研究的。一般的做法是首先估计扩频波形, 然后通过解扩获取信息序列, 或者通过侦察信号的互相关实现盲解扩。其中, 扩频序列的估计方法主要有Massey 算法[1]、三阶相关算法[2,3]、基于特征值分析的盲估计算法[4,5,6,7,8,9]等。Massey 算法可以估计线性反馈移位寄存器 (LFSR) 序列 (如m序列) , 但是该方法不能适应低信噪比, 而且无法估计非线性序列;三阶相关算法利用m序列的三阶相关特性可以分析估计m 序列, 但是目前的研究也仅限于m 序列;基于特征值分析[6]的盲估计算法不受扩频序列类型的限制, 信噪比适应能力也比较强, 但是这种方法存在估计两段扩频序列存在正负模糊以及前后顺序不能确定等缺点。

本文通过对信号子空间的分析, 研究扩频波形与信息序列同时盲估计的有效方法, 用2倍PN周期的时间窗对接收信号进行分割, 并组成矩阵, 然后对矩阵进行SVD分解, 最后利用分解得到的左右奇异值矢量进行PN序列和信息序列盲估计。该方法避免文献[7]形成信号子空间相关矩阵运算, 降低了计算量;同时, 利用一个向量进行扩频序列估计, 避免传统方法[4,6]利用两个向量估计连接而产生部分序列反相问题, 进一步提高了盲估计的正确性。

1 信号模型

假设信号是扩频周期已知的短码扩频序列[5,6], 且DS-SS信号模型为:

$s(t)={\displaystyle \sum _{k=-\infty }^{+\infty }a}{}_{k}h(t-k{\rm T}{}_{\text{s}})(1)$

式中:Ts为符号周期;{ak=±1, k∈Ζ}为等概率随机分布的信息符号序列;h (t) 为一个完整周期的扩频序列基带信号与传输链路所有滤波器的卷积, 满足:

$h(t)={\displaystyle \sum _{i=0}^{{\rm N}-1}c}{}_{i}p(t-i{\rm T}{}_{\text{c}})(2)$

式中:{ci=±1, i=0, 1, 2, …, N-1}为扩频码序列;p (t) 为发射机滤波器与信道冲击响应、接收机滤波器的卷积;Tc为码片时宽;N为扩频增益。在无干扰情况下, 接收机接收信号为:

$y(t)=s(t)+n(t)(3)$

式中:n (t) 为功率谱密度为σ${}_n^2$的高斯白噪声。

不失一般性, 本文做如下假设:

(1) 信息序列为均匀分布且互不相关;噪声为零均值高斯白噪声, 且与信号不相关。

(2) 失步时间t0可按照文献[4,10]的方法计算求得。

2 扩频波形与信息序列盲估计

2.1 信号矩阵组成

传统基于特征值分解的求解短码扩频序列方法[4,6]中信号分段方法如图1所示。由图可知, 该分段仅包含一个周期的PN序列, 当且仅当t0=0时, 所估计的序列才包含一个整周期的PN序列估计, 而在实际的取窗过程中, 很难保证t0=0。当t0≠0时, 传统的分段方法包含2段PN序列估计;只有当正确组合时才得到PN序列估计。组合过程中存在相位模糊问题, 同时, 传统的方法不能盲估计信息序列。

本文的信号分段方法如图2所示, 无论t0是否为0, 该分段总能确保包含一个整周期的PN序列估计, 因此可以避免组合过程中相位模糊问题, 并能同时盲估计信息序列。 令信号通过持续时间为2Ts, 且不重叠观测窗[7], 所得观测样本矢量为xk=[xk-1, xk-2, …, xk-2N]T, 假设采样周期等于码片宽度, 则xk为2N×1维矢量, 表示为:

$\mathit{x}{}_k=\mathit{s}{}_k+\mathit{n}{}_k(4)$

假设0≤t0<Ts, 则有:

$\mathit{s}{}_k=a{}_k\mathit{h}{}_0+a{}_{k+1}\mathit{h}+a{}_{k+2}\mathit{h}{}_1(5)$

式中:

$\begin{array}{l}\mathit{h}{}_0=\{\begin{array}{cc}\mathit{h}(t)‚& k{\rm T}{}_{\text{s}}+t{}_0\le t<(k+1){\rm T}{}_{\text{s}}\\ 0‚& (k+1){\rm T}{}_{\text{s}}\le t<(k+2){\rm T}{}_{\text{s}}+t{}_0\end{array}\\ \mathit{h}=\{\begin{array}{cc}0‚& k{\rm T}{}_{\text{s}}+t{}_0\le t<(k+1){\rm T}{}_{\text{s}}\\ \mathit{h}(t)‚& (k+1){\rm T}{}_{\text{s}}\le t<(k+2){\rm T}{}_{\text{s}}\\ 0‚& (k+2){\rm T}{}_{\text{s}}\le t<(k+2){\rm T}{}_{\text{s}}+t{}_0\end{array}(6)\\ \mathit{h}{}_1=\{\begin{array}{cc}\mathit{h}(t),& (k+2){\rm T}{}_{\text{s}}\le t<(k+2){\rm T}{}_{\text{s}}+t{}_0\\ 0,& k{\rm T}{}_{\text{s}}+t{}_0<t<(k+2){\rm T}{}_{\text{s}}\end{array}\end{array}$

另假设窗个数为M, 则有:

$\begin{array}{l}\mathit{X}=\frac{1}{\sqrt{{\rm M}}}\left[\begin{array}{c}x{}_1\\ x{}_2\\ ⋮\\ x{}_{{\rm M}}\end{array}\right]=\frac{1}{\sqrt{{\rm M}}}\\ \left[\begin{array}{c}a{}_1\\ a{}_3\\ ⋮\\ a{}_{2{\rm M}-1}\end{array}\right]\mathit{h}{}_0+\left[\begin{array}{c}a{}_2\\ a{}_4\\ ⋮\\ a{}_{2{\rm M}}\end{array}\right]\mathit{h}+\left[\begin{array}{c}a{}_3\\ a{}_5\\ ⋮\\ a{}_{2{\rm M}+1}\end{array}\right]\mathit{h}{}_1+\left[\begin{array}{c}n{}_1\\ n{}_2\\ ⋮\\ n{}_{{\rm M}}\end{array}\right]=\\ \frac{1}{\sqrt{{\rm M}}}(\mathit{a}{}_0\mathit{h}{}_0+\mathit{a}{}_1\mathit{h}+\mathit{a}{}_2\mathit{h}{}_1+\mathit{n})(7)\end{array}$

式中:a0= (a0, a1, …, aM-1) T, a1= (a1, a2, …, aM) T, a2= (a2, a3, …, aM+1) T为观测信息序列矢量。设每个数据窗内的采样点数为$2{\rm N}=\frac{2{\rm T}{}_{\text{s}}}{{\rm T}{}_{\text{c}}}$, 则X为M×2N维矩阵。

2.2 信息序列和PN序列估计分析

根据奇异值分解定理, 可得X的奇异值分解:

$\mathit{X}=\mathit{U}\mathit{\Delta }\mathit{V}{}^{\text{Η}}(8)$

式中:U和V为酉矩阵;“H”表示矩阵的共轭转置;Δ为准对角矩阵:

$\mathit{\Delta }=\text{d}\text{i}\text{a}\text{g}(a{}_1,a{}_2\cdots ,a{}_{{\rm M}})(9)$

式中:ai (i=1, 2, …, M) 为X的非零奇异值, 且a1≥a2≥…≥aM>0。

2.2.1 信息序列估计证明

下面, 证明矩阵Χ的SVD分解可得到短码扩频信号的信息序列和扩频序列的估计值。

证明:

首先, 定义协方差矩阵:

$\mathit{R}{}_1=E\{\mathit{X}\mathit{X}{}^{\text{Η}}\}(10)$

式中:E{·}为期望。将式 (8) 、式 (9) 代入式 (10) 可得R1的特征值分解:

$\begin{array}{l}\mathit{R}{}_1=\mathit{U}\mathit{\Delta }\mathit{V}{}^{\text{Η}}\mathit{V}\mathit{\Delta }{}^{\text{Η}}\mathit{U}{}^{\text{Η}}=\mathit{U}\mathit{\Delta }\mathit{\Delta }{}^{\text{Η}}\mathit{U}{}^{\text{Η}}\\ =\mathit{U}{\displaystyle \sum \mathit{U}}{}^{\text{Η}}(11)\end{array}$

非零特征值λi满足:

$\lambda {}_i=a{}_i^2,i=1,2,\cdots ,{\rm M}(12)$

因信息序列ak为±1的均匀分布随机序列, 则aiaj (i≠j) 也是±1均匀分布的随机序列, 因此有:

$E\left\{\frac{1}{{\rm M}}\mathit{a}{}_i^{\text{Η}}\mathit{a}{}_j{}^{}\right\}=E\left\{\frac{1}{{\rm M}}{\displaystyle \sum _{\begin{array}{l}i,j=1；\\ i\ne j\end{array}}^{{\rm M}}a}{}_{i}a{}_j\right\}=0(13)$

又因窗个数为M, 则:

$\parallel \mathit{a}{}_0\parallel {}^2=\parallel \mathit{a}{}_1\parallel {}^2=\parallel \mathit{a}{}_2\parallel {}^2={\displaystyle \sum _{i=1}^{m}a}{}_i^2={\rm M}(14)$

由式 (6) 得双周期内的各部分扩频序列积分为:

$\begin{array}{l}{\displaystyle {\int }_0^{2{\rm T}{}_{\text{s}}}\mathit{h}}{}_0(t)\mathit{h}{}_1(t)\text{d}t=0;\\ {\displaystyle {\int }_0^{2{\rm T}{}_{\text{s}}}\mathit{h}}{}_0(t)\mathit{h}{}_{}(t)\text{d}t=0\\ {\displaystyle {\int }_0^{2{\rm T}{}_{\text{s}}}\mathit{h}}{}_{}(t)\mathit{h}{}_1(t)\text{d}t=0\\ {\displaystyle {\int }_0^{2{\rm T}{}_{\text{s}}}\mathit{h}}{}_0(t)\mathit{h}{}_0(t)\text{d}t={\displaystyle {\int }_{{\rm T}{}_{\text{s}}+t{}_0}^{2{\rm T}{}_{\text{s}}}\mathit{h}}(t)\mathit{h}(t)\text{d}t\\ ={\displaystyle {\int }_{{\rm T}{}_{\text{s}}+t{}_0}^{2{\rm T}{}_{\text{s}}}\mathit{h}}{}_0(t)\mathit{h}{}_0(t)\text{d}t={\rm T}{}_{\text{s}}-t{}_0\\ {\displaystyle {\int }_0^{2{\rm T}{}_{\text{s}}}\mathit{h}}{}_1(t)\mathit{h}{}_1(t)\text{d}t={\displaystyle {\int }_0^{t{}_0}\mathit{h}}(t)\mathit{h}(t)\text{d}t\\ ={\displaystyle {\int }_0^{t{}_0}\mathit{h}}{}_1(t)\mathit{h}{}_1(t)\text{d}t=t{}_0\\ {\displaystyle {\int }_0^{2{\rm T}{}_{\text{s}}}\mathit{h}}(t)\mathit{h}(t)\text{d}t={\displaystyle {\int }_{t{}_0}^{{\rm T}{}_{\text{s}}+t{}_0}\mathit{h}}(t)\mathit{h}(t)\text{d}t={\rm T}{}_{\text{s}}\end{array}(15)$

将式 (7) 代入式 (10) 得:

$\begin{array}{l}\mathit{R}{}_1=E\{\frac{1}{{\rm M}}\parallel \mathit{h}{}_0\parallel {}^2\mathit{a}{}_0\mathit{a}{}_0^{\text{Η}}+\frac{1}{{\rm M}}\parallel \mathit{h}\parallel {}^2\mathit{a}{}_1\mathit{a}{}_1^{\text{Η}}+\\ \frac{1}{{\rm M}}\parallel \mathit{h}{}_1\parallel {}^2\mathit{a}{}_2\mathit{a}{}_2^{\text{Η}}+\sigma {}_n^2\mathit{{\rm I}}\}\\ =\frac{1}{{\rm M}}\parallel \mathit{h}{}_0\parallel {}^2\mathit{a}{}_0(\mathit{a}{}_0^{\text{Η}})+\frac{1}{{\rm M}}\parallel \mathit{h}\parallel {}^2\mathit{a}{}_1(\mathit{a}{}_1^{\text{Η}})+\\ \frac{1}{{\rm M}}\parallel \mathit{h}{}_1\parallel {}^2\mathit{a}{}_2(\mathit{a}{}_2^{\text{Η}})+\sigma {}_n^2(\mathit{{\rm I}})(16)\end{array}$

用a0右乘式 (16) 并根据式 (14) , 式 (15) 可得:

$\begin{array}{l}\mathit{R}{}_1\mathit{a}{}_0=E\{\frac{1}{{\rm M}}\parallel \mathit{h}{}_0\parallel {}^2\mathit{a}{}_0\mathit{a}{}_0^{\text{Η}}+\frac{1}{{\rm M}}\parallel \mathit{h}\parallel {}^2\mathit{a}{}_1\mathit{a}{}_1^{\text{Η}}+\\ \frac{1}{{\rm M}}\parallel \mathit{h}{}_1\parallel {}^2\mathit{a}{}_2\mathit{a}{}_2^{\text{Η}}+\sigma {}_n^2\mathit{{\rm I}}\}\mathit{a}{}_0\\ =\frac{1}{{\rm M}}\parallel \mathit{h}{}_0\parallel {}^2\mathit{a}{}_0(\mathit{a}{}_0^{\text{Η}}\mathit{a}{}_0)+\frac{1}{{\rm M}}\parallel \mathit{h}\parallel {}^2\mathit{a}{}_1(\mathit{a}{}_1^{\text{Η}}\mathit{a}{}_0)+\\ \frac{1}{{\rm M}}\parallel \mathit{h}{}_1\parallel {}^2\mathit{a}{}_2(\mathit{a}{}_2^{\text{Η}}\mathit{a}{}_0)+\sigma {}_n^2(\mathit{{\rm I}}\mathit{a}{}_0)\\ =\left(\frac{1}{{\rm M}}\parallel \mathit{h}{}_0\parallel {}^2\parallel \mathit{a}{}_0\parallel {}^2+\sigma {}_n^2\right)\mathit{a}{}_0\\ =(\parallel \mathit{h}{}_0\parallel {}^2+\sigma {}_n^2)\mathit{a}{}_0(17)\end{array}$

可见, a0是R1的特征向量, 对应特征值为‖h0‖2+σ${}_n^2$。同理, a1, a2也是R1的特征向量, 对应特征值分别为‖h‖2+σ${}_n^2$和‖h1‖2+σ${}_n^2$;R1其他特征向量对应的非零特征值为σ${}_n^2$。

设u1′, u2′, u3′分别为R1前3个较大值对应的归一化特征向量, 则由上面的分析可知, 信息序列向量a1, a0, a2也是R1前3个最大特征值对应的特征向量, 它们仅存在一个常系数差别。根据式 (16) 以及‖u1′‖2 =‖u2′‖2 =‖u3′‖2 = 1, 可得a1, a0和a2的估计值为:

$\widehat{\mathit{a}}{}_1=\pm \sqrt{{\rm M}}\mathit{u}{}_{1^{\prime }}‚\widehat{\mathit{a}}=\pm \sqrt{{\rm M}}\mathit{u}{}_{2^{\prime }}‚\widehat{\mathit{a}}{}_2=\pm \sqrt{{\rm M}}\mathit{u}{}_{3^{\prime }}(18)$

式中:“±”说明$\widehat{\mathit{a}}$1, $\widehat{\mathit{a}}$0和$\widehat{\mathit{a}}$2存在正负模糊。由此得证矩阵R1的特征分解可包含信息序列估计向量。

2.2.2 扩频序列估计证明

定义如下矩阵:

$\mathit{R}{}_2=\mathit{E}\{\mathit{X}{}^{\text{Η}}\mathit{X}\}(19)$

将式 (7) 代入式 (19) 可得:

$\begin{array}{l}\mathit{R}{}_2=E\{\frac{1}{{\rm M}}(\mathit{a}{}_0\mathit{h}{}_0+\mathit{a}{}_1\mathit{h}+\mathit{a}{}_2\mathit{h}{}_1+\mathit{n}){}^{\text{Η}}(\mathit{a}{}_0\mathit{h}{}_0+\\ \mathit{a}{}_1\mathit{h}+\mathit{a}{}_2\mathit{h}{}_1+\mathit{n})\}\\ =\frac{1}{{\rm M}}\parallel \mathit{a}{}_0\parallel {}^2\mathit{h}{}_0^{\text{Η}}\mathit{h}{}_0+\frac{1}{{\rm M}}\parallel \mathit{a}{}_1\parallel {}^2\mathit{h}{}_{}^{\text{Η}}\mathit{h}{}_{}+\\ \frac{1}{{\rm M}}\parallel \mathit{a}{}_2\parallel {}^2\mathit{h}{}_1^{\text{Η}}\mathit{h}{}_1+\sigma {}_n^2\mathit{{\rm I}}\\ =(\mathit{h}{}_0^{\text{Η}}\mathit{h}{}_0+\mathit{h}{}_{}^{\text{Η}}\mathit{h}+\mathit{h}{}_1^{\text{Η}}\mathit{h}{}_1)+\sigma {}_n^2\mathit{{\rm I}}(20)\end{array}$

因为R2具有和R1相同的三个较大特征值, 设v1″, v2″和v3″分别为R2前3个大特征值对应的归一化特征向量。当t0<Ts/2且λ1=‖h‖2+σ${}_n^2$=Tsσ2s+σ${}_n^2$时, 满足:

$\begin{array}{l}\lambda {}_1\mathit{v}{}_{1^{″}}=\mathit{R}{}_2\mathit{v}{}_{1^{″}}\\ =((\mathit{h}{}_0^{\text{Η}}\mathit{h}{}_0+\mathit{h}{}^{\text{Η}}\mathit{h}+\mathit{h}{}_1^{\text{Η}}\mathit{h}{}_1)+\sigma {}_n^2{\rm I})\mathit{v}{}_{1^{″}}\\ =({\rm T}{}_{\text{s}}\sigma {}_{\text{s}}^2+\sigma {}_n^2)\mathit{v}{}_{1^{″}}(21)\end{array}$

结合式 (20) 得:

$\mathit{h}=k{}_1\mathit{v}{}^1˝(22)$

同理可得:

$\mathit{h}{}_0=k{}_2\mathit{v}{}_2˝,\mathit{h}{}_0=k{}_3\mathit{v}{}_3˝(23)$

式中:k1, k2和k3为比例系数。由此得证, 矩阵R2的特征分解可包含扩频序列估计向量。

2.3 信息和扩频序列联合估计结论

对信息序列的估计由上述分析可知, 无论失步时间t0是否大于Ts/2, 对应于最大奇异值的左奇异矩阵的归一化向量总为a1= (a1, a2, …, aM) T估计。但是, ‖h‖2+σ${}_n^2$和‖h1‖2+σ${}_n^2$的大小取决于t0, 当t0≤Ts/2时, 则对应次最大值的左奇异矩阵的归一化向量为$\widehat{\mathit{a}}{}_0=(a{}_0‚a{}_1‚a{}_2\cdots ‚a{}_{{\rm M}-1}){}^{\text{Τ}}$估计, 对应于第三大值的左奇异矩阵的归一化向量为$\widehat{\mathit{a}}{}_2=(a{}_2,a{}_5,\cdots ,a{}_{2{\rm M}+1})$估计;反之, 则对应次最大值的左奇异矩阵的归一化向量为$\widehat{\mathit{a}}{}_2=(a{}_2,a{}_5,\cdots ,a{}_{2{\rm M}+1})$, 对应于第三大值的左奇异矩阵的归一化向量为$\widehat{\mathit{a}}{}_0=(a{}_0‚a{}_1‚a{}_2‚\cdots ‚a{}_{{\rm M}-1}){}^{\text{Τ}}$估计。

对于扩频序列的估计, 由上述分析可知, 矩阵X分解后, 对应于最大奇异值的右奇异矩阵的向量包含整个扩频序列h的估计值。同时, 由式 (6) 定义可知, 当失步时间t0已知时, h估计向量去掉尾部t0长度后前推Ts为后得到的向量即为扩频序列估计。

同理也可利用右奇异矩阵第3列向量的后t0位 (h1) 与第2列向量的前Ts-t0位 (h0) 组合构成扩频序列估计, 仿真如图3 (a) 所示。但是, 由上述分析可知, 第2, 3列向量之间存在相位模糊, 所以在连接过程中存在反相的问题, 仿真见图3 (b) 所示, 为了便于阐述问题, 仿真时, 假设无噪声情况。因此, 利用第1列向量前推法获得扩频序列估计, 在同一向量空间不存在反相问题。

3 仿真分析

为了验证本文提出的算法, 设置了2个仿真试验。第一个实验通过仿真将本文算法的PN序列和信息序列估计性能进行了验证;第二个试验通过和其他方法进行仿真比较, 验证本文算法的优越性。仿真中直扩信号采用63位的m序列进行扩频, 码片速率10 MHz, 符号速率为157.5 kHz, 信号为BPSK调制, 信噪比为-7 dB, 仿真信号长度为100个扩频周期, 即50个信息码 (窗) 。采样频率为10 MHz, 每个扩频周期采样点数为63。接收端滤波器忽略不计。当t0=20, 矩阵分解后, 右特征矩阵第一列特征向量的一部分即为所求得扩频序列, 仿真分析见图4 (a) 。图中右侧较小码值的个数即为失步时间的码片数, 由此前推N为即得的扩频序列估计。同样, 为进一步验证, 图4 (b) 为无噪声时的比较示意图, 从图中可知, 该方法能够较好得到扩频序列估计。对应的信息序列估计如图5所示。图5 (a) ～ (c) 分别为信息序列$\widehat{\mathit{a}}$0、$\widehat{\mathit{a}}$1和$\widehat{\mathit{a}}$2估计值与真实值的比较示意图, 在t0<Ts/2时, 图5 (a) 为a1= (a2, a4, …, a2M) T估计示意图;图5 (b) 为a0= (a1, a3, a5, …, a2M-1) T估计示意图;图5 (c) 为$\widehat{\mathit{a}}$0一个数据符号延迟的估计, 则$\widehat{\mathit{a}}$1和$\widehat{\mathit{a}}$0的组合构成正确的信息序列。由图可知, 信息序列的估计值的极性误码都为0, 同时, 满足$\widehat{\mathit{a}}$2为$\widehat{\mathit{a}}$0的一个信息码延迟, 从而得证理论推导的正确性。

图3的仿真结果分别验证了第2.3节的分析结果。由图可知, 当利用文献[4,6]等所描述的算法, 链接PN序列估计的两部分时, 存在相位模糊现象, 如图5 (b) 所示。图6为采用本文算法以及对应文献[4,6]和互相关解扩算法时进行PN和信息序列估计。对于每个SNR值, 采用1 000次蒙特卡罗仿真, 且随机采用失步时间时, 估计矢量与真实矢量间的极性错误比值。由图可知, 本文算法在较低信噪比环境下仍能有效估计扩频波形和信息序列, 其中信息序列的估计性能接近互相关解扩的性能;扩频波形的估计性能优于文献[4,6]算法, 这是因为文献[4,6]算法在链接PN估计序列时需要去掉两段估计模糊。

4 结 语

本文提出了一种基于双PN周期的短码DS-SS信号扩频波形与信息序列联合估计的方法。该方法通过对奇异值分解所得左右奇异矩阵进行分析, 实现了扩频波形与信息序列的同时估计。该方法不需对接收信号进行相关矩阵积累, 降低了计算量。同时, 该方法利用同一向量空间进行扩频向量估计, 避免了传统的子空间估计方法所采用的部分长度估计再连接时存在的反相问题。

摘要：研究了DS-SS信号的扩频波形及信息序列盲估计问题, 提出一种基于双PN周期分解的扩频波形与信息序列盲估计算法, 该算法在低信噪比条件下可以同时完成扩频波形和信息序列的估计, 计算量小, 且可避免传统方法利用单PN周期分解时扩频序列连接的相位模糊等问题, 提高了盲估计的正确率。仿真结果验证了算法的有效性。

关键词：直接序列扩频,盲解扩,扩频序列,信息序列

参考文献

[1]ADAMS E R, GOUDA M, HILL P C J.Statistical tech-niques for blind detection&discrimination of M-sequencecodes in DSSS systems[C]//Proceedings of 1998 IEEE 5thInternational Symposium on Spread Spectrum Techniquesand Applications.[S.l.]:IEEE, 1998, 3:853-857.

[2]BOUDER C, AZOU S, BUREL G.A robust synchroniza-tion procedure for blind estimation of the symbol period andthe timing offset in spread spectrum transmissions[C]//Proceedings of IEEE International Symposium on SpreadSpectrum Techniques and Applications (ISSSTA) .Prague, Czech Republic:IEEE, 2002:238-241.

[3]BUREL G, BOUDER C.Blind estimation of the pseudo-random sequence of a direct sequence spread spectrum signal[C]//Proceedings of IEEE 21st Century Military Comm.Conf.Los Angeles, USA:IEEE, 2000:1-8.

[4]HILL P C J, RIDLEY M E.Blind estimation of direct-se-quence spread spectrum-sequence chip codes[C]//Proceed-ings of IEEE 6th Int.Symp.on Spread-Spectrum Technique&Application.New Jersey, USA:IEEE, 2000:305-309.

[5]YAO Ying-wei, VINCENT Poor H.Bind detection of syn-chronous CDMA in on-gaussian channels[J].IEEE Trans-actions on Signal Processing, 2004, 52, (1) :271-279.

[6]ZHANG T Q, MU A P, ZHANG C.Analyze the eigen-structure of DS-SS signals under narrow band interferences[J].Digital Signal Processing, 2006, 16:746-753.

[7]ZHANG T Q, MU A P.A modified eigen-structure analy-zer to lower SNR DS-SS signals under narrow band interfer-ences[J].Digital Signal Processing, 2007, 17:746-753.

[8]QUI P Y, HUANG Z T, JIANG W L, et al.Improvedblind-spreading sequence estimation algorithm for direct se-quence spread spectrum signals[J].IET Signal Process, 2008, l2:139-146.

[9]詹亚锋, 曹志刚, 马正新.DSSS信号的扩频序列估计[J].电子与信息学报, 2005, 27 (2) :169-173.

信号周期 第6篇

扩频码是直扩信号的重要参数之一, 在非合作条件下, 获取扩频码是解调直扩信号的前提。长码直扩信号不仅保密性能好, 还能抑制部分符号间的相互干扰, 在卫星通信中有着广泛应用, 因此长码直扩信号的参数估计有待深入研究。

目前针对直扩信号的扩频码估计, 文献[1-5]使用子空间方法;文献[6-7]利用Fast-ICA算法估计短码CDMA信号的扩频码;文献[8-9]利用Fast-ICA算法估计短码和非周期长码CDMA信号的扩频码。这些方法均不能直接应用于周期长码直扩信号的扩频码估计。文献[10]提出了一种周期长码直扩信号的虚拟多用户模型, 并利用二阶统计量实现了扩频码估计。文献[11]在虚拟多用户模型基础上, 利用改进二阶统计盲辨识 (SOBI) 算法来估计扩频码。但是文献[10]和文献[11]都是在信号已完成同步的前提下进行的。

在虚拟多用户模型的基础上, 提出了一种新的PLC直扩信号模型, 在新模型的基础上利用FastICA算法同时实现盲同步和扩频码估计。

1 信号模型

假设码片宽度为Tc, 扩频码周期为T0=PTc, 信息符号周期为Ts=MTc, 1周期扩频码调制G个信息符号, 那么有T0=GTs, 即P=GM, 其中P、G和M均为正整数。扩频码为{h (t) =±1, 0≤t

不失一般性, 做如下假设: (1) 信息序列{a[k], k=0, 1, …}服从独立同分布; (2) 噪声为零均值高斯白噪声, 与信号不相关; (3) 接收信号的扩频码周期数不小于G, 实际中容易实现。

假设接收信号的时延为td=d Tc, 其中d为整数。为确保一个时间窗内包含完整的扩频码, 采用的分段方法如图1所示。将接收信号以T0+T0/G为宽度, 重叠T0/G的时间窗进行分段。对接收信号以Tc为采样周期进行采样, 那么对应时间窗内的P+M维观测样本向量xk可以表示为:

1.1 同步信号模型

当时延td=0时, 即观测窗与信息符号对齐时, 则有:

上式等价于

式中,

当有 (L+1) 个窗口时, 则有

式中, 为信息序列组成的矩阵。

1.2 异步信号模型

当时延td≠0时, 即观测窗与信息符号没有对齐时, 则有:

上式等价于

式中,

当有 (L+1) 个窗口时, 则有

式中, 为信息序列组成的矩阵。

2 利用Fast-ICA估计时延和扩频码

由假设 (1) 可知, 式 (11) 和式 (19) 中的信息序列矩阵A中的每一列的各个元素是相互独立的, 因此式 (11) 和式 (19) 与盲源分离中的线性混合模型一致, 所以可使用盲源分离算法完成扩频码的估计。具体推导如下:

步骤1:白化观测信号

结合公式 (19) , 观测信号的相关矩阵可以表示为

其中, , σn2为噪声功率。由于hn G、hu (u=1, 2, …, G) 和hb1两两正交, 所以有:

同步信号下的特征值分解可以经过类似的推导得出。观测信号白化后的信号为:

对于异步信号, Z为 (G+2) × (L+1) 维;对于同步信号, Z为 (G+1) × (L+1) 维。

步骤2:利用Fast-ICA估计混合矩阵H

通过Fast-ICA可以利用白化后的信号Z求得分离矩阵W (具体步骤见文献[12]) , 使得:

式中, 是信息序列矩阵A的估计。当信号已同步时, W为 (G+1) × (G+1) 维;当信号未同步时, W为 (G+2) × (G+2) 维。在实际应用中, 由于大多数情况下事先不知道信号是否已经同步而且同步情况可以看作是异步情况的特例, 所以可将信号都当作异步信号处理。

假如没有噪声, 可以认为:

式中, Q为 (G+2) × (G+2) 的位置与符号模糊矩阵 (其每一行每一列均只有一个元素为+1或-1) , 因此矩阵Q满足:

由式 (11) 、式 (19) 、式 (28) 、式 (29) 、式 (30) 和式 (31) , 并且不考虑噪声影响可知:

当H=UsΛs1/2WHQ时, 式 (32) 成立。即使考虑噪声, 将上面的等号换成约等于号, 结论仍然成立。由于模糊矩阵Q是未知的, 所以实际上是利用:

来估计混合矩阵H。估计的混合矩阵中存在着位置和符号模糊, 可以根据扩频码起始点和PN码序列的性质来解模糊。

步骤3:解模糊, 估计时延和扩频序列

为向量cu的扩频码起始点, 则{τu, u=1, 2, …, G+2}中有G个扩频码起始点理论上将组成等差数列{d+1, d+1+M, …, d+1+ (G-1) M}。综合考虑向量cu内元素的幅度和扩频码起始点的等差数列特性, 可以从中选取G个向量, 其对应的扩频码起始点近似组成等差数列。将选取的向量按扩频码起始点由小到大的顺序排列, 记为{gi, i=1, 2, …, G};对应的起始点按从小到大的顺序排列, 记为{λi, i=1, 2, …, G}。时延d的估计可以表示为:

第i段扩频码的估计可以表示为:

式中, gi (λi:λi+M-1) 表示向量gi中从第λi个元素到第λi+M-1个元素组成的向量, sgn (·) 是符号函数。这样就利用扩频码起始点的等差特性和大小关系实现了位置解模糊。估计的扩频码可以表示为:

其中{bi=±1, i=1, 2, …, G}, 利用扩频码序列的相关特性可以确定参数bi, 实现符号解模糊[11]。

3 仿真结果

仿真中不考虑符号解模糊。1周期扩频码调制G个信息符号, 可以认为1周期扩频码被分为G段。借鉴文献[7]中的定义, 估计的扩频码的误极性率可以认为是G段扩频码的误极性率的平均值。

仿真1:时延估计和扩频码估计仿真

假设接收信号为采用BPSK调制的PLC直扩信号的基带信号, 扩频码长度为63, G=3, 信号的采样周期等于Tc, 时延td=11Tc, 信号长度为16 000个信息符号, 信噪比为-10 d B。仿真中得到的向量ci与其在无噪声情况下的理论值的对比如图3、图4、图5和图6所示。从图中可以看出:图3、图4和图5中向量的扩频码起始点构成等差数列。由图3可知:估计的时延为12-1=11, 与实际的时延相等, 这验证了前面理论推导的有效性。

仿真2:误极性率仿真

假设接收信号为采用BPSK调制的PLC直扩信号的基带信号, 扩频码长度为600。信号采样周期等于Tc, 时延为0~599的随机整数, 信号长度为4 000个扩频码周期。仿真实验结果是在每周期扩频码调制的信息符号数G分别取3、5和10以及不同信噪比下进行1 000次蒙特卡罗实验得到的。提出的基于Fast-ICA的周期长码直扩信号扩频码估计方法与文献[10]提出的算法估计出的扩频码的误极性率的比较如图7所示。sos-i表示当G取i时文献[10]提出的算法的误极性率曲线, ica-i表示当G取i时提出算法的误极性率曲线, i分别取3、5和10。从图7可以看出, 每周期扩频码调制的信息符号数G增加, 两种算法的误极性率都增加。这是因为当扩频码周期一定时, 每周期扩频码调制的信息符号数的增加意味着扩频增益的减小。从图7还可以看出, 以估计的扩频码的误极性率为0所需的最小信噪比作为比较标准, 文献[10]提出的算法比本文提出的方法优越1 d B或2 d B。这是因为文献[10]提出的算法是在信号已同步的基础上进行的, 而提出的方法不需要信号已同步。现有的周期长码直扩信号的扩频码估计算法 (如文献[10-11]中的算法) 大多是在信号已同步的基础上进行的, 提出的基于Fast-ICA的周期长码直扩信号扩频码估计方法则可以同时完成时延与扩频码的估计, 这是其优势所在。

4 结束语

在周期长码直扩信号的虚拟多用户模型的基础上建立了一种新的信号模型, 在新模型的基础上利用Fast-ICA算法进行盲同步和扩频码估计。传统的周期长码直扩信号的盲扩频码估计大多是在完成同步的基础上进行的, 本文中的方法可以同时实现盲同步与盲扩频码估计。提出的基于Fast-ICA的周期长码直扩信号盲扩频码估计方法在通信侦察中有广阔的应用前景。

参考文献

[1]BUREL G, BOUDER C.Blind Estimation of the Pseudorandom Sequence of A Direct Sequence Spread Spectrum Signal[J].IEEE MILCOM, 2000, 2 (10) :967-970.

[2]张红波, 吕明.基于子空间跟踪的扩频码盲估计算法[J].系统工程与电子技术, 2006, 28 (10) :1470-1532.

[3]吕明, 张红波, 唐斌.基于E-PASTd的盲扩频码序列估计算法[J].电子科技大学学报, 2007, 36 (5) :886-888.

[4]张天骐, 林孝康, 周正中.基于神经网络的低信噪比直扩信号扩频码的盲估计方法[J].电路与系统学报, 2007, 12 (2) :118-123.

[5]任啸天, 徐晖, 黄知涛.短码DS-SS信号扩频序列及信息序列联合盲估计方法[J].通信学报, 2012, 33 (4) :169-175.

[6]付卫红, 杨小牛, 刘乃安.基于盲源分离的CDMA多用户检测与伪码估计[J].电子学报, 2008, 36 (7) :1319-1323.

[7]任啸天, 徐晖, 黄知涛, 等.基于Fast-ICA同、异步系统短码CDMA信号扩频序列与信息序列盲估计[J].电子学报, 2011, 39 (12) :2726-2732.

[8]任啸天, 徐晖, 黄知涛, 等.基于Fast-ICA的CDMA信号扩频序列优化盲估计[J].电子学报, 2012, 40 (8) :1532-1538.

[9]陆凤波, 黄知涛, 姜文利.基于Fast-ICA的CDMA信号扩频序列盲估计及性能分析[J].通信学报, 2011, 32 (8) :136-142.

[10]DANESHMAND S, AGHAEINIA H, TOHIDIAN M, et al.Blind Estimation of Signal in Periodic Long-code DSSS Communications[C]∥SARNOFF’09 Proceedings of the32nd International Conference on Sarnoff symposium, 2009:1-6.

[11]任啸天, 徐晖, 王翔, 等.一种针对短码、周期长码直扩信号的扩频序列盲估计方法[J].宇航学报, 2011, 32 (12) :2596-2600.