定频率相量提取方法
技术领域:
本发明涉及数字信号处理领域,特别涉及一种获取已知频率的正弦信号幅值和相位的数字信号处理方法。
背景技术:
正弦信号幅值和相位的计算是数字信号处理领域中的基本计算,该计算在振动信号处理、地球物理信号处理、雷达声纳信号处理及电力系统信号处理等绝大多数工程领域内获得了广泛的应用。
例如,振动的模态分析中,要计算振动的幅值和相位;地球物理学中地震勘探领域里,需要测量反射的地震波信号的幅值和相位;测距雷达在提取回波中频信号的功率谱,即需要测量回波的幅值;声纳中的接收换能器需要测量声波,也要分析声波信号的幅值;电力系统中系统的状态估计、稳定控制、故障分析及继电保护等领域中,相量的计算是十分重要的。
目前的计算方法有:过零检测法、离散傅里叶变换法(DFT)、基于DFT的若干算法,但是这些方法都有一些不足。
过零检测法:精度不高,易受谐波、噪声和非周期分量的影响,实时性不好。
DFT法:当被测信号中含有衰减高次谐波、衰减非整次谐波和衰减直流分量时,DFT法测量结果误差较大;当实际的频率偏离额定频率时,由于数据窗长的约束,会产生频率泄露,测量结果不稳定。
基于DFT的若干算法:针对DFT法受各类谐波和衰减直流分量影响的问题,目前有若干方法,但它们都对采样点提出了特殊要求,算法缺乏适应性;针对DFT的频率泄露问题,若干频率自适应采样方法被提出,这些方法本质上都有不可忽略的缺陷。
如电力系统同步相量测量中,若干依据信号频率来实时确定采样时刻的自适应方法,但它们也导致了变化的采样时刻与数据再同步所要求的确定时刻不易协调的矛盾,因此都难以满足高密度数据同步的要求。后来,又有学者提出基于定间隔采样的相量校正算法,虽然消除了频率泄露的影响,但该算法的动态延时很长,也无法用于瞬时暂态稳定控制。
另外,在地震勘探信号处理领域内,反射地震波往往不能转化为一个线性模型,因此直接使用目前这些线性模型都会带来误差;测距雷达在提取回波中频信号的功率谱时,由于受到系统噪声和杂波的影响,功率谱的幅值达不到系统设计的门限,造成雷达系统无法测距;而在电力系统中,特别是故障后的暂态信号,除了含有工频分量外,还有大量整次谐波、非整次谐波和衰减直流分量,目前的算法在提取工频相量时受这些分量的影响较大。
发明内容:
本发明要解决的技术问题有:提供一种不受整次谐波、非整次谐波和直流分量影响的幅值、相位计算方法;允许使用任意长度的数据窗(采样率确定时,至少应满足采样定理),并能解决频率泄露问题;计算量小,速度快,满足实时性要求。
本发明的目的是提供一种能克服频率泄露影响的、快速的、准确的计算正弦信号幅值和相位的方法,该方法可以测量信号中的任意一个或多个频率对应的幅值和相位。
为实现上述目的,本发明采用如下技术方案:
一种正弦信号幅值和相位的计算方法,先介绍单一频率的相量提取,然后介绍多个频率的相量提取。对于单一频率,具体包括以下步骤:
步骤一、通过信息获取或数据采集得到原始信号;
步骤二、对原始信号低通滤波,设滤波后的信号中仅含有p个频率分量;
步骤三、若步骤二中滤波后得到的信号是连续信号,还须对其进行采样,得到采样值x,设采样率为fs(满足采样定理);
步骤四、利用步骤三中得到的采样值x构成采样阵X;具体步骤如下:
取N个连续采样点x(i),i=1,2,…,N(N>2p),按以下形式构成采样阵:
注意:参数L不能小于2p,并且一般选则N/2≤L≤2N/3;而采样阵的行数没有特殊要求,因此该方法对数据窗长没有约束,若有冗余数据,则增加行数可以提高算法的抗噪声性能,通过后面的具体实施方案可以看到这一点。
步骤五、构造一个幅值为1、初相位为0的参考信号y=cos(2πf0t),令其频率与待测频率相同,并用对其进行采样(采样率与步骤三中相同,采样起始点为t=0),同理构造采样阵Y:
步骤六、将步骤五中的采样阵Y与步骤三中采样阵的Moore-Penrose广义逆X+相乘,得到N-L+1阶方阵:
A=YX+ (3);
这一步骤是通过构造出的仅含工频分量的矩阵Y与矩阵X+相乘,来实现仅保留信号X中的f0频率分量而滤除其余分量的目的。
步骤七、对步骤六中得到的方阵A进行矩阵变换,得到一个2阶方阵A′,保证它们工频对应的特征根相等;具体步骤如下:
首先说明两点:
③矩阵Y的秩为2。因为矩阵是由单一频率信号的采样值构成的,易知其秩为2。
④矩阵A的秩为2。由矩阵乘法知:r(A)=r(YX-1)≤r(Y)=2;又因为矩阵Y和X-1都含有频率f0,故其乘积A也含有该频率,因此有r(A)≥2;从而有r(A)=2。
其次,具体得到A′的步骤为:
⑤将矩阵A通过QR分解得到:
A=QR (4);
其中Q为正交矩阵,R仅有前两行有非零元,具体形式为:
其中n=N-L+1。
⑥可将矩阵R分解为:
R=ΛR′ (5);
其中Λ为对角阵,R′为上三角矩阵,且它们分别具有以下形式:
利用矩阵乘法可求得Λ和R′中各元素值:
故由R可求得唯一的ΛR′分解。
⑦用满秩矩阵R′对矩阵A进行相似变换,得到A′′,即:
A′′=R′AR′-1=R′QRR′-1=R′QΛ(R′R′-1)=R′QΛ (8);
由此可知,A′′仅有前两列有非零元,即具有如下形式:
⑧令矩阵A′′的二阶顺序主子式为矩阵A′,即:
根据特征方程的定义,分别写出A′′和A′的特征方程:
A′′:[λ2-(c11+c22)λ+c11c22-c12c21]λn-2=0 (11);
A′:λ2-(c11+c22)λ+c11c22-c12c21=0 (12);
由此可知,2阶方阵A′的特征根恰好是n阶方阵A′′的两个非零
特征根,且这两个特征根是频率f0对应的一对共轭复特征根。
步骤八、求解步骤七中2阶方阵A′的特征根λ1,2;
解得2阶方阵A′的特征根为:
λ1,2=a±jb (13);
其中
步骤九、根据步骤八得到的λ1,λ2,计算出X的幅值,即:
|X|=1/|λ1|=1/|λ2| (15);
步骤十、对参考信号y一个延时ΔT得到y′,然后得到采样阵Y′,重复步骤五~步骤九,得到共轭特征根λ1′,λ2′(λ1′,2=a′±jb′),结合λ1,λ2,计算得到X的相角;具体步骤如下:
由于仅进行一次计算得到的是一对共轭特征根λ1,λ2,其模值的倒数代表被测信号X的幅值,而辐角的倒数代表被测信号X与参考信号Y的夹角,但无法判断X与Y的超前或滞后关系。因此,需要进行两次测量,通过比较夹角的变化来确定信号X的相角。
④设信号X的相位为θ,根据信号Y的初相位为0,得到信号X与
信号Y的夹角|θ|为:
⑤根据信号Y′的初相位为2πf
0ΔT(即Y′超前Y的弧度),得到信号X与信号Y′的夹角|θ′|为:
⑥通过计算①、②中计算出的|θ|和|θ′|,查下表得到X的相位为θ。
对于多个频率(设待测k个相量的频点分别为:f1,f2,…,fk),处理过程与单一频率相似,不同的地方是:
①步骤五中,构造一个含有多个频率的参考信号,它们的幅值均为1、初相位均为0,即:y=cos(2πf1t)+cos(2πf2t)+…+cos(2πfkt);
②步骤七中,矩阵Y,R,Λ和A′均为2k阶方阵;
③步骤八、九和十中,求出的特征根均为k对共轭复根;
④特别地,当2k→N-L+1时,本方法趋向于未降阶模型。
发明主要提出了两种新方法:构造一个与待测频率同频的参考信号,对其采样得到采样阵,通过矩阵乘法滤除其它频率分量和直流分量的影响;一般地,相乘后得到的矩阵的阶数较高,但其中待测频率信号幅值和相位信息可以用一个低阶矩阵表示,本发明提出了一种从高阶方阵出抽取出一个低阶方阵的方法,保留了计算需要的频点信息,有效地提高了计算速度。
与现有技术相比,本发明主要具有以下三个优点:
1、消除了传统算法中对数据窗长的限制,没有传统DFT幅值相位测量中的频率泄露问题,并且在数据窗较短时依然可以得到准确的结果;
2、本方法消除了待测频率以外的所有谐波及直流分量的影响;
3、本方法计算量小,仅需要两次矩阵求逆。
附图说明:
图1是实现本发明所介绍方法的流程图;
图2是信号X、Y和Y′的矢量关系图;
图3是SNR=300dB,数据窗为15ms时,在40ms内滑窗计算的结果图;
图4是SNR=200dB,数据窗为15ms时,在40ms内滑窗计算的结果图;
图5是SNR=200dB,数据窗为22ms时,在40ms内滑窗计算的结果图。
具体实施方式:
下面结合附图对本发明做进一步详细描述:
下面结合附图的流程图,对本发明在电力系统工频相量计算中的应用做进一步详细描述。
一、在电力系统中,特别是故障后的暂态信号,除了含有工频分量外,还有大量整次谐波、非整次谐波和衰减直流分量,因此假设原始信号如式(1)所示,下面介绍用本发明的方法实现提取该信号中频率f0=50Hz的幅值和相位。
式中,ω=2πf0,τ=30ms,f0=50Hz,再加入方差为零,信噪比为300dB的白噪声。
二、低通滤波器的截止频率设为300Hz,采样时选10kHz作为采样频率。对x(t)进行低通滤波、采样保持和A/D转换后,得到采样序列x(k),k=1,2,…。
三、取150个连续采样点(对应15ms的数据窗)x(i),i=1,2,…,150,构造采样阵X,形式如下:
四、设置与待测频率相同、幅值为1、初相位为0的参考信号y=cos(2πf0t);同样地,用4kHz作为采样频率对其进行采样得到序列y(k),k=1,2,…;取前150个采样点构成采样阵Y,形式如下:
五、用采样阵X的Moore-Penrose广义逆X+右乘采样阵Y,得到63阶方阵:
A=YX+ (4);
六、从63阶方阵A中抽取2阶方阵A′,具体步骤如下:
①计算对方阵A进行QR分解后的矩阵R;
设63阶方阵A,Q和R的形式为:
由矩阵乘法及正交阵的性质依次可得:
其中j=3,4,…,63 (8);
②计算R=ΛR′分解后的矩阵R′;
设63阶方阵Λ和R′的形式为:
由矩阵乘法依次可得:
⑤用满秩矩阵R′对矩阵A进行相似变换,得到A′′:
⑥令A′为A′′的二阶主子式,即:
七、计算2阶方阵A′的特征根,得:
λ1,2=a±jb (13);
其中
八、计算信号x中频率为f0分量的幅值|X|为:
九、将参考信号y延时一个采样间隔得到y′,即ΔT=1/10kHz=0.1ms,再按四中方法构造对应的采样阵Y′;实际上,此时Y′就是用四中参考信号第2~151个采样点(共150个点)构造的采样阵,即:
十、用Y′代替Y,重复五、六、七,得到另一对共轭特征根λ1′,2=a′±jb′;
十一、计算信号x中频率为f0分量的相角θ;具体步骤如下:
①计算θ和θ′
②判断|θ|和|θ′|满足以下四个等式中的哪一个:
|θ|-|θ′|=2πf0ΔT (18)
|θ′|-|θ|=2πf0ΔT (19)
|θ′|+|θ|=2πf0ΔT (20)
|θ′|+|θ|=2π(1-f0ΔT) (21)
若式(18)或(20)成立,则θ=|θ|;否则,θ=-|θ|。
以上即为本发明提出的正弦信号幅值和相位的计算流程。
以采样间隔为步长,在两个周波内滑窗得到的测量结果如图三所示,其中幅值误差0。
若将信噪比变为200dB,其他均不变,则在两个周波内滑窗得到的测量结果如图四所示,其中最大幅值误差接近10%。
若保持信噪比为200dB,同时将数据窗变为22ms,即使用220个采样点,并设L=128,则在两个周波内滑窗得到的测量结果如图五所示,其中最大幅值误差不超过0.1%。
由仿真结果可知本发明提出的算法:
消除了待测频率以外所有谐波及直流分量的影响;对数据窗长没有限制,且计算结果的误差随数据窗长的增加而降低;计算量小,可用于实时监控。
以上所述,仅是本发明的较佳实施例而已,并非对本发明作任何形式上的限制,虽然本发明已以较佳实施例揭露如上,然而并非用以限定本发明,任何熟悉本专业的技术人员,在不脱离本发明技术方案范围内,当可利用上述揭示的方法及技术内容作出些许的更动或修饰为等同变化的等效实施例,但凡是未脱离本发明技术方案的内容,依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰,仍属于本发明技术方案的范围内。