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