1 /*
2 * Copyright (c) 2007 - 2015 Joseph Gaeddert
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23 //
24 // Flipped Nyquist/root-Nyquist filter designs
25 //
26 // References:
27 // [Beaulieu:2001]
28 // [Assalini:2004]
29 //
30
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <string.h>
34 #include <math.h>
35
36 #include "liquid.internal.h"
37
38 // Design flipped Nyquist/root-Nyquist filter
39 // _type : filter type (e.g. LIQUID_FIRFILT_FEXP)
40 // _root : square-root Nyquist filter?
41 // _k : samples/symbol
42 // _m : symbol delay
43 // _beta : rolloff factor (0 < beta <= 1)
44 // _dt : fractional sample delay
45 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_fnyquist(liquid_firfilt_type _type,int _root,unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)46 void liquid_firdes_fnyquist(liquid_firfilt_type _type,
47 int _root,
48 unsigned int _k,
49 unsigned int _m,
50 float _beta,
51 float _dt,
52 float * _h)
53 {
54 // validate input
55 if ( _k < 1 ) {
56 fprintf(stderr,"error: liquid_firdes_fnyquist(): k must be greater than 0\n");
57 exit(1);
58 } else if ( _m < 1 ) {
59 fprintf(stderr,"error: liquid_firdes_fnyquist(): m must be greater than 0\n");
60 exit(1);
61 } else if ( (_beta < 0.0f) || (_beta > 1.0f) ) {
62 fprintf(stderr,"error: liquid_firdes_fnyquist(): beta must be in [0,1]\n");
63 exit(1);
64 } else;
65
66 unsigned int i;
67
68 // derived values
69 unsigned int h_len = 2*_k*_m+1; // filter length
70
71 float H_prime[h_len]; // frequency response of Nyquist filter (real)
72 float complex H[h_len]; // frequency response of Nyquist filter
73 float complex h[h_len]; // impulse response of Nyquist filter
74
75 // compute Nyquist filter frequency response
76 switch (_type) {
77 case LIQUID_FIRFILT_FEXP:
78 liquid_firdes_fexp_freqresponse(_k, _m, _beta, H_prime);
79 break;
80 case LIQUID_FIRFILT_FSECH:
81 liquid_firdes_fsech_freqresponse(_k, _m, _beta, H_prime);
82 break;
83 case LIQUID_FIRFILT_FARCSECH:
84 liquid_firdes_farcsech_freqresponse(_k, _m, _beta, H_prime);
85 break;
86 default:
87 fprintf(stderr,"error: liquid_firdes_fnyquist(), unknown/unsupported filter type\n");
88 exit(1);
89 }
90
91 // copy result to fft input buffer, computing square root
92 // if required
93 for (i=0; i<h_len; i++)
94 H[i] = _root ? sqrtf(H_prime[i]) : H_prime[i];
95
96 // compute ifft
97 fft_run(h_len, H, h, LIQUID_FFT_BACKWARD, 0);
98
99 // copy shifted, scaled response
100 for (i=0; i<h_len; i++)
101 _h[i] = crealf( h[(i+_k*_m+1)%h_len] ) * (float)_k / (float)(h_len);
102 }
103
104 // Design fexp Nyquist filter
105 // _k : samples/symbol
106 // _m : symbol delay
107 // _beta : rolloff factor (0 < beta <= 1)
108 // _dt : fractional sample delay
109 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_fexp(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)110 void liquid_firdes_fexp(unsigned int _k,
111 unsigned int _m,
112 float _beta,
113 float _dt,
114 float * _h)
115 {
116 // compute resonse using generic function
117 liquid_firdes_fnyquist(LIQUID_FIRFILT_FEXP, 0, _k, _m, _beta, _dt, _h);
118 }
119
120 // Design fexp square-root Nyquist filter
121 // _k : samples/symbol
122 // _m : symbol delay
123 // _beta : rolloff factor (0 < beta <= 1)
124 // _dt : fractional sample delay
125 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_rfexp(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)126 void liquid_firdes_rfexp(unsigned int _k,
127 unsigned int _m,
128 float _beta,
129 float _dt,
130 float * _h)
131 {
132 // compute resonse using generic function
133 liquid_firdes_fnyquist(LIQUID_FIRFILT_FEXP, 1, _k, _m, _beta, _dt, _h);
134 }
135
136 // flipped exponential frequency response
liquid_firdes_fexp_freqresponse(unsigned int _k,unsigned int _m,float _beta,float * _H)137 void liquid_firdes_fexp_freqresponse(unsigned int _k,
138 unsigned int _m,
139 float _beta,
140 float * _H)
141 {
142 // TODO : validate input
143
144 unsigned int i;
145
146 unsigned int h_len = 2*_k*_m + 1;
147
148 float f0 = 0.5f*(1.0f - _beta) / (float)_k;
149 float f1 = 0.5f*(1.0f ) / (float)_k;
150 float f2 = 0.5f*(1.0f + _beta) / (float)_k;
151
152 float B = 0.5f/(float)_k;
153 float gamma = logf(2.0f)/(_beta*B);
154
155 // compute frequency response of Nyquist filter
156 for (i=0; i<h_len; i++) {
157 float f = (float)i / (float)h_len;
158 if (f > 0.5f) f = f - 1.0f;
159
160 // enforce even symmetry
161 f = fabsf(f);
162
163 if ( f < f0 ) {
164 // pass band
165 _H[i] = 1.0f;
166 } else if (f > f0 && f < f2) {
167 // transition band
168 if ( f < f1) {
169 _H[i] = expf(gamma*(B*(1-_beta) - f));
170 } else {
171 _H[i] = 1.0f - expf(gamma*(f - (1+_beta)*B));
172 }
173 } else {
174 // stop band
175 _H[i] = 0.0f;
176 }
177 }
178 }
179
180 // Design fsech Nyquist filter
181 // _k : samples/symbol
182 // _m : symbol delay
183 // _beta : rolloff factor (0 < beta <= 1)
184 // _dt : fractional sample delay
185 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_fsech(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)186 void liquid_firdes_fsech(unsigned int _k,
187 unsigned int _m,
188 float _beta,
189 float _dt,
190 float * _h)
191 {
192 // compute resonse using generic function
193 liquid_firdes_fnyquist(LIQUID_FIRFILT_FSECH, 0, _k, _m, _beta, _dt, _h);
194 }
195
196 // Design fsech square-root Nyquist filter
197 // _k : samples/symbol
198 // _m : symbol delay
199 // _beta : rolloff factor (0 < beta <= 1)
200 // _dt : fractional sample delay
201 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_rfsech(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)202 void liquid_firdes_rfsech(unsigned int _k,
203 unsigned int _m,
204 float _beta,
205 float _dt,
206 float * _h)
207 {
208 // compute resonse using generic function
209 liquid_firdes_fnyquist(LIQUID_FIRFILT_FSECH, 1, _k, _m, _beta, _dt, _h);
210 }
211
212 // flipped exponential frequency response
liquid_firdes_fsech_freqresponse(unsigned int _k,unsigned int _m,float _beta,float * _H)213 void liquid_firdes_fsech_freqresponse(unsigned int _k,
214 unsigned int _m,
215 float _beta,
216 float * _H)
217 {
218 // TODO : validate input
219
220 unsigned int i;
221
222 unsigned int h_len = 2*_k*_m + 1;
223
224 float f0 = 0.5f*(1.0f - _beta) / (float)_k;
225 float f1 = 0.5f*(1.0f ) / (float)_k;
226 float f2 = 0.5f*(1.0f + _beta) / (float)_k;
227
228 float B = 0.5f/(float)_k;
229 float gamma = logf(sqrtf(3.0f) + 2.0f) / (_beta*B);
230
231 // compute frequency response of Nyquist filter
232 for (i=0; i<h_len; i++) {
233 float f = (float)i / (float)h_len;
234 if (f > 0.5f) f = f - 1.0f;
235
236 // enforce even symmetry
237 f = fabsf(f);
238
239 if ( f < f0 ) {
240 // pass band
241 _H[i] = 1.0f;
242 } else if (f > f0 && f < f2) {
243 // transition band
244 if ( f < f1) {
245 _H[i] = 1.0f / coshf(gamma*(f - B*(1-_beta)));
246 } else {
247 _H[i] = 1.0f - 1.0f / coshf(gamma*(B*(1+_beta) - f));
248 }
249 } else {
250 // stop band
251 _H[i] = 0.0f;
252 }
253 }
254 }
255
256 // Design farcsech Nyquist filter
257 // _k : samples/symbol
258 // _m : symbol delay
259 // _beta : rolloff factor (0 < beta <= 1)
260 // _dt : fractional sample delay
261 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_farcsech(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)262 void liquid_firdes_farcsech(unsigned int _k,
263 unsigned int _m,
264 float _beta,
265 float _dt,
266 float * _h)
267 {
268 // compute resonse using generic function
269 liquid_firdes_fnyquist(LIQUID_FIRFILT_FARCSECH, 0, _k, _m, _beta, _dt, _h);
270 }
271
272 // Design farcsech square-root Nyquist filter
273 // _k : samples/symbol
274 // _m : symbol delay
275 // _beta : rolloff factor (0 < beta <= 1)
276 // _dt : fractional sample delay
277 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_rfarcsech(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)278 void liquid_firdes_rfarcsech(unsigned int _k,
279 unsigned int _m,
280 float _beta,
281 float _dt,
282 float * _h)
283 {
284 // compute resonse using generic function
285 liquid_firdes_fnyquist(LIQUID_FIRFILT_FARCSECH, 1, _k, _m, _beta, _dt, _h);
286 }
287
288 // hyperbolic arc-secant
liquid_asechf(float _z)289 float liquid_asechf(float _z)
290 {
291 if (_z <= 0.0f || _z > 1.0f) {
292 fprintf(stderr,"warning: liquid_asechf(), input out of range\n");
293 return 0.0f;
294 }
295
296 float z_inv = 1.0f / _z;
297
298 return logf( sqrtf(z_inv - 1.0f)*sqrtf(z_inv + 1.0f) + z_inv );
299 }
300
301 // flipped exponential frequency response
liquid_firdes_farcsech_freqresponse(unsigned int _k,unsigned int _m,float _beta,float * _H)302 void liquid_firdes_farcsech_freqresponse(unsigned int _k,
303 unsigned int _m,
304 float _beta,
305 float * _H)
306 {
307 // TODO : validate input
308
309 unsigned int i;
310
311 unsigned int h_len = 2*_k*_m + 1;
312
313 float f0 = 0.5f*(1.0f - _beta) / (float)_k;
314 float f1 = 0.5f*(1.0f ) / (float)_k;
315 float f2 = 0.5f*(1.0f + _beta) / (float)_k;
316
317 float B = 0.5f/(float)_k;
318 float gamma = logf(sqrtf(3.0f) + 2.0f) / (_beta*B);
319 float zeta = 1.0f / (2.0f * _beta * B);
320
321 // compute frequency response of Nyquist filter
322 for (i=0; i<h_len; i++) {
323 float f = (float)i / (float)h_len;
324 if (f > 0.5f) f = f - 1.0f;
325
326 // enforce even symmetry
327 f = fabsf(f);
328
329 if ( f < f0 ) {
330 // pass band
331 _H[i] = 1.0f;
332 } else if (f > f0 && f < f2) {
333 // transition band
334 if ( f < f1) {
335 _H[i] = 1.0f - (zeta/gamma)*liquid_asechf(zeta*(B*(1+_beta) - f));
336 } else {
337 _H[i] = (zeta/gamma)*liquid_asechf(zeta*(f - B*(1-_beta)));
338 }
339 } else {
340 // stop band
341 _H[i] = 0.0f;
342 }
343 }
344 }
345
346