88
99namespace dsplib {
1010
11- // -------------------------------------------------------------------------------------------------
11+ namespace {
12+
1213template <class T >
1314static void _conv (const T* restrict x, const T* restrict h, T* restrict r, int nh, int nx) {
1415 const int nr = nx - nh + 1 ;
@@ -21,6 +22,92 @@ static void _conv(const T* restrict x, const T* restrict h, T* restrict r, int n
2122 }
2223}
2324
25+ arr_real _lowpass_fir (int n, real_t wn, const arr_real& win) {
26+ if (win.size () != (n + 1 )) {
27+ DSPLIB_THROW (" Window must be n+1 elements" );
28+ }
29+
30+ const bool is_odd = (n % 2 == 1 );
31+ const int M = win.size ();
32+ const int L = M / 2 ;
33+ const real_t fc = wn / 2 ;
34+ const arr_real w = win.slice (0 , L);
35+
36+ auto tt = arange (L) - real_t (n) / 2 ;
37+ auto h = zeros (M);
38+ h.slice (0 , L) = sin (2 * pi * fc * tt) / tt * w;
39+ if (!is_odd) {
40+ h (L) = 2 * pi * fc;
41+ h.slice (L + 1 , M) = flip (h.slice (0 , L));
42+ } else {
43+ h.slice (L, M) = flip (h.slice (0 , L));
44+ }
45+
46+ h /= sum (h);
47+ return h;
48+ }
49+
50+ arr_real _highpass_fir (int n, real_t wn, const arr_real& win) {
51+ wn = 1 - wn;
52+ int t1 = 0 ;
53+ if (n % 2 == 1 ) {
54+ n += 1 ;
55+ t1 = 1 ;
56+ }
57+
58+ auto h = _lowpass_fir (n, wn, win);
59+ auto hh = arr_real (h.slice (t1, n, 2 ));
60+ h.slice (t1, n, 2 ) = -hh;
61+ return h;
62+ }
63+
64+ arr_real _bandpass_fir (int n, real_t wn1, real_t wn2, const arr_real& win) {
65+ wn1 = wn1 / 2 ;
66+ wn2 = wn2 / 2 ;
67+ auto wp = (wn2 - wn1) / 2 ;
68+ auto wc = wn1 + wp;
69+ auto h = _lowpass_fir (n, 2 * wp, win);
70+ auto t = arange (h.size ()) - real_t (n) / 2 ;
71+ h = 2 * h * cos (2 * pi * wc * t);
72+ return h;
73+ }
74+
75+ arr_real _bandstop_fir (int n, real_t wn1, real_t wn2, const arr_real& win) {
76+ if (n % 2 == 1 ) {
77+ n += 1 ;
78+ }
79+
80+ auto h = (-1 ) * _bandpass_fir (n, wn1, wn2, win);
81+ h (n / 2 ) = h (n / 2 ) + 1 ;
82+ return h;
83+ }
84+
85+ bool _equal (const real_t & x1, const real_t & x2) {
86+ return (abs (x1 - x2) < 2 * eps ());
87+ }
88+
89+ bool _is_symmetric (const dsplib::arr_real& h) {
90+ const int n = h.size ();
91+ for (int i = 0 ; i < n / 2 ; i++) {
92+ if (!_equal (h[i], h[n - i - 1 ])) {
93+ return false ;
94+ }
95+ }
96+ return true ;
97+ }
98+
99+ bool _is_antisymmetric (const dsplib::arr_real& h) {
100+ const int n = h.size ();
101+ for (int i = 0 ; i < n / 2 ; i++) {
102+ if (!_equal (h[i], -h[n - i - 1 ])) {
103+ return false ;
104+ }
105+ }
106+ return true ;
107+ }
108+
109+ } // namespace
110+
24111// -------------------------------------------------------------------------------------------------
25112template <>
26113arr_real FirFilter<real_t >::conv(const arr_real& x, const arr_real& h) {
@@ -29,7 +116,6 @@ arr_real FirFilter<real_t>::conv(const arr_real& x, const arr_real& h) {
29116 return r;
30117}
31118
32- // -------------------------------------------------------------------------------------------------
33119template <>
34120arr_cmplx FirFilter<cmplx_t >::conv(const arr_cmplx& x, const arr_cmplx& h) {
35121 arr_cmplx r (x.size () - h.size () + 1 );
@@ -48,12 +134,10 @@ FftFilter::FftFilter(const arr_cmplx& h)
48134 _x = complex (zeros (fft_len));
49135}
50136
51- // -------------------------------------------------------------------------------------------------
52137FftFilter::FftFilter (const arr_real& h)
53138 : FftFilter(complex (h)) {
54139}
55140
56- // -------------------------------------------------------------------------------------------------
57141arr_cmplx FftFilter::process (const arr_cmplx& x) {
58142 const int nr = (x.size () + _nx) / _n * _n;
59143 arr_cmplx r (nr);
@@ -80,76 +164,11 @@ arr_cmplx FftFilter::process(const arr_cmplx& x) {
80164 return r;
81165}
82166
83- // -------------------------------------------------------------------------------------------------
84167arr_real FftFilter::process (const arr_real& x) {
85168 // TODO: real optimization
86169 return real (process (complex (x)));
87170}
88171
89- // ----------------------------------------------------------------------------------------------
90- static arr_real _lowpass_fir (int n, real_t wn, const arr_real& win) {
91- if (win.size () != (n + 1 )) {
92- DSPLIB_THROW (" Window must be n+1 elements" );
93- }
94-
95- const bool is_odd = (n % 2 == 1 );
96- const int M = win.size ();
97- const int L = M / 2 ;
98- const real_t fc = wn / 2 ;
99- const arr_real w = win.slice (0 , L);
100-
101- auto tt = arange (L) - real_t (n) / 2 ;
102- auto h = zeros (M);
103- h.slice (0 , L) = sin (2 * pi * fc * tt) / tt * w;
104- if (!is_odd) {
105- h (L) = 2 * pi * fc;
106- h.slice (L + 1 , M) = flip (h.slice (0 , L));
107- } else {
108- h.slice (L, M) = flip (h.slice (0 , L));
109- }
110-
111- h /= sum (h);
112- return h;
113- }
114-
115- // ----------------------------------------------------------------------------------------------
116- static arr_real _highpass_fir (int n, real_t wn, const arr_real& win) {
117- wn = 1 - wn;
118- int t1 = 0 ;
119- if (n % 2 == 1 ) {
120- n += 1 ;
121- t1 = 1 ;
122- }
123-
124- auto h = _lowpass_fir (n, wn, win);
125- auto hh = arr_real (h.slice (t1, n, 2 ));
126- h.slice (t1, n, 2 ) = -hh;
127- return h;
128- }
129-
130- // ----------------------------------------------------------------------------------------------
131- static arr_real _bandpass_fir (int n, real_t wn1, real_t wn2, const arr_real& win) {
132- wn1 = wn1 / 2 ;
133- wn2 = wn2 / 2 ;
134- auto wp = (wn2 - wn1) / 2 ;
135- auto wc = wn1 + wp;
136- auto h = _lowpass_fir (n, 2 * wp, win);
137- auto t = arange (h.size ()) - real_t (n) / 2 ;
138- h = 2 * h * cos (2 * pi * wc * t);
139- return h;
140- }
141-
142- // ----------------------------------------------------------------------------------------------
143- static arr_real _bandstop_fir (int n, real_t wn1, real_t wn2, const arr_real& win) {
144- if (n % 2 == 1 ) {
145- n += 1 ;
146- }
147-
148- auto h = (-1 ) * _bandpass_fir (n, wn1, wn2, win);
149- h (n / 2 ) = h (n / 2 ) + 1 ;
150- return h;
151- }
152-
153172// ----------------------------------------------------------------------------------------------
154173arr_real fir1 (int n, real_t wn, FilterType ftype, const arr_real& win) {
155174 assert (n > 0 );
@@ -166,7 +185,6 @@ arr_real fir1(int n, real_t wn, FilterType ftype, const arr_real& win) {
166185 DSPLIB_THROW (" Not supported for current filter type" );
167186}
168187
169- // ----------------------------------------------------------------------------------------------
170188arr_real fir1 (int n, real_t wn1, real_t wn2, FilterType ftype, const arr_real& win) {
171189 assert (n > 0 );
172190 assert (wn1 < wn2);
@@ -200,33 +218,6 @@ arr_real fir1(int n, real_t wn1, real_t wn2, FilterType ftype) {
200218 return fir1 (n, wn1, wn2, ftype, window::hamming (nn));
201219}
202220
203- // ----------------------------------------------------------------------------------------------
204- static bool _equal (const real_t & x1, const real_t & x2) {
205- return (abs (x1 - x2) < 2 * eps ());
206- }
207-
208- // ----------------------------------------------------------------------------------------------
209- static bool _is_symmetric (const dsplib::arr_real& h) {
210- const int n = h.size ();
211- for (int i = 0 ; i < n / 2 ; i++) {
212- if (!_equal (h[i], h[n - i - 1 ])) {
213- return false ;
214- }
215- }
216- return true ;
217- }
218-
219- // ----------------------------------------------------------------------------------------------
220- static bool _is_antisymmetric (const dsplib::arr_real& h) {
221- const int n = h.size ();
222- for (int i = 0 ; i < n / 2 ; i++) {
223- if (!_equal (h[i], -h[n - i - 1 ])) {
224- return false ;
225- }
226- }
227- return true ;
228- }
229-
230221// ----------------------------------------------------------------------------------------------
231222FirType firtype (const dsplib::arr_real& h) {
232223 const int n = h.size ();
0 commit comments