jamie@141: /* jamie@141: * Copyright (C) 2012 Jamie Bullock jamie@140: * jamie@141: * Permission is hereby granted, free of charge, to any person obtaining a copy jamie@141: * of this software and associated documentation files (the "Software"), to jamie@141: * deal in the Software without restriction, including without limitation the jamie@141: * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or jamie@141: * sell copies of the Software, and to permit persons to whom the Software is jamie@141: * furnished to do so, subject to the following conditions: jamie@107: * jamie@141: * The above copyright notice and this permission notice shall be included in jamie@141: * all copies or substantial portions of the Software. jamie@107: * jamie@141: * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR jamie@141: * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, jamie@141: * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE jamie@141: * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER jamie@141: * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING jamie@141: * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS jamie@141: * IN THE SOFTWARE. jamie@107: * 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 jamie@107: jamie@154: #ifndef M_PI jamie@154: #define M_PI 3.14159265358979323846264338327 jamie@154: #endif jamie@154: jamie@107: #include "xtract_window_private.h" jamie@107: jamie@146: void gauss(double *window, const int N, const double sd) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1; jamie@146: double num, jamie@107: den, jamie@107: exponent; jamie@107: jamie@140: for (n = 0; n < N; n++) jamie@140: { jamie@107: jamie@146: num = n - M / 2.0; jamie@146: den = sd * M / 2.0; jamie@140: jamie@146: exponent = -0.5 * pow(num / den, 2); jamie@107: jamie@107: window[n] = exp(exponent); jamie@107: jamie@107: } jamie@107: } jamie@107: jamie@146: void hamming(double *window, const int N) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1; jamie@107: jamie@107: for (n = 0; n < N; n++) jamie@150: window[n] = 0.53836 - (0.46164 * cos(2.0 * M_PI * (double)n / M)); jamie@107: jamie@107: } jamie@107: jamie@146: void hann(double *window, const int N) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1; jamie@107: jamie@107: for (n = 0; n < N; n++) jamie@150: window[n] = 0.5 * (1.0 - cos(2.0 * M_PI * (double)n / M)); jamie@107: jamie@107: } jamie@107: jamie@146: void bartlett(double *window, const int N) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1; jamie@107: jamie@107: for (n = 0; n < N; n++) jamie@146: window[n] = 2.0 / M * (M / 2.0 - fabs(n - M / 2.0)); jamie@107: jamie@107: } jamie@107: jamie@146: void triangular(double *window, const int N) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1; jamie@107: jamie@107: for (n = 0; n < N; n++) jamie@146: window[n] = 2.0 / N * (N / 2.0 - fabs(n - M / 2.0)); jamie@107: } jamie@107: jamie@146: void bartlett_hann(double *window, const int N) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1, jamie@140: a0 = 0.62, jamie@140: a1 = 0.5, jamie@140: a2 = 0.38; jamie@146: double term1 = 0.0, jamie@146: term2 = 0.0; jamie@107: jamie@140: for (n = 0; n < N; n++) jamie@140: { jamie@107: jamie@146: term1 = a1 * fabs(n / M - 0.5); jamie@150: term2 = a2 * cos(2.0 * M_PI * (double)n / M); jamie@107: jamie@107: window[n] = a0 - term1 - term2; jamie@107: } jamie@107: } jamie@107: jamie@146: void blackman(double *window, const int N) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1, jamie@140: a0 = 0.42, jamie@140: a1 = 0.5, jamie@140: a2 = 0.08; jamie@146: double term1 = 0.0, jamie@146: term2 = 0.0; jamie@107: jamie@140: for (n = 0; n < N; n++) jamie@140: { jamie@140: jamie@150: term1 = a1 * cos(2.0 * M_PI * (double)n / M); jamie@150: term2 = a2 * cos(4.0 * M_PI * (double)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@146: double besselI0(double x) jamie@140: { jamie@107: jamie@146: double temp; jamie@146: double sum = 1.0; jamie@146: double u = 1.0; jamie@146: double halfx = x/2.0; jamie@140: int n = 1; jamie@107: jamie@140: do jamie@140: { jamie@107: jamie@146: temp = halfx/(double)n; jamie@140: u *=temp * temp; jamie@140: sum += u; jamie@140: n++; jamie@107: jamie@140: } jamie@140: while (u >= BIZ_EPSILON * sum); jamie@107: jamie@140: return(sum); jamie@107: jamie@107: } jamie@107: jamie@146: void kaiser(double *window, const int N, const double alpha) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1; jamie@146: double num; jamie@107: jamie@140: for (n = 0; n < N; n++) jamie@140: { jamie@107: jamie@146: num = besselI0(alpha * sqrt(1.0 - pow((2.0 * n / M - 1), 2))); jamie@107: window[n] = num / besselI0(alpha); jamie@140: jamie@107: } jamie@107: } jamie@107: jamie@146: void blackman_harris(double *window, const int N) jamie@140: { jamie@107: jamie@107: int n; jamie@146: const double M = N - 1, jamie@140: a0 = 0.35875, jamie@140: a1 = 0.48829, jamie@140: a2 = 0.14128, jamie@140: a3 = 0.01168; jamie@146: double term1 = 0.0, jamie@146: term2 = 0.0, jamie@146: term3 = 0.0; jamie@107: jamie@140: for (n = 0; n < N; n++) jamie@140: { jamie@107: jamie@150: term1 = a1 * cos(2.0 * M_PI * n / M); jamie@150: term2 = a2 * cos(4.0 * M_PI * n / M); jamie@150: term3 = a3 * cos(6.0 * M_PI * n / M); jamie@107: jamie@107: window[n] = a0 - term1 + term2 - term3; jamie@107: } jamie@107: }