OXYGENHUB
◆ main openmsx / src / utils / halfband.hh
3.2 KiB · C++ · last changed chore(): initial commit 1 month ago source blame history raw
1 
#ifndef HALFBAND_HH
2 
#define HALFBAND_HH
3
4 
#include <cassert>
5 
#include <span>
6
7 
// Remove the upper-half of the spectrum (low-pass filter)
8 
// and then halve the sample-rate.
9 
//
10 
// The input must contain twice the number of output samples plus some extra.
11 
//   in.size() == 2 * out.size() + HALF_BAND_EXTRA
12 
//
13 
// About overlapping input and output (calculate in-place):
14 
// The output buffer is smaller than the input. Overlap is allowed, but only
15 
// when the beginning of the output is aligned with the beginning of the input.
16 
static constexpr size_t HALF_BAND_EXTRA = 13; // depends on the number of coefficients in the filter
17 
inline void halfBand(std::span<const float> in, std::span<float> out)
18 
{
19 
	// The low-pass filter is implemented via a 15-point FIR filter.
20 
	//   [c1 0 c2 0 c3 0 c4 1 c4 0 c3 0 c2 0 c1]
21 
	// Note:
22 
	// * The filter is symmetrical.
23 
	// * The middle coefficient is 1.
24 
	// * Several coefficients are 0.
25 
	//  -> This allows for an efficient implementation (only 4 multiplications per output)
26 
	// * The sum of the coefficients is 2.
27 
	//  -> This implies an amplification by a factor 2. That's undesired,
28 
	//     but saves 1 multiplication (and usually not a big problem).
29 
	//
30 
	// After filtering we drop every 2nd sample (decimate).
31 
	// (obviously we don't need to compute the outputs that will be dropped).
32 
	//
33 
	// Schematic representation:
34 
	//           x0  x1  x2  x3  x4  x5  x6  x7  x8  x9  x10 x11 x12 x13 x14 x15 x16 x17 x18 ...
35 
	//     y0 = [c1   0  c2   0  c3   0  c4   1  c4   0  c3   0  c2   0  c1]
36 
	//     y1 =         [c1   0  c2   0  c3   0  c4   1  c4   0  c3   0  c2   0  c1]
37 
	//     y2 =                 [c1   0  c2   0  c3   0  c4   1  c4   0  c3   0  c2   0  c1]
38
39 
	// See script below for how to calculate these coefficients.
40 
	static constexpr float c1 = -0.076167110507847f;
41 
	static constexpr float c2 =  0.120096364504897f;
42 
	static constexpr float c3 = -0.215886324727246f;
43 
	static constexpr float c4 =  0.671957070730196f;
44
45 
	size_t n = out.size();
46 
	assert(in.size() == (2 * n + HALF_BAND_EXTRA));
47
48 
	for (size_t i = 0; i < n; ++i) {
49 
		out[i] =    in[2 * i + 7]
50 
		    + c1 * (in[2 * i + 0] + in[2 * i + 14])
51 
		    + c2 * (in[2 * i + 2] + in[2 * i + 12])
52 
		    + c3 * (in[2 * i + 4] + in[2 * i + 10])
53 
		    + c4 * (in[2 * i + 6] + in[2 * i +  8]);
54 
	}
55 
}
56
57 
#endif
58
59 
// The filter coefficients were obtained via this octave (matlab) script:
60 
/* ----- 8< ----- 8< ----- 8< ----- 8< ----- 8< ----- 8< ----- 8< -----
61
62 
pkg load signal
63
64 
ntaps = 4 * 3 + 3;
65 
beta = 1;
66 
resolution = 1024;        % relates to the 'resolution' of the plot
67 
samplefreq = 100;
68
69
70 
N= ntaps-1;
71 
n= -N/2:N/2;
72 
sinc= sin(n*pi/2)./(n*pi+eps);      % truncated impulse response; eps= 2E-16
73 
sinc(N/2 +1)= 1/2;                  % value for n --> 0
74 
win= kaiser(ntaps, beta);           % window function
75 
b= sinc.*win';
76 
b(2:2:ntaps/2)=0;
77 
b((ntaps+5)/2:2:ntaps)=0;
78
79 
b2=b;
80 
b2((ntaps+1)/2)=0;
81 
b2 = b2 / (2 * sum(b2));
82 
b2((ntaps+1)/2)=0.5;
83 
b=b2;
84
85 
x = b;
86 
x(N+1:resolution) = 0;
87 
y = fft(x);
88
89 
xx = [0 : samplefreq/(resolution) : samplefreq/2];
90 
plot(xx, 20 * log(abs(y(1:resolution/2+1))));
91 
xlabel("frequency [kHz]");
92 
ylabel("attenuation [dB]");
93
94 
% print result
95 
format long
96 
2 * b
97
98 
----- 8< ----- 8< ----- 8< ----- 8< ----- 8< ----- 8< ----- 8< ----- */