一个简单的算法：复数求相位角

zouzhichao

先贴出复数结构定义体，很容易理解：


求相位，那不是很简单么？先贴第一个函数


当然，当你看到第二个函数，就会发现第一个函数有问题了：


这个时候，不会有问题了吧？其实不然，那你看看第三个函数：


当然，还有第四个函数，这个与第三个类似：


其实还有第五个函数，比第三第四又好一些：


源码：

1. typedef struct tagCOMPLEX {
   
2.         float real;
   
3.         float image;
   
4. }COMPLEX;
   
5. 
   
6. float cal(const COMPLEX* dat) {
   
7.         return atan(dat->image / dat->real);
   
8. }
   
9. 
   
10. float cal(const COMPLEX* dat) {
    
11.         float phase = atan(dat->image / dat->image);
    
12.         if (dat->real < 0.0f) {
    
13.                 phase += 3.1415926f;
    
14.         }
    
15.         if (phase > 3.1415926f) {
    
16.                 phase -= 3.1415926f * 2.0f;
    
17.         }
    
18.         return phase;
    
19. }
    
20. 
    
21. float cal(const COMPLEX* dat) {
    
22.         float amp, cos_x;
    
23.         amp = dat->real * dat->real + dat->image * dat->image;
    
24.         amp = inv_sqrt(amp);
    
25.         cos_x = dat->real * amp;
    
26.         amp = acos(cos_x);
    
27.         return dat->image > 0? amp : -amp;
    
28. }
    
29. 
    
30. float cal(const COMPLEX* dat) {
    
31.         float amp, sin_x;
    
32.         amp = dat->real * dat->real + dat->image * dat->image;
    
33.         amp = inv_sqrt(amp);
    
34.         sin_x = dat->image * amp;
    
35.         amp = asin(sin_x);
    
36.         if (dat->real < 0.0f) {
    
37.                 amp += 3.1415926f;
    
38.         }
    
39.         if (amp > 3.1415926f) {
    
40.                 amp -= 3.1415926f * 2.0f;
    
41.         }
    
42.         return amp;
    
43. }
    
44. 
    
45. float cal(const COMPLEX* dat) {
    
46.         float amp, cos_x, sin_x;
    
47.         amp = dat->real * dat->real + dat->image * dat->image;
    
48.         amp = inv_sqrt(amp);
    
49.         cos_x = dat->real * amp;
    
50.         sin_x = dat->image * amp;
    
51.         if (fabs(cos_x) < fabs(sin_x)) {
    
52.                 amp = acos(cos_x);
    
53.                 return sin_x > 0? amp : -amp;
    
54.         }
    
55.         amp = asin(sin_x);
    
56.         if (cos_x < 0.0f) {
    
57.                 amp += 3.1415926f;
    
58.         }
    
59.         if (amp > 3.1415926f) {
    
60.                 amp -= 3.1415926f * 2.0f;
    
61.         }
    
62.         return amp;
    
63. }

复制代码

有必要说明一下，后面几个函数用到了一个叫inv_sqrt的倒数平方根函数，这个函数有个快速算法，故居很多人都有收藏，就不贴出来

zouzhichao

一个快速浮点开方程序,据说能比float sqrt(x)快4倍

takashiki

你确定第一个是正确的？这个“算法”算出1+i和-1-i的相位一样??? 都是pi/4???

zouzhichao

takashiki 发表于 2017-3-12 14:25
你确定第一个是正确的？这个“算法”算出1+i和-1-i的相位一样??? 都是pi/4???

liujingwei

有个函数叫 atan2

zouzhichao

liujingwei 发表于 2017-3-12 14:39
有个函数叫 atan2

这是一个好函数

jlhgold

liujingwei 发表于 2017-3-12 14:39
有个函数叫 atan2

顶这个方位角，哈哈哈

NJ8888

多块？F407如果普通sqrt，168MHz时钟3微秒左右得到结果，arm_sqrt调用浮点硬件330ns左右

XA144F

atan2能算90度的整倍数么？这时的输入值接近无限大了！所以还是换成正余弦算避免90度的问题。

xuyapple

第二个函数：float phase = atan(dat->image / dat->image);
没写错吗？

lidapang

随机验证过可靠性么？

zouzhichao

xuyapple 发表于 2017-3-13 11:09
第二个函数：float phase = atan(dat->image / dat->image);
没写错吗？

写错了，image／real

takashiki

XA144F 发表于 2017-3-12 23:22
atan2能算90度的整倍数么？这时的输入值接近无限大了！所以还是换成正余弦算避免90度的问题。 ...

标准的IEEE 754中是有无穷大的，因此atan2能算90度的整倍数。比如：atan2(1, 0) = PI / 2

xycfwrj

给一个我写的基于dsp的近似计算的代码，其中_rcpsp是求倒数，需要替换
                        float den = fx + fy;
                        float num = fx - fy;
                        rcp = _rcpsp(den);   //= 1.f/den，也可以用快速算法
                        rcp = rcp * (2.f - den * rcp);  //

                        theta = 1.f - num * rcp;

                        if (fx<0.f)
                                theta = theta + 2.f;
                        if (fy<0.f)
                                theta = -theta;
                        theta = theta * PI/4;