|
|
Попробуйте вот так вот, здесь промежут. результаты 64-х разрядные, т.е. в линии задержки. Комментировать не буду, думаю разберетесь.
ytf1,ytf2,ypf - беззнаковые на самом деле. Добавок на расширение разрядности отмечен комментарием /***fractional part****/
#include "stdafx.h"
#include <math.h>
#include <stdio.h>
/**************************/
int B0=1072470442;
int B1=-2144940883;
int B2=1072470442;
int A1=-2144939378;
int A2=1071200564;
/**************************/double PI=3.14159265358979323846;
double FDISKR=1500.;
double F0=0.5;
double AM0=1024.*1024.*1024.;
double AM1=16384.;
/***************************/
int XIN[65536];
void main(void)
{
FILE *filout;
__int64 yp,xs,xst1,xst2,yt1,yt2;
__int64 ytf1,ytf2,ypf;
double x,y;
int i,j,k,l,yout;
for(i=0;i<65536;i++)
{
j=i%3000;
k=j/750;
l=j%750;
x=2.*PI*F0*((double)l)/FDISKR;
if(k&1) x=(PI/2)-x;
y=sin(x);
if(k>=2) y=-y;
if(i<32768) y*=AM0;
else y*=AM1;
XIN[i]=(int)(floor(y+0.5));
}
yt1=0;
yt2=0;
ytf1=0;
ytf2=0;
xst1=0;
xst2=0;
filout=fopen("filtout.prn","wb");
x=0.;
y=0.;
for(i=0;i<65536;i++)
{
xs=(__int64)XIN[i];
yp=xs*(__int64)B0;
yp+=(xst1*(__int64)B1);
yp+=(xst2*(__int64)B2);
yp-=(yt1*(__int64)A1);
yp-=(yt2*(__int64)A2);
/****fractional part ***/
ypf=ytf1*(__int64)A1;
ypf+=ytf2*(__int64)A2;
ypf+=0x20000000UI64;
ypf>>=30;
yp-=ypf;
ytf2=ytf1;
ytf1=yp&0x3fffffffUI64;
/*************************/
xst2=xst1;
xst1=xs;
yt2=yt1;
yt1=yp>>30;
yout=(int)((yp+0x20000000UI64)>>30); //filter output
fprintf(filout,"%d\r\n",yout);
if(i<32768&&i>=29768) x+=((double)yout)*((double)yout);
if(i>=62536) y+=((double)yout)*((double)yout);
}
fclose(filout);
x=sqrt((x/3000.));
y=sqrt((y/3000.));
if(y!=0.) printf("HIGH LEVEL=%lf LOW LEVEL=%lf RELATION=%lf\n",x,y,(x/y));
else printf("HIGH LEVEL=%lf, LOW LEVEL=0\n",x);}
|
|
