用于DSP解码的相位计算方法 【技术领域】
本发明涉及一种相位计算方法，特别是一种用于DSP解码的相位计算方法。
背景技术
CORDIC(Coordinate Rotations Digital Computer)算法被广泛应用于数字信号处理(Digital Signal Processor，DSP)算法的硬件实现中。在收音机的立体声处理流程中，当经频率调制(FM，Frequency Modulation)的立体声中频信号被采样后，即输入至DSP芯片进行解码。在DSP解码过程中，相位计算将占用较多的时钟周期。相位计算即是由两路信号(I和Q)进行反三角函数计算，计算得出其相应相位。目前，所述相位计算皆采用CORDIC算法，
CORDIC算法由J.Volder于1959年提出，首先用于导航系统，使得矢量的旋转和定向运算不需要做查三角函数表、乘法、开方及反三角函数等复杂运算。CORDIC算法的核心思想是通过将一个复数与一常数序列相乘，以达到旋转该复数相位的目的。而为了用移位代替乘法，CORDIC算法一般选择的常数序列为2的整数幂序列。
CORDIC算法常用旋转模式，而正余弦计算则采用旋转模式。所述旋转模式是将某个向量旋转一个角度，具体公式(公式(1))与结果(公式(2))分别是：
xi+1=xi-yi*di*2-iyi+1=yi+xi*di*2-izi+1=zi-di*tan-1(2-i)]]>公式(1)
当Zi＜0时，di＝-1，否则di＝1，以及
xn=An[x0cosz0-y0sinz0]yn=An[y0cosz0+x0sinz0]zn=0An=Πn1+2-2i]]>公式(2)
上述CORDIC算法的缺点是：当提高相位计算精度时，相位的计算速度即会极大地下降。而相位计算的速度将极大地影响DSP的解码效率。
【发明内容】
有鉴于此，本发明的目的在于提供两种用于DSP解码的相位计算方法，使相位的计算速度和计算精度都得到极大的提升。
本发明是通过以下技术方案实现的：
一种用于DSP解码的相位计算方法，包括以下步骤：
A)设置查表分段数L的值，制作[0，π/4]区域的反正切表，并生成待查表的表格table[L+1]；
B)、计算索引index，查表得到相位值table[index]。
进一步地，所述步骤B)后进一步包括步骤C)、将坐标区域[-π，π]分为八个区域，将其他七个区域的相位映射到区域[0，π/4]，计算相应区域的相位。
进一步地，所述步骤A)的生成待查表的表格table[L+1]的公式是：
table[i]=(int){0.5+[0x100000*atan(iL)*4π]},]]>
其中，L为查表分段数，i为0，1，2，3...L。
进一步地，所述八个区域是[0，π/4]、(π/4，π/2)、(π/2，π3/4)、[π3/4，π]、(-π，-π3/4)、(-π3/4，-π/2)、[-π/2，-π/4]、[-π/4，0]。
进一步地，计算(π/4，π/2)的相位的公式是：
计算(π/2，π3/4)的相位的公式是：
计算[π3/4，π]的相位的公式是：θ＝π-α；
计算(-π，-π3/4)的相位的公式是：θ＝-π+α；
计算(-π3/4，-π/2)的相位的公式是：
计算[-π/2，-π/4]的相位的公式是：
计算[-π/4，0]的相位的公式是：θ＝-α，
其中，α为[0，π/4]的相位。
进一步地，所述步骤B)的计算索引的公式是：其中，X为待计算相位的点的横坐标，Y为待计算相位的点的纵坐标，L为查表分段数。
进一步地，所述查表分段数是自然数。
进一步地，所述查表分段数是2的n次幂，其中，n是整数。
本发明的另一技术方案如下：
一种用于DSP解码的相位计算方法，包括以下步骤：
A)设置查表分段数L的值，制作[0，π/4]区域的反正切表，并生成待查表的表格table[L+1]；
B)、计算索引index，查表得到相位值table[index]；
BC)、查表得到第一相位值table[index]和第二相位值table[index+1]，计算第三相位值。
进一步地，所述步骤BC)后进一步包括步骤C)、将坐标区域[-π，π]分为八个区域，将其他七个区域的相位映射到区域[0，π/4]，计算相应区域的相位。
进一步地，所述步骤A)地生成待查表的表格table[L+1]的公式是：
table[i]=(int){0.5+[0x100000*atan(iL)*4π]},]]>
其中，L为查表分段数，i为0，1，2，3...L。
进一步地，所述八个区域是[0，π/4]、(π/4，π/2)、(π/2，π3/4)、[π3/4，π]、(-π，-π3/4)、(-π3/4，-π/2)、[-π/2，-π/4]、[-π/4，0]。
进一步地，计算(π/4，π/2)的相位的公式是：
计算(π/2，π3/4)的相位的公式是：
计算[π3/4，π]的相位的公式是：θ＝π-α；
计算(-π，-π3/4)的相位的公式是：θ＝-π+α；
计算(-π3/4，-π/2)的相位的公式是：
计算[-π/2，-π/4]的相位的公式是：
计算[-π/4，0]的相位的公式是：θ＝-α，
其中，α为[0，π/4]的相位。
进一步地，所述步骤B)的计算索引的公式是：其中，X为待计算相位的点的横坐标，Y为待计算相位的点的纵坐标，L为查表分段数。
进一步地，所述步骤BC)的计算第三相位值的公式是：surplus＝Y*L-index*X，
其中，X为待计算相位的点的横坐标，Y为待计算相位的点的纵坐标，L为查表分段数，index是索引，table是待查表的表格。
进一步地，所述查表分段数是自然数。
进一步地，所述查表分段数是2的n次幂，其中，n是整数。
本发明的相位计算方法是不带插值查表法和带插值查表法，使相位的计算速度和计算精度都得到极大的提升。一次相位计算，不带插值查表法和带插值查表法占用的时钟周期约为CORDIC算法的一半。本发明的带插值查表法的计算误差呈规则分布并且足够小。不带插值查表法在查表分段数较大的情况下，计算误差也能达到足够小。
【附图说明】
图1是点(X，Y)的相位示意图。
图2是Δtanθ/Δθ-Δθ的曲线图。
图3是tanθ的坐标分区域图。
【具体实施方式】
请参阅图1，本发明的相位计算方法即是已知任一点的坐标(X，Y)，计算该点的相位θ的方法。
相位的计算可以使用查表法。查表法分为不带插值查表法和带插值查表法。以下对上述两种方法的原理、步骤和计算结果进行描述和对比。
请参阅图2，图2的横坐标是Δθ，纵坐标是Δtanθ/Δθ。在[0，π/4]范围内，Δtanθ/Δθ为单调增加。当Δθ≈0时，Δθ≈Δtanθ；当Δθ≈π/4时，Δtanθ/Δθ≈2。
对于不带插值查表法，由于介于两个关键点之间的任意点的函数值需用两个关键点中的某个点的函数值近似，因此，参阅图2可知，当Δθ≈π/4时，函数的斜率最大，点误差也最大。如果将0～π/4分为L段，L可为任意自然数。但为利于DSP处理，兼顾计算精度和占用存储器容量，一般将L取值为2的n次幂，即查表分段数L＝1，2，4，8，16，32，64，128，256，512，1024，...。表长为L+1，表中值分别为f(0)，f(1)，f(2)，...f(L-1)，f(L)，则不带插值查表法的误差的计算公式见公式(3)：
error=f(L-1)-f(L-2)2]]>公式(3)
因此，不带插值查表法的最大误差的计算公式见公式(4)：
errormax=π4*2*L]]>公式(4)
由此可见，不带插值查表法的最大误差和查表分段数L成反比，查表分段数L越小，最大误差越大。
不带插值查表法的计算步骤如下：
A)设置查表分段数L的值，制作[0，π/4]区域的atan(反正切)表，并生成待查表的表格(table表)；
根据公式(5)，制作待查表的表格：
table[i]=(int){0.5+[0x100000*atan(iL)*4π]}]]>公式(5)
其中，L为查表分段数，表长为L+1，i为0，1，2，3...L。
B)、计算索引index，查表得到相位值phase＝table[index]；
索引index的计算公式(6)如下：
index=Y*LX]]>公式(6)
其中，X为待计算相位的点的横坐标，Y为待计算相位的点的纵坐标，L为查表分段数。
上述步骤A)、B)是在区域[0，π/4]计算相位。因此，查表所得相位值phase是[0，π/4]上的相位。
C)、将坐标区域[-π，π]分为八个区域，将其他七个区域的相位映射到区域[0，π/4]，计算相应区域的相位。
根据图2所示的tanθ的单调性，将整个坐标分成八个区域，如图3所示。所述八个区域分别为[0，π/4]、(π/4，π/2)、(π/2，π3/4)、[π3/4，π]、(-π，-π3/4)、(-π3/4，-π/2)、[-π/2，-π/4]、[-π/4，0]。由于查表区域在区域1[0，π/4]，其他七个区域的角度需映射到区域1后，再进行查表。上述分法的优点在于：
(1)查表分段数可设置为较小值；
(2)由tan值的单调性可知，在[0，π/4]范围内误差比较小。
如区域1[0，π/4]的相位(角度)为α＝phase，可通过以下的公式计算其他七个区域相应的相位(角度)θ：
区域2(π/4，π/2)，角度
区域3(π/2，π3/4)，角度
区域4[π3/4，π]，角度θ＝π-α；
区域5(-π，-π3/4)，角度θ＝-π+α；
区域6(-π3/4，-π/2)，角度
区域7[-π/2，-π/4]，角度
区域8[-π/4，0]，角度θ＝-α。
当然，如果待计算相位的坐标点本身即位于区域1[0，π/4]，则不需执行步骤C)。
带插值查表法的计算方法在上述步骤B)、C)之间，添加如下步骤：
BC)、查表得到相位值table[index]和table[index+1]，根据公式(7)和公式(8)计算最终的相位值。
surplus＝Y*L-index*X 公式(7)
phase=table[index]+(table[index+1]-table[index])*surplusX]]>公式(8)
由上可知，带插值查表法与不带插值查表法相比，添加步骤BC)。
以下列举一具体实施例，对本发明作进一步阐述。
本具体实施例以coolFlux的DSP为硬件基础，作计算仿真。由于DSP计算仿真是定点处理，为减小误差，将π/2设置为2097152作为基准。当然，π/2也可以设置为其它任意值，本具体实施例并不限于此值。
现假设某一点的坐标值X＝1000，Y＝368，求其相位，具体的计算方法如下：
A1)、设置查表分段数L＝16，根据公式(5)制作待查表的表格，见表1；
表1
i＝0 i＝1 i＝2 i＝3 i＝4 i＝5 0 83335 166025 247456 327068 404378 i＝6 i＝7 i＝8 i＝9 i＝10 i＝11 478991 550604 619011 684085 745779 804107 i＝12 i＝13 i＝14 i＝15 i＝16 859131 910953 959702 1005524 1048576
如表1所示，0、83335、166025、247456、327068、404378、478991、550604、619011、684085、745779、804107、859131、910953、959702、1005524、1048576分别为i＝0、1、2、3、4、5、6、7、8、9、10、11、12、13、14、15、16时的待查表的值。
B1)、计算索引查表得到相位值phase＝table[index]＝table[5]＝404378；
404378即为不带插值查表法所得的相位值。
对于带插值查表法，继续执行步骤BC1)、surplus＝Y*L-index*X＝368＊16-5＊1000＝888；
phase=table[index]+(table[index+1]-table[index])*surplusX=]]>
table[5]+[(table[6]-table[5])*surplus]/1000=404378+[(478991-]]>
404378)*888]/1000=470634.]]>
由于X＝1000，Y＝368，点(X，Y)位于区域1[0，π/4]，因此，不需再执行步骤C)。
相位的实际精确值为470778，带插值查表法与不带插值查表法相比，精度要高得多。在本具体实施例中，查表分段数L＝16，不带插值查表法可以将查表分段数L的值设置得足够大，以减小误差。但L值越大，占用的存储容量越大，相位的计算速度也会随之下降。由此可见，带插值查表法在L较小的情况下，仍能保持较高的精度。
另外，再列举一些测试值，对比CORDIC算法、不带插值查表法和带插值查表法的优劣性。
表2为CORDIC算法的仿真结果；表3为不带插值查表法的仿真结果；表4为带插值查表法的仿真结果。其中，X，Y分别表示所求相位的输入横坐标X、纵坐标Y；期望值为相位的实际精确值；MIPS为仿真计算所需的机器语言指令数；仿真值为仿真结果；Δ为仿真误差。
表2CORDIC算法的仿真结果
表3不带插值查表法的仿真结果
表4带插值查表法的仿真结果
由表2、表3和表4可知，查表法与CORDIC算法相比，仿真计算所需的机器语言指令数MIPS要少得多。带插值查表法与不带插值查表法相比，仿真计算所需的MIPS稍多(增加步骤E))，但计算精度可以达到足够高。综合比较，带插值查表法使相位的计算速度和计算精度都得到极大的提升。
另外，由于查表分段数L可以设置为任意值，当然，可以减小至256，128，...16，8，4，2，1。对于不带插值查表法，当L＝16时，误差Δ为208，当表L＝8时，Δ为1317，误差太大。因此，当L值较小时，不建议使用该方法。
本发明的查表法可应用于FM立体声处理中。在软件解码器中，相位计算即是由两路信号(I和Q)进行反三角函数计算，计算得出其相应相位。相位计算的速度将极大地影响DSP的解码效率。FM立体声处理中的I和Q对应于公式和表中的X和Y，也即是横坐标X和纵坐标Y。当然，本发明并不限于FM立体声处理，还可以用于其它需进行相位计算的应用中。
以上介绍的仅仅是基于本发明的较佳实施例，并不能以此来限定本发明的范围。任何对本发明实施步骤作本技术领域内熟知的等同改变或替换均不超出本发明的揭露以及保护范围。