jamie@107: /* libxtract feature extraction library
jamie@107:  *  
jamie@107:  * Copyright (C) 2006 Jamie Bullock
jamie@107:  *
jamie@107:  * This program is free software; you can redistribute it and/or modify
jamie@107:  * it under the terms of the GNU General Public License as published by
jamie@107:  * the Free Software Foundation; either version 2 of the License, or
jamie@107:  * (at your option) any later version.
jamie@107:  *
jamie@107:  * This program is distributed in the hope that it will be useful,
jamie@107:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
jamie@107:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
jamie@107:  * GNU General Public License for more details.
jamie@107:  *
jamie@107:  * You should have received a copy of the GNU General Public License
jamie@107:  * along with this program; if not, write to the Free Software
jamie@107:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, 
jamie@107:  * USA.
jamie@107:  */
jamie@107: 
jamie@107: /* window.c: defines window generation functions (formulae courtesy of Wikipedia (http://en.wikipedia.org/wiki/Window_function) */
jamie@107: 
jamie@107: #include <math.h>
jamie@107: 
jamie@107: #include "xtract_window_private.h"
jamie@107: 
jamie@107: void gauss(float *window, const int N, const float sd){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1;
jamie@107:     float num,
jamie@107:           den,
jamie@107:           exponent;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++) {
jamie@107: 
jamie@107:         num = n - M / 2.f;
jamie@107:         den = sd * M / 2.f;
jamie@107:         
jamie@107:         exponent = -0.5 * powf(num / den, 2);
jamie@107: 
jamie@107:         window[n] = exp(exponent);
jamie@107: 
jamie@107:     }
jamie@107: }
jamie@107: 
jamie@107: void hamming(float *window, const int N){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++)
jamie@107:         window[n] = 0.53836 - (0.46164 * cosf(2.0 * PI * (float)n / M));
jamie@107: 
jamie@107: }
jamie@107: 
jamie@107: void hann(float *window, const int N){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++)
jamie@107:         window[n] = 0.5 * (1.0 - cosf(2.0 * PI * (float)n / M));
jamie@107: 
jamie@107: }
jamie@107: 
jamie@107: void bartlett(float *window, const int N){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++)
jamie@107:         window[n] = 2.f / M * (M / 2.f - fabsf(n - M / 2.f));
jamie@107: 
jamie@107: }
jamie@107: 
jamie@107: void triangular(float *window, const int N){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++)
jamie@107:         window[n] = 2.f / N * (N / 2.f - fabsf(n - M / 2.f));
jamie@107: }
jamie@107: 
jamie@107: void bartlett_hann(float *window, const int N){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1,
jamie@107:           a0 = 0.62,
jamie@107:           a1 = 0.5,
jamie@107:           a2 = 0.38;
jamie@107:     float term1 = 0.f,
jamie@107:           term2 = 0.f;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++){
jamie@107: 
jamie@107:         term1 = a1 * fabsf(n / M - 0.5);
jamie@107:         term2 = a2 * cosf(2.0 * PI * (float)n / M);
jamie@107: 
jamie@107:         window[n] = a0 - term1 - term2;
jamie@107:     }
jamie@107: }
jamie@107: 
jamie@107: void blackman(float *window, const int N){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1,
jamie@107:           a0 = 0.42,
jamie@107:           a1 = 0.5,
jamie@107:           a2 = 0.08;
jamie@107:     float term1 = 0.f,
jamie@107:           term2 = 0.f;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++) {
jamie@107:     
jamie@107:         term1 = a1 * cosf(2.0 * PI * (float)n / M);
jamie@107:         term2 = a2 * cosf(4.0 * PI * (float)n / M);
jamie@107: 
jamie@107:         window[n] = a0 - term1 + term2;
jamie@107:     }
jamie@107: }
jamie@107: 
jamie@107: #define BIZ_EPSILON 1E-21 // Max error acceptable 
jamie@107: 
jamie@107: /* Based on code from mplayer window.c, and somewhat beyond me */
jamie@107: float besselI0(float x){
jamie@107: 
jamie@107:   float temp;
jamie@107:   float sum   = 1.0;
jamie@107:   float u     = 1.0;
jamie@107:   float halfx = x/2.0;
jamie@107:   int      n     = 1;
jamie@107: 
jamie@107:   do {
jamie@107: 
jamie@107:     temp = halfx/(float)n;
jamie@107:     u *=temp * temp;
jamie@107:     sum += u;
jamie@107:     n++;
jamie@107: 
jamie@107:   } while (u >= BIZ_EPSILON * sum);
jamie@107: 
jamie@107:   return(sum);
jamie@107: 
jamie@107: }
jamie@107: 
jamie@107: void kaiser(float *window, const int N, const float alpha){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1;
jamie@107:     float num;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++) {
jamie@107: 
jamie@107:         num = besselI0(alpha * sqrtf(1.0 - powf((2.0 * n / M - 1), 2)));
jamie@107:         window[n] = num / besselI0(alpha);
jamie@107:         
jamie@107:     }
jamie@107: }
jamie@107: 
jamie@107: void blackman_harris(float *window, const int N){
jamie@107: 
jamie@107:     int n;
jamie@107:     const float M = N - 1,
jamie@107:           a0 = 0.35875,
jamie@107:           a1 = 0.48829,
jamie@107:           a2 = 0.14128,
jamie@107:           a3 = 0.01168;
jamie@107:     float term1 = 0.f,
jamie@107:           term2 = 0.f,
jamie@107:           term3 = 0.f;
jamie@107: 
jamie@107:     for (n = 0; n < N; n++) {
jamie@107: 
jamie@107:         term1 = a1 * cosf(2.0 * PI * n / M);
jamie@107:         term2 = a2 * cosf(4.0 * PI * n / M);
jamie@107:         term3 = a3 * cosf(6.0 * PI * n / M);
jamie@107: 
jamie@107:         window[n] = a0 - term1 + term2 - term3;
jamie@107:     }
jamie@107: }