00001
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "dsp.h"
00047
00048 #if defined(FORCE_ALL_GENERICS) || \
00049 defined(FORCE_GENERIC_WIN32_KAISER) || \
00050 !defined(TARGET_SPECIFIC_WIN32_KAISER)
00051
00052
00053 static dsp32_t dsp32_op_kaiser_i0(S64 ax, int *power)
00054 {
00055 S64 ans, y, k, num_temp;
00056
00057 const S64 A0 = (S64) (1.0 * (1 << 28));
00058 const S64 A1 = (S64) ((3.5156229 * (1 << 28)) / (3.75*3.75));
00059 const S64 A2 = (S64) ((3.0899424 * (1 << 28)) / (3.75*3.75*3.75*3.75));
00060 const S64 A3 = (S64) ((1.2067492 * (1 << 28)) / (3.75*3.75*3.75*3.75*3.75*3.75));
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));
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));
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
00077 *power = 3;
00078
00079
00080
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
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
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
00131
00132 k = 0;
00133 num_temp = ax >> DSP32_QB;
00134 while(num_temp)
00135 {
00136 num_temp >>= 2;
00137 k++;
00138 }
00139
00140 *power -= k;
00141
00142
00143 num_temp = ax >> (k << 1);
00144
00145
00146 num_temp = dsp32_op_sqrt(num_temp);
00147
00148
00149 ans = (y*ans) / num_temp;
00150 }
00151
00152 return ans;
00153 }
00154
00155
00156
00157
00158
00159
00160
00161 void dsp32_win_kaiser(dsp32_t *vect1, int size, int alpha)
00162 {
00163 int n, power_num, power_den;
00164 S64 pi_alpha, pi_alpha_div;
00165 S64 temp64;
00166 dsp32_t s, t, temp32, num, den;
00167
00168
00169 s = (DSP32_Q(1.) / (size - 1)) << 1;
00170
00171 pi_alpha = (S64) (CST_PI*(1LL << DSP32_QB))*alpha;
00172
00173 pi_alpha_div = pi_alpha >> (DSP32_QB >> 1);
00174
00175 den = dsp32_op_kaiser_i0(pi_alpha, &power_den);
00176
00177 t = 0;
00178
00179 for(n=0; n<(size >> 1); n++)
00180 {
00181
00182 temp32 = t - DSP32_Q(1);
00183
00184 temp32 = (((S64) temp32) * ((S64) temp32)) >> DSP32_QB;
00185
00186 temp32 = dsp32_op_sqrt(DSP32_Q(1.) - temp32);
00187
00188
00189 #if (DSP32_QB & 1)
00190 temp64 = (pi_alpha_div * ((S64) temp32)) >> ((DSP32_QB >> 1) + 1);
00191 #else
00192 temp64 = (pi_alpha_div * ((S64) temp32)) >> (DSP32_QB >> 1);
00193 #endif
00194
00195 num = dsp32_op_kaiser_i0(temp64, &power_num);
00196
00197
00198 power_num = power_den - power_num;
00199
00200 temp64 = ((S64) num) << (DSP32_QB);
00201 temp64 = temp64 / (S64) den;
00202 temp64 >>= power_num;
00203
00204 vect1[n] = temp64;
00205 vect1[size-n-1] = temp64;
00206
00207
00208 t += s;
00209 }
00210
00211
00212 if (size & 1)
00213 vect1[size >> 1] = DSP32_Q(1.);
00214 }
00215
00216 #endif