模拟信号求解相位差（2）

本文基于《正弦量频率相位测量的新方法——高金峰》所提出的方法。此方法对上篇博客中所提的方法进行了改进，计算精度有所提升，但是对于被测信号进行了限制，具体见下文分析

1.理论基础

假设输入的两个模拟信号的表达式分别为：
μ1(t)=U1msin(ωt+θ1) \mu_1(t)=\boldsymbol{U_{1m}}\sin(\omega t+\theta_1)
μ2(t)=U2msin(ωt+θ2) \mu_2(t)=\boldsymbol{U_{2m}}\sin(\omega t+\theta_2)
从表达式中可得：两个信号的频率为 ω2π \cfrac{\omega}{2\pi}，相位差为 θ=θ1−θ2 \theta=\theta_1-\theta_2
设观察时刻 t=0 t=0，则以 t=0 t=0为中心，以时间间隔为 △t \vartriangle t对 μ1(t) \mu_1(t)和 μ2(t) \mu_2(t)进行同时采样可得采样序列：
⋯μ1(−3△t),μ1(−2△t),μ1(−△t),μ1(0),μ1(△t),μ1(2△t),μ1(3△t),⋯ \cdots\mu_1(-3\vartriangle t),\mu_1(-2\vartriangle t),\mu_1(-\vartriangle t),\mu_1(0),\mu_1(\vartriangle t),\mu_1(2\vartriangle t),\mu_1(3\vartriangle t),\cdots
⋯μ2(−3△t),μ2(−2△t),μ2(−△t),μ2(0),μ2(△t),μ2(2△t),μ2(3△t),⋯ \cdots\mu_2(-3\vartriangle t),\mu_2(-2\vartriangle t),\mu_2(-\vartriangle t),\mu_2(0),\mu_2(\vartriangle t),\mu_2(2\vartriangle t),\mu_2(3\vartriangle t),\cdots
结合信号的表达式我们可以得到
μ1(0)=U1msin(θ1)(1)μ1(△t)=U1msin(ω△t+θ1)(2)μ1(−△t)=U1msin(−ω△t+θ1)(3)μ1(2△t)=U1msin(2ω△t+θ1)(4)μ1(−2△t)=U1msin(−2ω△t+θ1)(5) \mu_1(0)=\boldsymbol{U_{1m}}\sin(\theta_1) (1)\\\mu_1(\vartriangle t)=\boldsymbol{U_{1m}}\sin(\omega\vartriangle t+\theta_1)(2)\\\mu_1(-\vartriangle t)=\boldsymbol{U_{1m}}\sin(-\omega\vartriangle t+\theta_1)(3)\\\mu_1(2\vartriangle t)=\boldsymbol{U_{1m}}\sin(2\omega\vartriangle t+\theta_1)(4)\\\mu_1(-2\vartriangle t)=\boldsymbol{U_{1m}}\sin(-2\omega\vartriangle t+\theta_1)(5)
由和差化积公式
sina+sinb=2sina+b2cosa−b2sina−sinb=2cosa+b2sina−b2 \sin a+\sin b=2\sin\cfrac{a+b}{2}\cos\cfrac{a-b}{2}\\\sin a-\sin b=2\cos\cfrac{a+b}{2}\sin\cfrac{a-b}{2}
将 (2)+(3) (2)+(3)可得
2U1msinθ1cos(ω△t)=μ1(△t)+μ1(−△t)(6) 2\boldsymbol{U_{1m}}\sin\theta_1\cos(\omega\vartriangle t)=\mu_1(\vartriangle t)+\mu_1(-\vartriangle t)(6)
将 (2)−(3) (2)-(3)可得
2U1mcosθ1sin(ω△t)=μ1(△t)−μ1(−△t)(7) 2\boldsymbol{U_{1m}}\cos\theta_1\sin(\omega\vartriangle t)=\mu_1(\vartriangle t)-\mu_1(-\vartriangle t)(7)
将 (4)+(5) (4)+(5)可得
2U1msinθ1cos(2ω△t)=μ1(2△t)+μ1(−2△t)(8) 2\boldsymbol{U_{1m}}\sin\theta_1\cos(2\omega\vartriangle t)=\mu_1(2\vartriangle t)+\mu_1(-2\vartriangle t)(8)
由 (1) (1)和 (6)⇒cos(ω△t)=μ1(△t)+μ1(−△t)2μ1(0)⇒ω=arccos(μ1(△t)+μ1(−△t)2μ1(0))/△t(9) (6)\Rightarrow\cos(\omega\vartriangle t)=\cfrac{\mu_1(\vartriangle t)+\mu_1(-\vartriangle t)}{2\mu_1(0)}\\\Rightarrow\omega=\arccos(\cfrac{\mu_1(\vartriangle t)+\mu_1(-\vartriangle t)}{2\mu_1(0)})/\vartriangle t(9)
这样当 μ1(0)≠0 \mu_1(0)\ne0时，我们就可以通过三个采样点的数值计算出模拟信号的角频率 ω \omega
当 μ1(0)≠0 \mu_1(0)\ne0且 μ1(△t)−μ1(−△t)≠0 \mu_1(\vartriangle t)-\mu_1(-\vartriangle t)\ne0时，由 (6)(7)⇒ \cfrac{(6)}{(7)}\Rightarrow
tanθ1tan(ω△t)=μ1(△t)+μ1(−△t)μ1(△t)−μ1(−△t)⇒tanθ1=μ1(△t)+μ1(−△t)μ1(△t)−μ1(−△t)tan(ω△t)(10) \cfrac{\tan\theta_1}{\tan(\omega\vartriangle t)}=\cfrac{\mu_1(\vartriangle t)+\mu_1(-\vartriangle t)}{\mu_1(\vartriangle t)-\mu_1(-\vartriangle t)}\\\Rightarrow\tan\theta_1=\cfrac{\mu_1(\vartriangle t)+\mu_1(-\vartriangle t)}{\mu_1(\vartriangle t)-\mu_1(-\vartriangle t)}\tan(\omega\vartriangle t)(10)
由 (1) (1)和 (8)⇒cos(ω△t)=μ1(2△t)+μ1(−2△t)2μ1(0)(11) (8)\Rightarrow\cos(\omega\vartriangle t)=\cfrac{\mu_1(2\vartriangle t)+\mu_1(-2\vartriangle t)}{2\mu_1(0)}(11)
又由三角公式 tana=1−cos2a1+cos2a−−−−−−−−−√ \tan a=\sqrt{\cfrac{1-\cos 2a}{1+\cos 2a}}(取正值)可得
tanθ1=μ1(△t)+μ1(−△t)μ1(△t)−μ1(−△t)2μ1(0)−μ1(2△t)−μ1(−2△t)2μ1(0)+μ1(2△t)+μ1(−2△t)−−−−−−−−−−−−−−−−−−−−−−−−−−√(12) \tan\theta_1=\cfrac{\mu_1(\vartriangle t)+\mu_1(-\vartriangle t)}{\mu_1(\vartriangle t)-\mu_1(-\vartriangle t)}\sqrt{\cfrac{2\mu_1(0)-\mu_1(2\vartriangle t)-\mu_1(-2\vartriangle t)}{2\mu_1(0)+\mu_1(2\vartriangle t)+\mu_1(-2\vartriangle t)}}(12)
令 t11=μ1(△t)+μ1(−△t)μ1(△t)−μ1(−△t)t12=2μ1(0)−μ1(2△t)−μ1(−2△t)2μ1(0)+μ1(2△t)+μ1(−2△t)−−−−−−−−−−−−−−−−−−−−−−−−−−√t21=μ2(△t)+μ2(−△t)μ2(△t)−μ2(−△t)t22=2μ2(0)−μ2(2△t)−μ2(−2△t)2μ2(0)+μ2(2△t)+μ2(−2△t)−−−−−−−−−−−−−−−−−−−−−−−−−−√ t_{11}=\cfrac{\mu_1(\vartriangle t)+\mu_1(-\vartriangle t)}{\mu_1(\vartriangle t)-\mu_1(-\vartriangle t)}\\t_{12}=\sqrt{\cfrac{2\mu_1(0)-\mu_1(2\vartriangle t)-\mu_1(-2\vartriangle t)}{2\mu_1(0)+\mu_1(2\vartriangle t)+\mu_1(-2\vartriangle t)}}\\t_{21}=\cfrac{\mu_2(\vartriangle t)+\mu_2(-\vartriangle t)}{\mu_2(\vartriangle t)-\mu_2(-\vartriangle t)}\\t_{22}=\sqrt{\cfrac{2\mu_2(0)-\mu_2(2\vartriangle t)-\mu_2(-2\vartriangle t)}{2\mu_2(0)+\mu_2(2\vartriangle t)+\mu_2(-2\vartriangle t)}}
所以两个模拟信号的初始相位为
θ1=arctan(t11t12)θ2=arctan(t21t22) \theta_1=\arctan(t_{11}t_{12})\\\theta_2=\arctan(t_{21}t_{22})
最终我们得到两个信号的相位差为 θ=∣θ1−θ2∣ \theta=\mid\theta_1-\theta_2\mid

2.流程图

Created with Raphaël 2.1.0 开始对两路信号进行A/D 获得采样序列采样序列中两两的间隔是否过小? 改变采样时间间隔两个中心点的值是否为0? 计算t11、t12、t21、t22 计算两个模拟信号的初始相位计算角频率计算相位差输出相位差和角频率结束 yes no yes no

3.MATLAB仿真实现

clear all;F=1;%被采样信号频率1Hz
Ts=0.1;%采样间隔%对下列信号进行处理
u1=inline('2*sin(2*pi*t)');
u2=inline('3*sin(2*pi*t+pi/3)');n=-2:2;
nTs=n*Ts;%采样序列
f1=u1(F*nTs);
f2=u2(F*nTs);if(f1(3)==0)angle1=0;
elset11=(f1(4)+f1(2))/(f1(4)-f1(2));t12=sqrt((2*f1(3)-f1(5)-f1(1))/(2*f1(3)+f1(5)+f1(1)));angle1=atan(t11*t12);%信号1的相位                                                                                              w=(f1(2)+f1(4))/(2*f1(3));
endif(f2(3)==0)angle2=0;
elset21=(f2(4)+f2(2))/(f2(4)-f2(2));t22=sqrt((2*f2(3)-f2(5)-f2(1))/(2*f2(3)+f2(5)+f2(1)));angle2=atan(t21*t22);%信号2的相位                                                                                                  w=(f2(2)+f2(4))/(2*f2(3));
endangle=abs(angle1-angle2)%两个信号的相位差                                                   W=acos(w)/Ts;%求出角速度W
F=W/(2*pi)%求出模拟信号的频率
T=1/F;%求出周期

4.思考

4.1对于信号的要求

从实现原理中我们可以看出，此方法对于初始相位为0的模拟信号无法进行求解，我用自己的理解已将此部分的解决方法放到代码中，但正确性还待验证。

4.2异步采样

此方法的前提是对于两路信号采样的时间相同，即同步采样。而现实情况中实现理想情况的同步采样是较难实现的。文中给出了校验的方法。在频率 ω \omega确定的形况下，只要选定一个信号的 t=0 t=0的位置，另一路信号 t=0 t=0的位置就可以由第一路信号 t=0 t=0的位置确定。若两路信号的采集时间相差 △t2 \cfrac{\vartriangle t}{2}，且第二路信号 t=0 t=0的采集点在第一路采集点之后，则计算两路信号相位差的公式由 θ=∣θ1−θ2∣ \theta=\mid\theta_1-\theta_2\mid变为 θ=∣θ1−θ2+ω△t2∣ \theta=\mid\theta_1-\theta_2+\cfrac{\omega\vartriangle t}{2}\mid
以上就是我对这篇文章的理解和实验，有不足之处还望大家多多指教。