windows.h

Go to the documentation of this file.

/*
        dsp/windows.h
        
        Copyright 2004-11 Tim Goetze <tim@quitte.de>
        
        http://quitte.de/dsp/
 
        select few common windowing algorithms.
 
*/
/*
        This program is free software; you can redistribute it and/or
        modify it under the terms of the GNU General Public License
        as published by the Free Software Foundation; either version 2
        of the License, or (at your option) any later version.
 
        This program is distributed in the hope that it will be useful,
        but WITHOUT ANY WARRANTY; without even the implied warranty of
        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
        GNU General Public License for more details.
 
        You should have received a copy of the GNU General Public License
        along with this program; if not, write to the Free Software
        Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
        02111-1307, USA or point your web browser to http://www.gnu.org.
*/
 
#ifndef _DSP_WINDOWS_H_ 
#define _DSP_WINDOWS_H_
 
namespace DSP {
        
/* prototypes for window value application ... */
typedef void (*window_sample_func_t) (sample_t &, sample_t);
 
/* ... which go as template parameters for the window calculation below */
inline void
store_sample (sample_t & d, sample_t s)
{
        d = s;
}
 
inline void
apply_window (sample_t &d, sample_t s)
{
        d *= s;
}
 
template <window_sample_func_t F>
void
hanning (sample_t * s, int n)
{
        /* could speed up by using DSP::Sine but we rarely use this window */
        for (int i = 0; i < n; ++i)
        {
                double f = (double) i / n - 1;
                F (s[i], .5 - .5 * cos (2 * M_PI * f));
        }
}
 
template <window_sample_func_t F>
void
blackman (sample_t * s, int n)
{
        float w = n;
 
        for (int i = 0; i < n; ++i)
        {
                float f = (float) i;
 
                double b = .42f - 
                                                .5f * cos (2.f * f * M_PI / w) + 
                                                .08 * cos (4.f * f * M_PI / w);
 
                F (s[i], b);
        }
}
 
template <window_sample_func_t F>
void
blackman_harris (sample_t * s, int n)
{
        double w1 = 2.f * M_PI / (n - 1);
        double w2 = 2.f * w1;
        double w3 = 3.f * w1;
 
        for (int i = 0; i < n; ++i)
        {
                double f = (double) i;
 
                double bh = .35875f - 
                                .48829f * cos (w1 * f) + 
                                .14128f * cos (w2 * f) - 
                                .01168f * cos (w3 * f);
 
                bh *= .761f;
 
                F (s[i], bh);
        }
}
 
/* helper for the kaiser window, courtesy of R. Dobson, courtesy of csound */
inline double 
besseli (double x)
{
        double ax, ans;
        double y;
 
        if ((ax = fabs (x)) < 3.75)     
        {
                y = x / 3.75;
                y *= y;
                ans = (1.0 + y * (3.5156229 +
                        y * (3.0899424 +
                                y * (1.2067492 +
                                        y * (0.2659732 +
                                                y * (0.360768e-1 +
                                                        y * 0.45813e-2))))));
        }
        else 
        {
                y = 3.75 / ax;
                ans = ((exp (ax) / sqrt (ax))
                        * (0.39894228 +
                                y * (0.1328592e-1 +
                                        y * (0.225319e-2 +
                                                y * (-0.157565e-2 +
                                                        y * (0.916281e-2 +
                                                                y * (-0.2057706e-1 +
                                                                        y * (0.2635537e-1 +
                                                                                y * (-0.1647633e-1 +
                                                                                        y * 0.392377e-2)))))))));
        }
 
        return ans;
}
 
template <window_sample_func_t F>
void
kaiser (sample_t * s, int n, double beta)
{
        double bb = besseli (beta);
        int si = 0;
 
        for (double i = -n / 2 + .1; si < n; ++si, ++i)
        {
                double k = besseli ((beta * sqrt (1 - pow ((2 * i / (n - 1)), 2)))) / bb;
 
                /* can you spell hack */
                if (!isfinite (k) || isnan(k))
                        k = 0;
 
                F (s[si], k);
        }
        /* asymmetrical hack: sort out first value!
        win[0] = win[len-1];
        */
}
 
}; /* namespace DSP */
 
#endif /* _DSP_WINDOWS_H_ */