Alexandre Lision | 8af73cb | 2013-12-10 14:11:20 -0500 | [diff] [blame] | 1 | /* Copyright (C) 2006 Jean-Marc Valin */ |
| 2 | /** |
| 3 | @file filterbank.c |
| 4 | @brief Converting between psd and filterbank |
| 5 | */ |
| 6 | /* |
| 7 | Redistribution and use in source and binary forms, with or without |
| 8 | modification, are permitted provided that the following conditions are |
| 9 | met: |
| 10 | |
| 11 | 1. Redistributions of source code must retain the above copyright notice, |
| 12 | this list of conditions and the following disclaimer. |
| 13 | |
| 14 | 2. Redistributions in binary form must reproduce the above copyright |
| 15 | notice, this list of conditions and the following disclaimer in the |
| 16 | documentation and/or other materials provided with the distribution. |
| 17 | |
| 18 | 3. The name of the author may not be used to endorse or promote products |
| 19 | derived from this software without specific prior written permission. |
| 20 | |
| 21 | THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
| 22 | IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
| 23 | OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE |
| 24 | DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, |
| 25 | INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES |
| 26 | (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR |
| 27 | SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
| 28 | HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, |
| 29 | STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN |
| 30 | ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| 31 | POSSIBILITY OF SUCH DAMAGE. |
| 32 | */ |
| 33 | |
| 34 | #ifdef HAVE_CONFIG_H |
| 35 | #include "config.h" |
| 36 | #endif |
| 37 | |
| 38 | #include "filterbank.h" |
| 39 | #include "arch.h" |
| 40 | #include <math.h> |
| 41 | #include "math_approx.h" |
| 42 | #include "os_support.h" |
| 43 | |
| 44 | #ifdef FIXED_POINT |
| 45 | |
| 46 | #define toBARK(n) (MULT16_16(26829,spx_atan(SHR32(MULT16_16(97,n),2))) + MULT16_16(4588,spx_atan(MULT16_32_Q15(20,MULT16_16(n,n)))) + MULT16_16(3355,n)) |
| 47 | |
| 48 | #else |
| 49 | #define toBARK(n) (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n)) |
| 50 | #endif |
| 51 | |
| 52 | #define toMEL(n) (2595.f*log10(1.f+(n)/700.f)) |
| 53 | |
| 54 | FilterBank *filterbank_new(int banks, spx_word32_t sampling, int len, int type) |
| 55 | { |
| 56 | FilterBank *bank; |
| 57 | spx_word32_t df; |
| 58 | spx_word32_t max_mel, mel_interval; |
| 59 | int i; |
| 60 | int id1; |
| 61 | int id2; |
| 62 | df = DIV32(SHL32(sampling,15),MULT16_16(2,len)); |
| 63 | max_mel = toBARK(EXTRACT16(sampling/2)); |
| 64 | mel_interval = PDIV32(max_mel,banks-1); |
| 65 | |
| 66 | bank = (FilterBank*)speex_alloc(sizeof(FilterBank)); |
| 67 | bank->nb_banks = banks; |
| 68 | bank->len = len; |
| 69 | bank->bank_left = (int*)speex_alloc(len*sizeof(int)); |
| 70 | bank->bank_right = (int*)speex_alloc(len*sizeof(int)); |
| 71 | bank->filter_left = (spx_word16_t*)speex_alloc(len*sizeof(spx_word16_t)); |
| 72 | bank->filter_right = (spx_word16_t*)speex_alloc(len*sizeof(spx_word16_t)); |
| 73 | /* Think I can safely disable normalisation that for fixed-point (and probably float as well) */ |
| 74 | #ifndef FIXED_POINT |
| 75 | bank->scaling = (float*)speex_alloc(banks*sizeof(float)); |
| 76 | #endif |
| 77 | for (i=0;i<len;i++) |
| 78 | { |
| 79 | spx_word16_t curr_freq; |
| 80 | spx_word32_t mel; |
| 81 | spx_word16_t val; |
| 82 | curr_freq = EXTRACT16(MULT16_32_P15(i,df)); |
| 83 | mel = toBARK(curr_freq); |
| 84 | if (mel > max_mel) |
| 85 | break; |
| 86 | #ifdef FIXED_POINT |
| 87 | id1 = DIV32(mel,mel_interval); |
| 88 | #else |
| 89 | id1 = (int)(floor(mel/mel_interval)); |
| 90 | #endif |
| 91 | if (id1>banks-2) |
| 92 | { |
| 93 | id1 = banks-2; |
| 94 | val = Q15_ONE; |
| 95 | } else { |
| 96 | val = DIV32_16(mel - id1*mel_interval,EXTRACT16(PSHR32(mel_interval,15))); |
| 97 | } |
| 98 | id2 = id1+1; |
| 99 | bank->bank_left[i] = id1; |
| 100 | bank->filter_left[i] = SUB16(Q15_ONE,val); |
| 101 | bank->bank_right[i] = id2; |
| 102 | bank->filter_right[i] = val; |
| 103 | } |
| 104 | |
| 105 | /* Think I can safely disable normalisation for fixed-point (and probably float as well) */ |
| 106 | #ifndef FIXED_POINT |
| 107 | for (i=0;i<bank->nb_banks;i++) |
| 108 | bank->scaling[i] = 0; |
| 109 | for (i=0;i<bank->len;i++) |
| 110 | { |
| 111 | int id = bank->bank_left[i]; |
| 112 | bank->scaling[id] += bank->filter_left[i]; |
| 113 | id = bank->bank_right[i]; |
| 114 | bank->scaling[id] += bank->filter_right[i]; |
| 115 | } |
| 116 | for (i=0;i<bank->nb_banks;i++) |
| 117 | bank->scaling[i] = Q15_ONE/(bank->scaling[i]); |
| 118 | #endif |
| 119 | return bank; |
| 120 | } |
| 121 | |
| 122 | void filterbank_destroy(FilterBank *bank) |
| 123 | { |
| 124 | speex_free(bank->bank_left); |
| 125 | speex_free(bank->bank_right); |
| 126 | speex_free(bank->filter_left); |
| 127 | speex_free(bank->filter_right); |
| 128 | #ifndef FIXED_POINT |
| 129 | speex_free(bank->scaling); |
| 130 | #endif |
| 131 | speex_free(bank); |
| 132 | } |
| 133 | |
| 134 | void filterbank_compute_bank32(FilterBank *bank, spx_word32_t *ps, spx_word32_t *mel) |
| 135 | { |
| 136 | int i; |
| 137 | for (i=0;i<bank->nb_banks;i++) |
| 138 | mel[i] = 0; |
| 139 | |
| 140 | for (i=0;i<bank->len;i++) |
| 141 | { |
| 142 | int id; |
| 143 | id = bank->bank_left[i]; |
| 144 | mel[id] += MULT16_32_P15(bank->filter_left[i],ps[i]); |
| 145 | id = bank->bank_right[i]; |
| 146 | mel[id] += MULT16_32_P15(bank->filter_right[i],ps[i]); |
| 147 | } |
| 148 | /* Think I can safely disable normalisation that for fixed-point (and probably float as well) */ |
| 149 | #ifndef FIXED_POINT |
| 150 | /*for (i=0;i<bank->nb_banks;i++) |
| 151 | mel[i] = MULT16_32_P15(Q15(bank->scaling[i]),mel[i]); |
| 152 | */ |
| 153 | #endif |
| 154 | } |
| 155 | |
| 156 | void filterbank_compute_psd16(FilterBank *bank, spx_word16_t *mel, spx_word16_t *ps) |
| 157 | { |
| 158 | int i; |
| 159 | for (i=0;i<bank->len;i++) |
| 160 | { |
| 161 | spx_word32_t tmp; |
| 162 | int id1, id2; |
| 163 | id1 = bank->bank_left[i]; |
| 164 | id2 = bank->bank_right[i]; |
| 165 | tmp = MULT16_16(mel[id1],bank->filter_left[i]); |
| 166 | tmp += MULT16_16(mel[id2],bank->filter_right[i]); |
| 167 | ps[i] = EXTRACT16(PSHR32(tmp,15)); |
| 168 | } |
| 169 | } |
| 170 | |
| 171 | |
| 172 | #ifndef FIXED_POINT |
| 173 | void filterbank_compute_bank(FilterBank *bank, float *ps, float *mel) |
| 174 | { |
| 175 | int i; |
| 176 | for (i=0;i<bank->nb_banks;i++) |
| 177 | mel[i] = 0; |
| 178 | |
| 179 | for (i=0;i<bank->len;i++) |
| 180 | { |
| 181 | int id = bank->bank_left[i]; |
| 182 | mel[id] += bank->filter_left[i]*ps[i]; |
| 183 | id = bank->bank_right[i]; |
| 184 | mel[id] += bank->filter_right[i]*ps[i]; |
| 185 | } |
| 186 | for (i=0;i<bank->nb_banks;i++) |
| 187 | mel[i] *= bank->scaling[i]; |
| 188 | } |
| 189 | |
| 190 | void filterbank_compute_psd(FilterBank *bank, float *mel, float *ps) |
| 191 | { |
| 192 | int i; |
| 193 | for (i=0;i<bank->len;i++) |
| 194 | { |
| 195 | int id = bank->bank_left[i]; |
| 196 | ps[i] = mel[id]*bank->filter_left[i]; |
| 197 | id = bank->bank_right[i]; |
| 198 | ps[i] += mel[id]*bank->filter_right[i]; |
| 199 | } |
| 200 | } |
| 201 | |
| 202 | void filterbank_psy_smooth(FilterBank *bank, float *ps, float *mask) |
| 203 | { |
| 204 | /* Low freq slope: 14 dB/Bark*/ |
| 205 | /* High freq slope: 9 dB/Bark*/ |
| 206 | /* Noise vs tone: 5 dB difference */ |
| 207 | /* FIXME: Temporary kludge */ |
| 208 | float bark[100]; |
| 209 | int i; |
| 210 | /* Assumes 1/3 Bark resolution */ |
| 211 | float decay_low = 0.34145f; |
| 212 | float decay_high = 0.50119f; |
| 213 | filterbank_compute_bank(bank, ps, bark); |
| 214 | for (i=1;i<bank->nb_banks;i++) |
| 215 | { |
| 216 | /*float decay_high = 13-1.6*log10(bark[i-1]); |
| 217 | decay_high = pow(10,(-decay_high/30.f));*/ |
| 218 | bark[i] = bark[i] + decay_high*bark[i-1]; |
| 219 | } |
| 220 | for (i=bank->nb_banks-2;i>=0;i--) |
| 221 | { |
| 222 | bark[i] = bark[i] + decay_low*bark[i+1]; |
| 223 | } |
| 224 | filterbank_compute_psd(bank, bark, mask); |
| 225 | } |
| 226 | |
| 227 | #endif |