1 /*
2 dsp/windows.h
3
4 Copyright 2004-11 Tim Goetze <tim@quitte.de>
5
6 http://quitte.de/dsp/
7
8 select few common windowing algorithms.
9
10 */
11 /*
12 This program is free software; you can redistribute it and/or
13 modify it under the terms of the GNU General Public License
14 as published by the Free Software Foundation; either version 2
15 of the License, or (at your option) any later version.
16
17 This program is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 GNU General Public License for more details.
21
22 You should have received a copy of the GNU General Public License
23 along with this program; if not, write to the Free Software
24 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
25 02111-1307, USA or point your web browser to http://www.gnu.org.
26 */
27
28 #ifndef _DSP_WINDOWS_H_
29 #define _DSP_WINDOWS_H_
30
31 namespace DSP {
32
33 /* prototypes for window value application ... */
34 typedef void (*window_sample_func_t) (sample_t &, sample_t);
35
36 /* ... which go as template parameters for the window calculation below */
37 inline void
store_sample(sample_t & d,sample_t s)38 store_sample (sample_t & d, sample_t s)
39 {
40 d = s;
41 }
42
43 inline void
apply_window(sample_t & d,sample_t s)44 apply_window (sample_t &d, sample_t s)
45 {
46 d *= s;
47 }
48
49 template <window_sample_func_t F>
50 void
hanning(sample_t * s,int n)51 hanning (sample_t * s, int n)
52 {
53 /* could speed up by using DSP::Sine but we rarely use this window */
54 for (int i = 0; i < n; ++i)
55 {
56 register double f = (double) i / n - 1;
57 F (s[i], .5 - .5 * cos (2 * M_PI * f));
58 }
59 }
60
61 template <window_sample_func_t F>
62 void
blackman(sample_t * s,int n)63 blackman (sample_t * s, int n)
64 {
65 register float w = n;
66
67 for (int i = 0; i < n; ++i)
68 {
69 register float f = (float) i;
70
71 register double b = .42f -
72 .5f * cos (2.f * f * M_PI / w) +
73 .08 * cos (4.f * f * M_PI / w);
74
75 F (s[i], b);
76 }
77 }
78
79 template <window_sample_func_t F>
80 void
blackman_harris(sample_t * s,int n)81 blackman_harris (sample_t * s, int n)
82 {
83 register double w1 = 2.f * M_PI / (n - 1);
84 register double w2 = 2.f * w1;
85 register double w3 = 3.f * w1;
86
87 for (int i = 0; i < n; ++i)
88 {
89 register double f = (double) i;
90
91 register double bh = .35875f -
92 .48829f * cos (w1 * f) +
93 .14128f * cos (w2 * f) -
94 .01168f * cos (w3 * f);
95
96 bh *= .761f;
97
98 F (s[i], bh);
99 }
100 }
101
102 /* helper for the kaiser window, courtesy of R. Dobson, courtesy of csound */
103 inline double
besseli(double x)104 besseli (double x)
105 {
106 double ax, ans;
107 double y;
108
109 if ((ax = fabs (x)) < 3.75)
110 {
111 y = x / 3.75;
112 y *= y;
113 ans = (1.0 + y * (3.5156229 +
114 y * (3.0899424 +
115 y * (1.2067492 +
116 y * (0.2659732 +
117 y * (0.360768e-1 +
118 y * 0.45813e-2))))));
119 }
120 else
121 {
122 y = 3.75 / ax;
123 ans = ((exp (ax) / sqrt (ax))
124 * (0.39894228 +
125 y * (0.1328592e-1 +
126 y * (0.225319e-2 +
127 y * (-0.157565e-2 +
128 y * (0.916281e-2 +
129 y * (-0.2057706e-1 +
130 y * (0.2635537e-1 +
131 y * (-0.1647633e-1 +
132 y * 0.392377e-2)))))))));
133 }
134
135 return ans;
136 }
137
138 template <window_sample_func_t F>
139 void
kaiser(sample_t * s,int n,double beta)140 kaiser (sample_t * s, int n, double beta)
141 {
142 double bb = besseli (beta);
143 int si = 0;
144
145 for (double i = -n / 2 + .1; si < n; ++si, ++i)
146 {
147 double k = besseli ((beta * sqrt (1 - pow ((2 * i / (n - 1)), 2)))) / bb;
148
149 /* can you spell hack */
150 if (!finite (k) || isnan(k))
151 k = 0;
152
153 F (s[si], k);
154 }
155 /* asymmetrical hack: sort out first value!
156 win[0] = win[len-1];
157 */
158 }
159
160 }; /* namespace DSP */
161
162 #endif /* _DSP_WINDOWS_H_ */
163