正余弦的近似计算

G-Spider

1.FSIN，FCOS指令 2.Taylor series expansion 泰勒展开 http://en.wikipedia.org/wiki/Sine http://www.cnblogs.com/sun11086/archive/2009/03/20/1417944.html http://blog.csdn.net/ben_jiang/archive/2009/04/28/4133108.aspx (SSE2版本) 3.the CORDIC algorithm cordic_C语言.rar (15.01 KB, 下载次数: 1)

2011-2-18 18:35 上传

点击文件名下载附件

..... 有好的想法讨论一下吧。

G-Spider

以上给的TestCord.c 计算的是误差，改一下目的更明确了。 ------------------------------------------------------------------------------

1. /****************************************************************************
2. *
3. * Name: TestCord.c
4. *
5. * Synopsis: Test program for the Cordic module
6. *
7. * Copyright 1999 Grant R. Griffin
8. *
9. * The Wide Open License (WOL)
10. *
11. * Permission to use, copy, modify, distribute and sell this software and its
12. * documentation for any purpose is hereby granted without fee, provided that
13. * the above copyright notice and this license appear in all source copies.
14. * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF
15. * ANY KIND. See http://www.dspguru.com/wol.htm for more information.
16. *
17. *****************************************************************************/
18. //win7自带计算器得到：
19. //sin(11°)=0.19080899537654481240514048795839
20. //cos(11°)=0.98162718344766395349650489981814
23. #include <stdio.h>
24. #include "cordic.h"
27. /****************************************************************************/
28. int main(void)
29. {
31. #define DEGREES_TO_RADS(x) (x / 180.0 * PI) /*x为角度,DEGREES_TO_RADS(x)转成弧度*/
33. int max_L; /*精度依赖于max_L的值，如取45*/
34. double Cos_out, Sin_out;
36. for (;;) {
38. do {
39. max_L = 0;
40. printf("\nEnter maximum L >= 3, (0 to quit) : ");
41. scanf("%i", &max_L);
42. if (max_L == 0) {
43. exit(0);
44. }
45. } while (max_L < 3);
47. cordic_construct(max_L);
49. cordic_get_cos_sin(DEGREES_TO_RADS(11.0), &Cos_out, &Sin_out);
50. /*以上函数用于计算cos(11°) ，sin(11°) */
52. printf("Cos_out = %20.15lf dB, Sin_out = %20.15lf dB\n",Cos_out,Sin_out);
54. cordic_destruct();
56. }
58. return 0;
59. }

复制代码

G-Spider

cordic.c中 ... mp_cordic_table = (CORDIC_TABLE *) calloc(max_L + 1, sizeof(CORDIC_TABLE)); ... K = 1.0f; for (L = 0; L <= max_L; L++) { mp_cordic_table[L].K = K; mp_cordic_table[L].phase_rads = (double) atan(K); //arctan(k)计算 K *= 0.5f; //因为C doesn't let you shift "doubles" } 对于固定的max_L ，如max_L=45（可以得到14～15位的精度）时，可以给一个固定的mp_cordic_table表空间而不用动态分配，mp_cordic_table[L]也是确定的，可以预先给出，所以也省去了arctan(k)和 K *= 0.5f的计算。

G-Spider

the CORDIC algorithm, which can calculate the trigonometric functions of sin, cos, magnitude and phase using multiplying factors which are powers of two. Therefore, the multiplies can be implemented in hardware as shift/adds. (Since C doesn't let you shift "doubles", the shift/add multiplies appear here as literal multiplies, but you can get idea.) 由于c中无法对double数据进行移位操作，则乘法运算不可避免，以即下面的程序片段

1. /* rotate using "1 + jK" factors */
2. for (L = 0; L < m_max_L; L++) {
3. K = K_phase_rads[L][0];
4. phase_rads = K_phase_rads[L][1];
5. tmp_I = I;
6. if (desired_phase_rads - acc_phase_rads < 0.0) {
7. /* do negative roation */
8. I += Q * K;
9. Q -= tmp_I * K;
10. acc_phase_rads -= phase_rads;
11. } else {
12. /* do positive rotation */
13. I -= Q * K;
14. Q += tmp_I * K;
15. acc_phase_rads += phase_rads;
16. }
17. }

复制代码

中的Q*K 和tmp_I*K不能通过移位实现，这样一来若取 m_max_L=46 比泰勒来的还要慢，在FPGA中是可行的，所以说“you can get idea”。 改写最初的程序： cordic.zip (118.32 KB, 下载次数: 0)

2011-2-19 13:50 上传

点击文件名下载附件

wayne

可以看看fastfunlib http://code.google.com/p/fastfunlib/

G-Spider

试了下fsincos指令，用在fft中应该可以。

1. .386
2. .model flat,stdcall
3. option casemap:none
5. ;-----------------------------------------------------------
6. .data
7. ConsPi_180 real8 0.017453292519943295769
8. ;***********************************************************
9. .code
10. ;-----------------------------------------------------------
11. _Degree2SinCos proc deg:real4,lpsin:dword,lpcos:dword
12. mov esi,lpcos
13. mov edi,lpsin
14. finit
15. fld real4 ptr[deg]
16. fld qword ptr [ConsPi_180]
17. fmul st(0),st(1)
18. fsincos
19. fstp QWORD ptr[esi]
20. fst QWORD ptr[edi]
21. ret
23. _Degree2SinCos endp
25. _Rad2SinCos proc rad:real8,lpsin:dword,lpcos:dword
26. mov esi,lpcos
27. mov edi,lpsin
28. finit
29. fld real8 ptr[rad]
30. fsincos
31. fstp QWORD ptr[esi]
32. fst QWORD ptr[edi]
33. ret
35. _Rad2SinCos endp
37. end

复制代码

c中调用一下：

1. #include <stdio.h>
2. #define PI 3.1415926535897932
4. #ifdef __cplusplus
5. #include <windows.h>
6. extern "C" {
7. #endif
8. void _stdcall _Degree2SinCos(float degree,double *sin,double *cos);
9. void _stdcall _Rad2SinCos(double Rad,double *sin,double *cos);
11. #ifdef __cplusplus
12. }
13. #endif
15. int main()
16. {
17. float deg=30.0;
18. double rad=30.0*PI/180.0;
19. double sin,cos;
21. _Degree2SinCos(deg,&sin,&cos);
22. printf("_Degree2SinCos:\ncos(%.4f)=%18.16lf\nsin(%.4f)=%18.16lf\n",deg,cos,deg,sin);
24. _Rad2SinCos(rad,&sin,&cos);
25. printf("_Rad2SinCos:\ncos(%.4f)=%18.16lf\nsin(%.4f)=%18.16lf\n",rad,cos,rad,sin);
26. //system("PAUSE");
27. getch();
28. return 0;
29. }

复制代码

_Degree2SinCos: cos(30.0000)=0.8660254037844386 sin(30.0000)=0.5000000000000000 _Rad2SinCos: cos(0.5236)=0.8660254037844387 sin(0.5236)=0.4999999999999999 唯一不爽的是上面代码生成的可执行文件，360竟然报毒，没法蛋定了，你懂的。 test_src.zip (22.79 KB, 下载次数: 0)

2011-2-22 12:15 上传

点击文件名下载附件

wayne

6# G-Spider python 3.2 昨天释出，我从官网 下载了安装的时候，360还很积极的拦截呢，说是

发现木马：Win32/Trojan.e8c 详细描述： 木马名称：Win32/Trojan.e8c 所在路径：G:\Documents and Settings\yoyo\桌面\stardict\python-3.2.msi

G-Spider

7# wayne “永久免费”背后的生财之道。。。你懂的。。。

G-Spider

上面的fsincos与math.h 中的sin(),cos()时间测试。 j=10~30 N=2^j，对于每次的j++ ，计算N0=N/4个sin，cos值的时间，随j增加，时间减少接近于40%。

1. Elapsed time:
2. _Rad2SinCos Std_SinCos (seconds)
3. 0.00000000 0.00000000
4. 0.00000000 0.00000000
5. 0.00000000 0.00000000
6. 0.00000000 0.00000000
7. 0.00000000 0.00000000
8. 0.00000000 0.00000000
9. 0.01600000 0.01500000
10. 0.01600000 0.01600000
11. 0.03100000 0.04700000
12. 0.06200000 0.07800000
13. 0.12500000 0.18800000
14. 0.23400000 0.37500000
15. 0.43800000 0.73400000
16. 0.87500000 1.46900000
17. 1.81200000 3.00000000
18. 3.86000000 6.25000000
19. 7.04700000 11.84300000
20. 14.06300000 23.81200000
21. 28.17200000 47.43800000
22. 56.42200000 94.78100000
23. 请按任意键继续. . .

复制代码

测试代码：

1. #include <stdio.h>
2. #include <math.h>
3. #include <time.h>
4. #define PI 3.1415926535897932
5. #define _2PI (2.0*PI)
6. #ifdef __cplusplus
7. #include <windows.h>
8. extern "C" {
9. #endif
10. void _stdcall _Degree2SinCos(float degree,double *sin,double *cos);
11. void _stdcall _Rad2SinCos(double Rad,double *sin,double *cos);
13. #ifdef __cplusplus
14. }
15. #endif
17. int main()
18. {
19. double cos0,sin0,cos1,sin1;
20. double t0,t1;
21. double i,N,N0,N1,theta;//1/N
22. int j;
24. printf("Elapsed time: \n_Rad2SinCos Std_SinCos (seconds)\n");
25. for(j=10;j<30;j++)
26. {
27. N=(unsigned int)1<<j;
28. N0=N/4.0;
29. N1=1.0/N;
30. //=======================================================
31. t0=clock();
32. for(i=0;i<=N0;i++)
33. {
34. theta=(double)i*N1*_2PI;
35. _Rad2SinCos(theta,&sin0,&cos0);
36. //printf("%18.16lf %18.16lf %18.16lf\n",theta,sin0,cos0);
37. }
38. printf("%.8f ",(clock()-t0)/CLOCKS_PER_SEC);
39. //=======================================================
40. t1=clock();
41. for(i=0;i<=N0;i++)
42. {
43. theta=(double)i*N1*_2PI;
44. sin1=sin(theta);
45. cos1=cos(theta);
46. //printf("%18.16lf %18.16lf %18.16lf\n",theta,sin1,cos1);
47. }
48. printf("%.8f \n",(clock()-t1)/CLOCKS_PER_SEC);
49. //=======================================================
50. }
52. system("PAUSE");
53. return 0;
54. }

复制代码

G-Spider

9# G-Spider 上面的弧度都是倍数关系，看来可以用倍角公式之类的解决。 sin(2x)=2sin(x)*cos(x); cos(2x)=1 -2(sinx)^2; sin(x+y)=sin(x)cos(y)+cos(x)sin(y); cos(x+y)=cos(x)cos(y)-sin(x)sin(y); 对于弧度x,2x,3x,4x,5x,6x,7x,........ 偶数项2x,4x,6x,..... 奇数项3x ,5x ,7x,....改成(2x+x), (4x+x),(6x+x),...... 于是 sin(a+x)=sin(a)cos(x)+cos(a)sin(x)=c0*sin(a)+c1*cos(a); cos(a+x)=cos(a)cos(x)-sin(a)sin(x)=c0*cos(a)-c1*sin(a); 奇数项的sin,cos值只与前一项的正余弦有关。 偶数项的sin,cos值由n/2的正余弦得到。 按奇偶交替(数组下标定位)应该可以得到结果。 这样对于弧度x,2x,3x,4x,........所有的sin,cos值均可由sin(x),cos(x)得出，不知精度如何，有时间试试。