/ libcsdr.h
libcsdr.h
1 /* 2 This software is part of libcsdr, a set of simple DSP routines for 3 Software Defined Radio. 4 5 Copyright (c) 2014, Andras Retzler <randras@sdr.hu> 6 All rights reserved. 7 8 Redistribution and use in source and binary forms, with or without 9 modification, are permitted provided that the following conditions are met: 10 * Redistributions of source code must retain the above copyright 11 notice, this list of conditions and the following disclaimer. 12 * Redistributions in binary form must reproduce the above copyright 13 notice, this list of conditions and the following disclaimer in the 14 documentation and/or other materials provided with the distribution. 15 * Neither the name of the copyright holder nor the 16 names of its contributors may be used to endorse or promote products 17 derived from this software without specific prior written permission. 18 19 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 20 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 21 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 22 DISCLAIMED. IN NO EVENT SHALL ANDRAS RETZLER BE LIABLE FOR ANY 23 DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 24 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 25 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 26 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 27 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 28 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 29 */ 30 31 #pragma once 32 #define MIN_M(x,y) (((x)>(y))?(y):(x)) 33 #define MAX_M(x,y) (((x)<(y))?(y):(x)) 34 35 /* 36 _____ _ 37 / ____| | | 38 | | ___ _ __ ___ _ __ | | _____ __ 39 | | / _ \| '_ ` _ \| '_ \| |/ _ \ \/ / 40 | |___| (_) | | | | | | |_) | | __/> < 41 \_____\___/|_| |_| |_| .__/|_|\___/_/\_\ 42 | | 43 |_| 44 */ 45 46 typedef struct complexf_s { float i; float q; } complexf; 47 48 //apply to pointers: 49 #define iof(complexf_input_p,i) (*(((float*)complexf_input_p)+2*(i))) 50 #define qof(complexf_input_p,i) (*(((float*)complexf_input_p)+2*(i)+1)) 51 #define absof(complexf_input_p,i) (sqrt((iof(complexf_input_p,i)*iof(complexf_input_p,i))+(qof(complexf_input_p,i)*qof(complexf_input_p,i)))) 52 #define argof(complexf_input_p,i) (atan2(qof(complexf_input_p,i),iof(complexf_input_p,i))) 53 #define cmult(cfo, cfi1, cfi2) {iof(cfo,0)=iof(cfi1,0)*iof(cfi2,0)-qof(cfi1,0)*qof(cfi2,0);qof(cfo,0)=iof(cfi1,0)*qof(cfi2,0)+iof(cfi2,0)*qof(cfi1,0);} 54 //(ai+aq*j)*(bi+bq*j)=ai*bi-aq*bq+(aq*bi+ai*bq)*j 55 #define cmultadd(cfo, cfi1, cfi2) { iof(cfo,0)+=iof(cfi1,0)*iof(cfi2,0)-qof(cfi1,0)*qof(cfi2,0);qof(cfo,0)+=iof(cfi1,0)*qof(cfi2,0)+iof(cfi2,0)*qof(cfi1,0); } 56 #define csetnull(cf) { iof(cf,0)=0.0; qof(cf,0)=0.0; } 57 #define e_powj(cf,w) { iof(cf,0)=cos(w); qof(cf,0)=sin(w); } 58 #define ccopy(dst,src) { iof(dst,0)=iof(src,0); qof(dst,0)=qof(src,0); } 59 60 //apply to values 61 #define iofv(complexf_input) (*((float*)&complexf_input)) 62 #define qofv(complexf_input) (*(((float*)&complexf_input)+1)) 63 64 //they dropped M_PI in C99, so we define it: 65 #define PI ((float)3.14159265358979323846) 66 67 #define TIME_TAKEN(start,end) ((end.tv_sec-start.tv_sec)+(end.tv_nsec-start.tv_nsec)/1e9) 68 69 //window 70 typedef enum window_s 71 { 72 WINDOW_BOXCAR, WINDOW_BLACKMAN, WINDOW_HAMMING 73 } window_t; 74 75 #define WINDOW_DEFAULT WINDOW_HAMMING 76 77 //FFT 78 //Note: these reference to things in this file (e.g. complexf): 79 #include "fft_fftw.h" 80 #include "fft_rpi.h" 81 82 // ================================================================================= 83 84 //filter design 85 void firdes_lowpass_f(float *output, int length, float cutoff_rate, window_t window); 86 void firdes_bandpass_c(complexf *output, int length, float lowcut, float highcut, window_t window); 87 float firdes_wkernel_blackman(float input); 88 float firdes_wkernel_hamming(float input); 89 float firdes_wkernel_boxcar(float input); 90 window_t firdes_get_window_from_string(char* input); 91 char* firdes_get_string_from_window(window_t window); 92 int firdes_filter_len(float transition_bw); 93 94 //demodulators 95 complexf fmdemod_quadri_cf(complexf* input, float* output, int input_size, float *temp, complexf last_sample); 96 complexf fmdemod_quadri_novect_cf(complexf* input, float* output, int input_size, complexf last_sample); 97 float fmdemod_atan_cf(complexf* input, float *output, int input_size, float last_phase); 98 void amdemod_cf(complexf* input, float *output, int input_size); 99 void amdemod_estimator_cf(complexf* input, float *output, int input_size, float alpha, float beta); 100 void limit_ff(float* input, float* output, int input_size, float max_amplitude); 101 102 //filters, decimators, resamplers, shift, etc. 103 float fir_one_pass_ff(float* input, float* taps, int taps_length); 104 int fir_decimate_cc(complexf *input, complexf *output, int input_size, int decimation, float *taps, int taps_length); 105 int fir_interpolate_cc(complexf *input, complexf *output, int input_size, int interpolation, float *taps, int taps_length); 106 int deemphasis_nfm_ff (float* input, float* output, int input_size, int sample_rate); 107 float deemphasis_wfm_ff (float* input, float* output, int input_size, float tau, int sample_rate, float last_output); 108 float shift_math_cc(complexf *input, complexf* output, int input_size, float rate, float starting_phase); 109 110 typedef struct dcblock_preserve_s 111 { 112 float last_input; 113 float last_output; 114 } dcblock_preserve_t; 115 dcblock_preserve_t dcblock_ff(float* input, float* output, int input_size, float a, dcblock_preserve_t preserved); 116 float fastdcblock_ff(float* input, float* output, int input_size, float last_dc_level); 117 118 typedef struct fastagc_ff_s 119 { 120 float* buffer_1; 121 float* buffer_2; 122 float* buffer_input; //it is the actual input buffer to fill 123 float peak_1; 124 float peak_2; 125 int input_size; 126 float reference; 127 float last_gain; 128 } fastagc_ff_t; 129 130 void fastagc_ff(fastagc_ff_t* input, float* output); 131 132 typedef struct rational_resampler_ff_s 133 { 134 int input_processed; 135 int output_size; 136 int last_taps_delay; 137 } rational_resampler_ff_t; 138 139 rational_resampler_ff_t rational_resampler_ff(float *input, float *output, int input_size, int interpolation, int decimation, float *taps, int taps_length, int last_taps_delay); 140 void rational_resampler_get_lowpass_f(float* output, int output_size, int interpolation, int decimation, window_t window); 141 142 float *precalculate_window(int size, window_t window); 143 void apply_window_c(complexf* input, complexf* output, int size, window_t window); 144 void apply_precalculated_window_c(complexf* input, complexf* output, int size, float *windowt); 145 void apply_precalculated_window_f(float* input, float* output, int size, float *windowt); 146 void apply_window_f(float* input, float* output, int size, window_t window); 147 void logpower_cf(complexf* input, float* output, int size, float add_db); 148 void accumulate_power_cf(complexf* input, float* output, int size); 149 void log_ff(float* input, float* output, int size, float add_db); 150 151 typedef struct fractional_decimator_ff_s 152 { 153 float where; 154 int input_processed; 155 int output_size; 156 int num_poly_points; //number of samples that the Lagrange interpolator will use 157 float* poly_precalc_denomiator; //while we don't precalculate coefficients here as in a Farrow structure, because it is a fractional interpolator, but we rather precaculate part of the interpolator expression 158 //float* last_inputs_circbuf; //circular buffer to store the last (num_poly_points) number of input samples. 159 //int last_inputs_startsat; //where the circular buffer starts now 160 //int last_inputs_samplewhere; 161 float* coeffs_buf; 162 float* filtered_buf; 163 int xifirst; 164 int xilast; 165 float rate; 166 float *taps; 167 int taps_length; 168 } fractional_decimator_ff_t; 169 fractional_decimator_ff_t fractional_decimator_ff_init(float rate, int num_poly_points, float* taps, int taps_length); 170 void fractional_decimator_ff(float* input, float* output, int input_size, fractional_decimator_ff_t* d); 171 172 typedef struct old_fractional_decimator_ff_s 173 { 174 float remain; 175 int input_processed; 176 int output_size; 177 } old_fractional_decimator_ff_t; 178 old_fractional_decimator_ff_t old_fractional_decimator_ff(float* input, float* output, int input_size, float rate, float *taps, int taps_length, old_fractional_decimator_ff_t d); 179 180 typedef struct shift_table_data_s 181 { 182 float* table; 183 int table_size; 184 } shift_table_data_t; 185 void shift_table_deinit(shift_table_data_t table_data); 186 shift_table_data_t shift_table_init(int table_size); 187 float shift_table_cc(complexf* input, complexf* output, int input_size, float rate, shift_table_data_t table_data, float starting_phase); 188 189 typedef struct shift_addfast_data_s 190 { 191 float dsin[4]; 192 float dcos[4]; 193 float phase_increment; 194 } shift_addfast_data_t; 195 shift_addfast_data_t shift_addfast_init(float rate); 196 shift_addfast_data_t shift_addfast_init(float rate); 197 float shift_addfast_cc(complexf *input, complexf* output, int input_size, shift_addfast_data_t* d, float starting_phase); 198 199 typedef struct shift_unroll_data_s 200 { 201 float* dsin; 202 float* dcos; 203 float phase_increment; 204 int size; 205 } shift_unroll_data_t; 206 float shift_unroll_cc(complexf *input, complexf* output, int input_size, shift_unroll_data_t* d, float starting_phase); 207 shift_unroll_data_t shift_unroll_init(float rate, int size); 208 209 int log2n(int x); 210 int next_pow2(int x); 211 void apply_fir_fft_cc(FFT_PLAN_T* plan, FFT_PLAN_T* plan_inverse, complexf* taps_fft, complexf* last_overlap, int overlap_size); 212 void gain_ff(float* input, float* output, int input_size, float gain); 213 float get_power_f(float* input, int input_size, int decimation); 214 float get_power_c(complexf* input, int input_size, int decimation); 215 216 void add_dcoffset_cc(complexf* input, complexf* output, int input_size); 217 float fmmod_fc(float* input, complexf* output, int input_size, float last_phase); 218 void fixed_amplitude_cc(complexf* input, complexf* output, int input_size, float amp); 219 220 void convert_u8_f(unsigned char* input, float* output, int input_size); 221 void convert_f_u8(float* input, unsigned char* output, int input_size); 222 void convert_s8_f(signed char* input, float* output, int input_size); 223 void convert_f_s8(float* input, signed char* output, int input_size); 224 void convert_f_s16(float* input, short* output, int input_size); 225 void convert_s16_f(short* input, float* output, int input_size); 226 void convert_f_i16(float* input, short* output, int input_size); 227 void convert_i16_f(short* input, float* output, int input_size); 228 void convert_f_s24(float* input, unsigned char* output, int input_size, int bigendian); 229 void convert_s24_f(unsigned char* input, float* output, int input_size, int bigendian); 230 231 232 int is_nan(float f); 233 234 //digital demod 235 236 typedef struct rtty_baudot_item_s 237 { 238 unsigned long long code; 239 unsigned char ascii_letter; 240 unsigned char ascii_figure; 241 } rtty_baudot_item_t; 242 243 typedef enum rtty_baudot_decoder_state_e 244 { 245 RTTY_BAUDOT_WAITING_STOP_PULSE = 0, 246 RTTY_BAUDOT_WAITING_START_PULSE, 247 RTTY_BAUDOT_RECEIVING_DATA 248 } rtty_baudot_decoder_state_t; 249 250 typedef struct rtty_baudot_decoder_s 251 { 252 unsigned char fig_mode; 253 unsigned char character_received; 254 unsigned short shr; 255 unsigned char bit_cntr; 256 rtty_baudot_decoder_state_t state; 257 } rtty_baudot_decoder_t; 258 259 #define RTTY_FIGURE_MODE_SELECT_CODE 0b11011 260 #define RTTY_LETTER_MODE_SELECT_CODE 0b11111 261 262 char rtty_baudot_decoder_lookup(unsigned char* fig_mode, unsigned char c); 263 char rtty_baudot_decoder_push(rtty_baudot_decoder_t* s, unsigned char symbol); 264 265 //PSK31 266 267 typedef struct psk31_varicode_item_s 268 { 269 unsigned long long code; 270 int bitcount; 271 unsigned char ascii; 272 } psk31_varicode_item_t; 273 274 char psk31_varicode_decoder_push(unsigned long long* status_shr, unsigned char symbol); 275 276 //Serial 277 278 typedef struct serial_line_s 279 { 280 float samples_per_bits; 281 int databits; //including parity 282 float stopbits; 283 int output_size; 284 int input_used; 285 float bit_sampling_width_ratio; 286 } serial_line_t; 287 288 void serial_line_decoder_f_u8(serial_line_t* s, float* input, unsigned char* output, int input_size); 289 void binary_slicer_f_u8(float* input, unsigned char* output, int input_size); 290 291 292 typedef enum pll_type_e 293 { 294 PLL_P_CONTROLLER=1, 295 PLL_PI_CONTROLLER=2 296 } pll_type_t; 297 298 typedef struct pll_s 299 { 300 pll_type_t pll_type; 301 //common: 302 float output_phase; 303 float dphase; 304 float frequency; 305 float alpha; 306 float beta; 307 float iir_temp; 308 } pll_t; 309 310 void pll_cc_init_pi_controller(pll_t* p, float bandwidth, float ko, float kd, float damping_factor); 311 void pll_cc_init_p_controller(pll_t* p, float alpha); 312 void pll_cc(pll_t* p, complexf* input, float* output_dphase, complexf* output_nco, int input_size); 313 314 typedef enum timing_recovery_algorithm_e 315 { 316 TIMING_RECOVERY_ALGORITHM_GARDNER, 317 TIMING_RECOVERY_ALGORITHM_EARLYLATE 318 } timing_recovery_algorithm_t; 319 320 #define TIMING_RECOVERY_ALGORITHM_DEFAULT TIMING_RECOVERY_ALGORITHM_GARDNER 321 322 typedef struct timing_recovery_state_s 323 { 324 timing_recovery_algorithm_t algorithm; 325 int decimation_rate; // = input_rate / output_rate. We should get an input signal that is N times oversampled. 326 int output_size; 327 int input_processed; 328 int use_q; //use both I and Q for calculating the error 329 int debug_phase; 330 int debug_every_nth; 331 char* debug_writefiles_path; 332 int last_correction_offset; 333 float earlylate_ratio; 334 float loop_gain; 335 float max_error; 336 } timing_recovery_state_t; 337 338 timing_recovery_state_t timing_recovery_init(timing_recovery_algorithm_t algorithm, int decimation_rate, int use_q, float loop_gain, float max_error, int debug_every_nth, char* debug_writefiles_path); 339 void timing_recovery_cc(complexf* input, complexf* output, int input_size, float* timing_error, int* sampled_indexes, timing_recovery_state_t* state); 340 timing_recovery_algorithm_t timing_recovery_get_algorithm_from_string(char* input); 341 char* timing_recovery_get_string_from_algorithm(timing_recovery_algorithm_t algorithm); 342 void octave_plot_point_on_cplxsig(complexf* signal, int signal_size, float error, int index, int correction_offset, char* writefiles_path, int points_size, ...); 343 void psk_modulator_u8_c(unsigned char* input, complexf* output, int input_size, int n_psk); 344 void duplicate_samples_ntimes_u8_u8(unsigned char* input, unsigned char* output, int input_size_bytes, int sample_size_bytes, int ntimes); 345 complexf psk31_interpolate_sine_cc(complexf* input, complexf* output, int input_size, int interpolation, complexf last_input); 346 void psk31_varicode_encoder_u8_u8(unsigned char* input, unsigned char* output, int input_size, int output_max_size, int* input_processed, int* output_size); 347 unsigned char differential_codec(unsigned char* input, unsigned char* output, int input_size, int encode, unsigned char state); 348 349 #if 0 350 typedef struct bpsk_costas_loop_state_s 351 { 352 float rc_filter_alpha; 353 float vco_phase_addition_multiplier; 354 float vco_phase; 355 float last_lpfi_output; 356 float last_lpfq_output; 357 float last_vco_phase_addition; 358 } bpsk_costas_loop_state_t; 359 360 bpsk_costas_loop_state_t init_bpsk_costas_loop_cc(float samples_per_bits); 361 void bpsk_costas_loop_cc(complexf* input, complexf* output, int input_size, bpsk_costas_loop_state_t* state); 362 #endif 363 364 typedef struct bpsk_costas_loop_state_s 365 { 366 float alpha; 367 float beta; 368 int decision_directed; 369 float current_freq; 370 float dphase; 371 float nco_phase; 372 float dphase_max; 373 int dphase_max_reset_to_zero; 374 } bpsk_costas_loop_state_t; 375 376 void plain_interpolate_cc(complexf* input, complexf* output, int input_size, int interpolation); 377 void bpsk_costas_loop_cc(complexf* input, complexf* output, int input_size, float* output_error, float* output_dphase, complexf* output_nco, bpsk_costas_loop_state_t* s); 378 void init_bpsk_costas_loop_cc(bpsk_costas_loop_state_t* s, int decision_directed, float damping_factor, float bandwidth); 379 380 void simple_agc_cc(complexf* input, complexf* output, int input_size, float rate, float reference, float max_gain, float* current_gain); 381 void firdes_add_peak_c(complexf* output, int length, float rate, window_t window, int add, int normalize); 382 int apply_fir_cc(complexf* input, complexf* output, int input_size, complexf* taps, int taps_length); 383 384 385 FILE* init_get_random_samples_f(); 386 void get_random_samples_f(float* output, int output_size, FILE* status); 387 void get_random_gaussian_samples_c(complexf* output, int output_size, FILE* status); 388 int deinit_get_random_samples_f(FILE* status); 389 float* add_ff(float* input1, float* input2, float* output, int input_size); 390 float total_logpower_cf(complexf* input, int input_size); 391 float normalized_timing_variance_u32_f(unsigned* input, float* temp, int input_size, int samples_per_symbol, int initial_sample_offset, int debug_print); 392 393 typedef enum matched_filter_type_e 394 { 395 MATCHED_FILTER_RRC, 396 MATCHED_FILTER_COSINE 397 } matched_filter_type_t; 398 399 #define MATCHED_FILTER_DEFAULT MATCHED_FILTER_RRC 400 401 int firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol); 402 int firdes_rrc_f(float* taps, int taps_length, int samples_per_symbol, float beta); 403 matched_filter_type_t matched_filter_get_type_from_string(char* input); 404 int apply_real_fir_cc(complexf* input, complexf* output, int input_size, float* taps, int taps_length); 405 void generic_slicer_f_u8(float* input, unsigned char* output, int input_size, int n_symbols); 406 void plain_interpolate_cc(complexf* input, complexf* output, int input_size, int interpolation);; 407 void normalize_fir_f(float* input, float* output, int length); 408 float* add_const_cc(complexf* input, complexf* output, int input_size, complexf x); 409 void pack_bits_1to8_u8_u8(unsigned char* input, unsigned char* output, int input_size); 410 unsigned char pack_bits_8to1_u8_u8(unsigned char* input); 411 void dbpsk_decoder_c_u8(complexf* input, unsigned char* output, int input_size); 412 int bfsk_demod_cf(complexf* input, float* output, int input_size, complexf* mark_filter, complexf* space_filter, int taps_length);