正余弦的近似计算
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
.....
有好的想法讨论一下吧。 以上给的TestCord.c 计算的是误差，改一下目的更明确了。
------------------------------------------------------------------------------/****************************************************************************
*
* Name: TestCord.c
*
* Synopsis: Test program for the Cordic module
*
* Copyright 1999Grant R. Griffin
*
* The Wide Open License (WOL)
*
* Permission to use, copy, modify, distribute and sell this software and its
* documentation for any purpose is hereby granted without fee, provided that
* the above copyright notice and this license appear in all source copies.
* THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF
* ANY KIND. See http://www.dspguru.com/wol.htm for more information.
*
*****************************************************************************/
//win7自带计算器得到：
//sin(11°)=0.19080899537654481240514048795839
//cos(11°)=0.98162718344766395349650489981814
#include <stdio.h>
#include "cordic.h"
/****************************************************************************/
int main(void)
{
#define DEGREES_TO_RADS(x)(x / 180.0 * PI) /*x为角度,DEGREES_TO_RADS(x)转成弧度*/
int max_L; /*精度依赖于max_L的值，如取45*/
doubleCos_out, Sin_out;
for (;;) {
do {
max_L = 0;
printf("\nEnter maximum L >= 3, (0 to quit) : ");
scanf("%i", &max_L);
if (max_L == 0) {
exit(0);
}
} while (max_L < 3);
cordic_construct(max_L);
cordic_get_cos_sin(DEGREES_TO_RADS(11.0), &Cos_out, &Sin_out);
/*以上函数用于计算cos(11°) ，sin(11°) */
printf("Cos_out = %20.15lf dB, Sin_out = %20.15lf dB\n",Cos_out,Sin_out);
cordic_destruct();
}
return 0;
} 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.K = K;
mp_cordic_table.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也是确定的，可以预先给出，所以也省去了arctan(k)和K *= 0.5f的计算。 the CORDIC algorithm, which can calculate thetrigonometric 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 youshift "doubles", the shift/add multiplies appear here as literalmultiplies, but you can get idea.)
由于c中无法对double数据进行移位操作，则乘法运算不可避免，以即下面的程序片段/* rotate using "1 + jK" factors */
for (L = 0; L < m_max_L; L++) {
K = K_phase_rads;
phase_rads = K_phase_rads;
tmp_I = I;
if (desired_phase_rads - acc_phase_rads < 0.0) {
/* do negative roation */
I += Q * K;
Q -= tmp_I * K;
acc_phase_rads -= phase_rads;
} else {
/* do positive rotation */
I -= Q * K;
Q += tmp_I * K;
acc_phase_rads += phase_rads;
}
}中的Q*K和tmp_I*K不能通过移位实现，这样一来若取 m_max_L=46 比泰勒来的还要慢，在FPGA中是可行的，所以说“you can get idea”。
改写最初的程序：
可以看看fastfunlib
http://code.google.com/p/fastfunlib/ 试了下fsincos指令，用在fft中应该可以。.386
.model flat,stdcall
option casemap:none
;-----------------------------------------------------------
.data
ConsPi_180 real8 0.017453292519943295769
;***********************************************************
.code
;-----------------------------------------------------------
_Degree2SinCos proc deg:real4,lpsin:dword,lpcos:dword
mov esi,lpcos
mov edi,lpsin
finit
fld real4 ptr
fld qword ptr
fmul st(0),st(1)
fsincos
fstp QWORD ptr
fstQWORD ptr
ret
_Degree2SinCos endp
_Rad2SinCos proc rad:real8,lpsin:dword,lpcos:dword
mov esi,lpcos
mov edi,lpsin
finit
fld real8 ptr
fsincos
fstp QWORD ptr
fstQWORD ptr
ret
_Rad2SinCos endp
endc中调用一下：#include <stdio.h>
#define PI 3.1415926535897932
#ifdef __cplusplus
#include <windows.h>
extern "C" {
#endif
void _stdcall _Degree2SinCos(float degree,double *sin,double *cos);
void _stdcall _Rad2SinCos(double Rad,double *sin,double *cos);
#ifdef __cplusplus
}
#endif
int main()
{
floatdeg=30.0;
double rad=30.0*PI/180.0;
double sin,cos;
_Degree2SinCos(deg,&sin,&cos);
printf("_Degree2SinCos:\ncos(%.4f)=%18.16lf\nsin(%.4f)=%18.16lf\n",deg,cos,deg,sin);
_Rad2SinCos(rad,&sin,&cos);
printf("_Rad2SinCos:\ncos(%.4f)=%18.16lf\nsin(%.4f)=%18.16lf\n",rad,cos,rad,sin);
//system("PAUSE");
getch();
return 0;
}_Degree2SinCos:
cos(30.0000)=0.8660254037844386
sin(30.0000)=0.5000000000000000
_Rad2SinCos:
cos(0.5236)=0.8660254037844387
sin(0.5236)=0.4999999999999999
唯一不爽的是上面代码生成的可执行文件，360竟然报毒，没法蛋定了，你懂的。
6# G-Spider
:lol
python 3.2 昨天释出，我从官网 下载了安装的时候，360还很积极的拦截呢，说是
发现木马：Win32/Trojan.e8c
详细描述：
木马名称：Win32/Trojan.e8c
所在路径：G:\Documents and Settings\yoyo\桌面\stardict\python-3.2.msi 7# wayne
“永久免费”背后的生财之道。。。你懂的。。。:*-^ 上面的fsincos与math.h 中的sin(),cos()时间测试。
j=10~30N=2^j，对于每次的j++ ，计算N0=N/4个sin，cos值的时间，随j增加，时间减少接近于40%。Elapsed time:
_Rad2SinCos Std_SinCos (seconds)
0.00000000 0.00000000
0.00000000 0.00000000
0.00000000 0.00000000
0.00000000 0.00000000
0.00000000 0.00000000
0.00000000 0.00000000
0.01600000 0.01500000
0.01600000 0.01600000
0.03100000 0.04700000
0.06200000 0.07800000
0.12500000 0.18800000
0.23400000 0.37500000
0.43800000 0.73400000
0.87500000 1.46900000
1.81200000 3.00000000
3.86000000 6.25000000
7.04700000 11.84300000
14.06300000 23.81200000
28.17200000 47.43800000
56.42200000 94.78100000
请按任意键继续. . .测试代码：#include <stdio.h>
#include <math.h>
#include <time.h>
#define PI 3.1415926535897932
#define _2PI(2.0*PI)
#ifdef __cplusplus
#include <windows.h>
extern "C" {
#endif
void _stdcall _Degree2SinCos(float degree,double *sin,double *cos);
void _stdcall _Rad2SinCos(double Rad,double *sin,double *cos);
#ifdef __cplusplus
}
#endif
int main()
{
double cos0,sin0,cos1,sin1;
double t0,t1;
double i,N,N0,N1,theta;//1/N
int j;
printf("Elapsed time: \n_Rad2SinCos Std_SinCos (seconds)\n");
for(j=10;j<30;j++)
{
N=(unsigned int)1<<j;
N0=N/4.0;
N1=1.0/N;
//=======================================================
t0=clock();
for(i=0;i<=N0;i++)
{
theta=(double)i*N1*_2PI;
_Rad2SinCos(theta,&sin0,&cos0);
//printf("%18.16lf %18.16lf %18.16lf\n",theta,sin0,cos0);
}
printf("%.8f ",(clock()-t0)/CLOCKS_PER_SEC);
//=======================================================
t1=clock();
for(i=0;i<=N0;i++)
{
theta=(double)i*N1*_2PI;
sin1=sin(theta);
cos1=cos(theta);
//printf("%18.16lf %18.16lf %18.16lf\n",theta,sin1,cos1);
}
printf("%.8f \n",(clock()-t1)/CLOCKS_PER_SEC);
//=======================================================
}
system("PAUSE");
return 0;
} 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)得出，不知精度如何，有时间试试。
页:
[1]
2