论坛首页 › DSP › 一个简单的算法：复数求相位角

zouzhichao 发表于 2017-3-12 14:02:26

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

本帖最后由 zouzhichao 于 2017-3-12 14:16 编辑

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


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


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


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


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


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


源码：
typedef struct tagCOMPLEX {
        float real;
        float image;
}COMPLEX;

float cal(const COMPLEX* dat) {
        return atan(dat->image / dat->real);
}

float cal(const COMPLEX* dat) {
        float phase = atan(dat->image / dat->image);
        if (dat->real < 0.0f) {
                phase += 3.1415926f;
        }
        if (phase > 3.1415926f) {
                phase -= 3.1415926f * 2.0f;
        }
        return phase;
}

float cal(const COMPLEX* dat) {
        float amp, cos_x;
        amp = dat->real * dat->real + dat->image * dat->image;
        amp = inv_sqrt(amp);
        cos_x = dat->real * amp;
        amp = acos(cos_x);
        return dat->image > 0? amp : -amp;
}

float cal(const COMPLEX* dat) {
        float amp, sin_x;
        amp = dat->real * dat->real + dat->image * dat->image;
        amp = inv_sqrt(amp);
        sin_x = dat->image * amp;
        amp = asin(sin_x);
        if (dat->real < 0.0f) {
                amp += 3.1415926f;
        }
        if (amp > 3.1415926f) {
                amp -= 3.1415926f * 2.0f;
        }
        return amp;
}

float cal(const COMPLEX* dat) {
        float amp, cos_x, sin_x;
        amp = dat->real * dat->real + dat->image * dat->image;
        amp = inv_sqrt(amp);
        cos_x = dat->real * amp;
        sin_x = dat->image * amp;
        if (fabs(cos_x) < fabs(sin_x)) {
                amp = acos(cos_x);
                return sin_x > 0? amp : -amp;
        }
        amp = asin(sin_x);
        if (cos_x < 0.0f) {
                amp += 3.1415926f;
        }
        if (amp > 3.1415926f) {
                amp -= 3.1415926f * 2.0f;
        }
        return amp;
}

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

zouzhichao 发表于 2017-3-12 14:20:49

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

takashiki 发表于 2017-3-12 14:25:10

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

zouzhichao 发表于 2017-3-12 14:34:10

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

liujingwei 发表于 2017-3-12 14:39:57

有个函数叫 atan2

zouzhichao 发表于 2017-3-12 14:44:58

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

这是一个好函数

jlhgold 发表于 2017-3-12 14:46:04

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

顶这个方位角，哈哈哈

NJ8888 发表于 2017-3-12 15:16:19

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

XA144F 发表于 2017-3-12 23:22:46

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

xuyapple 发表于 2017-3-13 11:09:15

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

lidapang 发表于 2017-3-13 11:17:31

随机验证过可靠性么？

zouzhichao 发表于 2017-3-13 11:18:32

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

写错了，image／real

takashiki 发表于 2017-3-14 07:05:58

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

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

xycfwrj 发表于 2017-4-17 18:19:20

本帖最后由 xycfwrj 于 2017-4-17 18:31 编辑

给一个我写的基于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;

查看完整版本: 一个简单的算法：复数求相位角