Definition in file win_dsp32_kaiser.c.
#include "dsp.h"
Go to the source code of this file.
Functions | |
static dsp32_t | dsp32_op_kaiser_i0 (S64 ax, int *power) |
void | dsp32_win_kaiser (dsp32_t *vect1, int size, int alpha) |
32-bit fixed point version of the kaiser windowing function. |
static dsp32_t dsp32_op_kaiser_i0 | ( | S64 | ax, | |
int * | power | |||
) | [static] |
Definition at line 53 of file win_dsp32_kaiser.c.
References dsp32_op_exp(), dsp32_op_sqrt(), DSP32_Q, and DSP32_QB.
Referenced by dsp32_win_kaiser().
00054 { 00055 S64 ans, y, k, num_temp; 00056 00057 const S64 A0 = (S64) (1.0 * (1 << 28)); // 1.0 << 12 00058 const S64 A1 = (S64) ((3.5156229 * (1 << 28)) / (3.75*3.75)); // 14.0625 0.2499999 << 28 00059 const S64 A2 = (S64) ((3.0899424 * (1 << 28)) / (3.75*3.75*3.75*3.75)); // 197.75391 0.0156252 << 28 00060 const S64 A3 = (S64) ((1.2067492 * (1 << 28)) / (3.75*3.75*3.75*3.75*3.75*3.75)); // 2780.9143 0.0004339 << 12 00061 const S64 A4 = (S64) ((0.2659732 * (1 << 28)) / (3.75*3.75*3.75*3.75*3.75*3.75*3.75*3.75)); // 39106.607 0.0000068 << 12 00062 const S64 A5 = (S64) ((0.0360768 * (1 << 28)) / (3.75*3.75*3.75*3.75*3.75*3.75*3.75*3.75*3.75*3.75)); // 549936.67 6.560D-08 << 12 00063 00064 const S64 B0 = DSP32_Q(0.39894228); 00065 const S64 B1 = DSP32_Q(0.01328592); 00066 const S64 B2 = DSP32_Q(0.00225319); 00067 const S64 B3 = DSP32_Q(-0.00157565); 00068 const S64 B4 = DSP32_Q(0.00916281); 00069 const S64 B5 = DSP32_Q(-0.02057706); 00070 00071 const S64 cst_inv_ln_2 = DSP32_Q((1./0.69314718055994530941723212145818)/(1.*(1 << (DSP32_QB >> 1)))); 00072 const S64 cst_ln_2 = DSP32_Q(0.69314718055994530941723212145818); 00073 00074 if (ax < (S64) (0xFLL << (DSP32_QB - 2))) 00075 { 00076 // calculate the exponent* 00077 *power = 3; 00078 00079 // Format Q4.12 to contain 3.75^2 00080 // Change format 00081 #if (DSP32_QB > 28) 00082 y = ax >> (DSP32_QB - 28); 00083 #else 00084 y = ax << (28 - DSP32_QB); 00085 #endif 00086 y = (y*y) >> 28; 00087 00088 ans = A5; 00089 ans = (((ans*y) >> 28) + A4); 00090 ans = (((ans*y) >> 28) + A3); 00091 ans = (((ans*y) >> 28) + A2); 00092 ans = (((ans*y) >> 28) + A1); 00093 ans = (((ans*y) >> 28) + A0); 00094 00095 #if (DSP32_QB > 31) 00096 ans = ans << (DSP32_QB - 31); 00097 #else 00098 ans = ans >> (31 - DSP32_QB); 00099 #endif 00100 } 00101 else 00102 { 00103 // ans is in the range [1; 0] 00104 ans = ((S64) (0xFLL << (DSP32_QB - 2))) / ax; 00105 00106 y = B5; 00107 y = (((y*ans) >> DSP32_QB) + B4); 00108 y = (((y*ans) >> DSP32_QB) + B3); 00109 y = (((y*ans) >> DSP32_QB) + B2); 00110 y = (((y*ans) >> DSP32_QB) + B1); 00111 y = (((y*ans) >> DSP32_QB) + B0); 00112 00113 /***** Computation of exp *****/ 00114 #if (DSP32_QB & 1) 00115 ans = (ax >> ((DSP32_QB >> 1) + 1)); 00116 #else 00117 ans = (ax >> (DSP32_QB >> 1)); 00118 #endif 00119 00120 k = ((cst_inv_ln_2*ans) >> DSP32_QB) + 1; 00121 num_temp = ax - k*cst_ln_2; 00122 if (num_temp > 0) 00123 num_temp = ax - (++k)*cst_ln_2; 00124 00125 ans = dsp32_op_exp(num_temp); 00126 00127 *power = k; 00128 /******************************/ 00129 00130 /***** Computation of sqrt *****/ 00131 // get ~ int(log4(num)) 00132 k = 0; 00133 num_temp = ax >> DSP32_QB; 00134 while(num_temp) 00135 { 00136 num_temp >>= 2; 00137 k++; 00138 } 00139 // subtract k to the power because of the division 00140 *power -= k; 00141 00142 // ax = ax / 4^k 00143 num_temp = ax >> (k << 1); 00144 00145 // now ax is in the range [1; 0] 00146 num_temp = dsp32_op_sqrt(num_temp); 00147 /******************************/ 00148 00149 ans = (y*ans) / num_temp; 00150 } 00151 00152 return ans; 00153 }