#define pi 3.141593

#define f 3 //原正弦波频率

#define cx pi/3 //原正弦波初相，每次采样的时候是随机的

double out[10000]; //存插值后各点值用的，大点好

double in[10000]; //采样值缓存，大点好

//使用Kaiser窗截断，并使用参数9

//num:采样点数

//insert_num:每2个采样点插入的点数

//T:周期数

//start_phase:生成的正弦波的起始相位

//a:生成正弦波的a系数

//b:生成正弦波的b系数

void test_sinc(int sample_num, int insert_num, double T, double start_phase, double a, double b)

{

insert_num ++;

for(int m = 0;m < sample_num;m++){//采样

in[m] = a*sin( (2*pi*m*T)/(sample_num) + start_phase) + b;

}

for(int u = 0, i =0;u< sample_num;u++){

out[i++] = in[u];

for(int k = 1;k < insert_num;k++){

double sum = 0;

double t = (double)u + (double)k/((double)insert_num);

for(int m = 0;m < sample_num;m++){

double tmp = pi*(t -(double)m);

sum += in[m]*sin(tmp)/tmp;

}

out[i++] = sum;

}

}

}

//生成hn系数数组的算法，len是数组总长度，insert_num是插入点数

void SincHn(double *data, int len, int insert_num)

{

int end = len /2;

const int start = -end;

if(len%2 != 0)

end++;

insert_num ++;

for(int i = start; i < end; i++){

if(i == 0){

*data++ = 1;

continue;

}

double tmp = pi * (double)i/ (double)insert_num;

*data++ = sin(tmp)/tmp;

}

}

unsigned long Fact(unsigned long k)

{

unsigned long sum = 1;

for(unsigned long i = 1; i <= k; i++)

sum *= i;

return sum;

}

double besseli_0(double x, int len)

{

double sum = 0;

for(int i = 0; i < len; i++){

sum += std::pow(-1, i) * std::pow(x/2, 2*i)/(Fact((unsigned long)i) * Fact((unsigned long)i+1));

}

return sum;

}

void Kaiser(double *data, int N, double beta)

{

double theta;

for( int n = 0; n < N; n++){

theta = beta * std::sqrt( 1 - std::pow( ( (2 * (double)n/(N -1)) - 1),2 ) );

double res = besseli_0(theta, 10) / besseli_0(beta, 10);

if(res < 0)

res = 0;

*data++ = res;

}

}

一键复制

编辑

Web IDE

原始数据

按行查看

历史