正余弦的近似计算

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
    
21. 
    
22. 
    
23. #include <stdio.h>
    
24. #include "cordic.h"
    
25. 
    
26. 
    
27. /****************************************************************************/
    
28. int main(void)
    
29. {
    
30.                  
    
31.   #define DEGREES_TO_RADS(x)  (x / 180.0 * PI) /*x为角度,DEGREES_TO_RADS(x)转成弧度*/
    
32. 
    
33.   int     max_L;                              /*精度依赖于max_L的值，如取45*/
    
34.   double  Cos_out, Sin_out;
    
35. 
    
36.   for (;;) {
    
37. 
    
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);
    
46. 
    
47.     cordic_construct(max_L);
    
48. 
    
49.     cordic_get_cos_sin(DEGREES_TO_RADS(11.0), &Cos_out, &Sin_out);
    
50.     /*以上函数用于计算cos(11°) ，sin(11°) */
    
51.    
    
52.     printf("Cos_out = %20.15lf dB, Sin_out = %20.15lf dB\n",Cos_out,Sin_out);
    
53.    
    
54.     cordic_destruct();
    
55.    
    
56.   }
    
57. 
    
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
   
4. 
   
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
    
22. 
    
23. _Degree2SinCos endp
    
24. 
    
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
    
34. 
    
35. _Rad2SinCos endp
    
36. 
    
37. end

复制代码

c中调用一下：

1. #include <stdio.h>
   
2. #define PI 3.1415926535897932
   
3. 
   
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);
   
10. 
    
11. #ifdef __cplusplus
    
12. }
    
13. #endif
    
14. 
    
15. int main()
    
16. {
    
17.     float  deg=30.0;
    
18.     double rad=30.0*PI/180.0;
    
19.     double sin,cos;
    
20.    
    
21.     _Degree2SinCos(deg,&sin,&cos);
    
22.     printf("_Degree2SinCos:\ncos(%.4f)=%18.16lf\nsin(%.4f)=%18.16lf\n",deg,cos,deg,sin);
    
23.    
    
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);
    
12. 
    
13. #ifdef __cplusplus
    
14. }
    
15. #endif
    
16. 
    
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;
    
23.    
    
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.     }
    
51. 
    
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)得出，不知精度如何，有时间试试。