diff --git a/libfreedv/CMakeLists.txt b/libfreedv/CMakeLists.txt index 94d558a20..f32460e4e 100644 --- a/libfreedv/CMakeLists.txt +++ b/libfreedv/CMakeLists.txt @@ -9,8 +9,14 @@ set(freedv_SOURCES freedv_filter.cpp freedv_vhf_framing.cpp fsk.cpp + gp_interleaver.cpp + HRA_112_112.cpp + interldpc.cpp kiss_fft.cpp linreg.cpp + mpdecode_core.cpp + ofdm.cpp + phi0.cpp ) set(freedv_HEADERS @@ -32,6 +38,7 @@ set(freedv_HEADERS fsk.h gp_interleaver.h hanning.h + HRA_112_112.h interldpc.h _kiss_fft_guts.h kiss_fft.h @@ -44,6 +51,7 @@ set(freedv_HEADERS mpdecode_core.h ofdm_internal.h os.h + phi0.h pilot_coeff.h pilots_coh.h rn_coh.h diff --git a/libfreedv/HRA_112_112.cpp b/libfreedv/HRA_112_112.cpp new file mode 100644 index 000000000..406e11225 --- /dev/null +++ b/libfreedv/HRA_112_112.cpp @@ -0,0 +1,27 @@ +/* + FILE....: HRA_112_112.c + + Static arrays for LDPC codec HRA_112_112, generated by ldpc_gen_c_h_file.m. +*/ + +#include +#include "HRA_112_112.h" + +namespace FreeDV +{ + +const uint16_t HRA_112_112_H_rows[] = { +22, 18, 15, 63, 16, 13, 1, 2, 29, 25, 28, 4, 36, 10, 38, 7, 60, 23, 11, 38, 28, 1, 12, 31, 57, 45, 57, 30, 23, 59, 67, 14, 16, 4, 14, 62, 15, 50, 7, 70, 64, 6, 42, 48, 9, 31, 19, 40, 49, 2, 25, 3, 41, 49, 36, 9, 29, 39, 31, 5, 17, 1, 29, 25, 11, 21, 18, 2, 8, 22, 39, 15, 8, 22, 13, 3, 19, 4, 21, 62, 34, 43, 6, 24, 17, 60, 8, 74, 6, 44, 60, 10, 33, 12, 26, 24, 45, 81, 69, 80, 41, 28, 23, 5, 10, 20, 52, 18, 13, 86, 3, 7, 59, 21, 65, 72, 34, 37, 26, 55, 47, 48, 34, 5, 44, 47, 68, 96, 82, 111, 61, 74, 30, 17, 55, 98, 81, 66, 89, 35, 74, 82, 91, 51, 55, 51, 30, 89, 61, 75, 40, 71, 73, 11, 56, 54, 19, 47, 94, 69, 64, 20, 64, 12, 54, 77, 42, 88, 36, 52, 90, 63, 70, 27, 32, 73, 91, 32, 56, 46, 9, 78, 51, 68, 88, 67, 20, 43, 40, 14, 66, 86, 39, 97, 38, 27, 50, 84, 54, 92, 61, 46, 67, 24, 58, 35, 58, 37, 98, 85, 73, 84, 48, 35, 57, 16, 26, 37, 65, 32, 72, 95, 107, 33, 77, 33, 85, 105, 106, 75, 56, 71, 79, 59, 52, 105, 79, 90, 93, 100, 88, 112, 86, 80, 65, 42, 106, 100, 93, 94, 99, 97, 93, 101, 111, 99, 83, 53, 85, 95, 108, 107, 41, 109, 84, 78, 104, 101, 69, 110, 98, 103, 80, 83, 77, 71, 76, 78, 87, 102, 104, 95, 96, 83, 87, 50, 110, 103, 112, 45, 58, 70, 94, 91, 89, 81, 101, 82, 63, 72, 100, 97, 76, 112, 53, 105, 49, 75, 109, 102, 66, 111, 68, 87, 92, 79, 96, 43, 90, 44, 110, 99, 102, 92, 103, 106, 62, 53, 27, 46, 108, 104, 107, 108, 109, 76 +}; + +const uint16_t HRA_112_112_H_cols[] = { +7, 8, 52, 12, 12, 42, 16, 69, 45, 14, 19, 23, 6, 32, 3, 5, 22, 2, 45, 50, 2, 1, 18, 84, 10, 7, 62, 11, 9, 21, 24, 63, 2, 5, 28, 13, 6, 15, 58, 39, 39, 22, 76, 13, 26, 68, 9, 10, 49, 38, 32, 11, 34, 44, 8, 7, 25, 67, 1, 17, 19, 36, 4, 41, 3, 26, 31, 15, 45, 40, 8, 4, 41, 20, 6, 53, 1, 42, 9, 20, 25, 17, 33, 41, 3, 19, 55, 17, 27, 14, 31, 88, 15, 26, 36, 16, 28, 24, 27, 16, 30, 56, 48, 43, 4, 5, 38, 37, 40, 46, 18, 18, 22, 50, 76, 34, 60, 83, 39, 73, 56, 92, 42, 52, 75, 35, 37, 33, 61, 67, 47, 75, 66, 70, 29, 92, 51, 95, 84, 21, 57, 28, 46, 66, 93, 11, 94, 55, 96, 20, 71, 48, 53, 43, 82, 90, 66, 90, 14, 44, 54, 62, 34, 58, 81, 53, 23, 43, 27, 93, 10, 86, 37, 80, 60, 49, 21, 79, 74, 72, 48, 61, 40, 76, 64, 29, 38, 79, 51, 54, 13, 49, 72, 30, 50, 86, 35, 80, 61, 56, 36, 59, 65, 91, 25, 47, 58, 59, 78, 47, 32, 24, 44, 86, 64, 57, 12, 23, 109, 107, 85, 63, 31, 65, 62, 68, 111, 78, 104, 89, 112, 87, 69, 105, 65, 94, 109, 78, 72, 104, 85, 108, 77, 106, 79, 74, 103, 96, 64, 105, 105, 102, 63, 35, 59, 108, 112, 81, 102, 57, 106, 83, 81, 77, 101, 55, 94, 96, 97, 106, 46, 101, 83, 85, 71, 107, 104, 87, 33, 67, 103, 95, 30, 91, 89, 103, 75, 51, 107, 87, 91, 89, 99, 68, 52, 109, 99, 88, 84, 112, 54, 70, 92, 100, 98, 74, 60, 100, 98, 110, 90, 73, 71, 95, 70, 100, 29, 69, 110, 93, 82, 97, 98, 77, 73, 99, 101, 108, 82, 102, 111, 110, 111, 97, 88, 80 +}; +const float HRA_112_112_input[] = { +-3.7496794787890972, 14.372112019392226, -7.5640452729302359, 6.9426063455159657, 5.3103644888713299, -6.9203550501252273, 8.4296575778653775, 13.495087143587781, 18.111520666852243, -9.9125748623510912, 10.601298534930972, -10.468591112149715, -9.0757329437720475, -14.471433733514324, 5.2048820572852641, -11.353785810284556, -9.4511008284496416, -9.5255219979484025, -2.0499245561876696, -9.8739646459388748, 22.03442141444015, -9.9745566449839878, -8.4276711655946226, -4.9811962116476307, -13.018434575859896, -5.3358535334627293, -5.6704294937789648, 14.243964608060018, -11.417925510314507, 9.1332657371467878, -14.380214782394296, 14.090409878618974, 6.5602278279998272, 15.53025696352436, -9.1752771765906616, -11.384503450560766, 12.240329442222599, -12.640059450058276, -11.824715154614376, -13.487656131954735, 15.38073452845444, -13.816294924566529, 6.3461114450644454, -2.5192445130977559, -11.916088712873863, 5.4360722876642518, 0.038031547223147381, -12.367220238860654, -2.747864039796549, -14.920508782249289, 16.487336720060863, -13.290002442259247, 19.142698450560925, -0.39443060583296108, 11.723442316413736, -3.6131702833965047, -4.6196487103817017, -11.794290650694531, -14.342351103186955, 2.8079943208330334, -15.290175151123936, 9.0801740558512414, 10.184385069676226, 8.400722260237572, 9.3504690108712936, -14.223531676384166, 11.752768386971752, 11.36995822251677, -15.285021241405444, -13.070613695054403, -11.869191325617697, 4.3191750845563401, 2.0836933404582791, -16.363829786416495, -5.7778094839806595, 11.06389861779129, 13.285433846434705, 9.2552396418849021, -11.065999403824057, -10.167040394420443, -7.0107225044503565, 2.3886881673282474, 5.0014484787306932, -9.2464083853314278, -12.043309174487364, -11.638411967211738, -16.302815497922911, 13.347129717938067, -4.1390259986125226, 0.7947480277507295, 11.538620744796759, -7.4410706619926028, 14.572449028311253, 12.392747919231169, -3.3027890746379289, -9.8431096813736687, 11.582657487369399, -7.85736442083219, -7.3780721969188443, -7.4006260265172212, -8.3937994980934327, -6.6804071011469555, 19.656301355404196, 1.1084340389939762, 3.6028635453146465, -4.5409495140900562, 7.3831459854578982, -5.5905999874445662, -13.852328482738232, 8.9999210644983041, 8.4742375282492315, 16.989947243749878, 7.5590035165610168, -6.154674423116183, 4.1119120658251855, 12.351217703790844, 11.070972687846792, 11.182587746846833, -4.9345619923565645, 9.0054892370887334, -10.841725474869696, 13.902796293412067, -6.7575171884905396, -5.8196703210757335, 1.9284357540668857, 9.9905382141440455, -13.983067199220674, -4.9130522479706453, -8.2369300184767908, 6.8953565265629644, 2.9285103862640871, -2.6303471135655325, -8.3563361642086047, 9.5712349244763715, 4.9728623009661161, -11.045088919587242, -5.7781337596219604, -17.732999074602972, 8.1353860976076646, -11.066240843831284, -1.7079574457159534, -16.411685365171998, -9.0471090651358299, -10.959376227315447, 8.5840398495674126, 6.6373658260736024, 11.422094029020409, 14.85785089306844, 13.185747281780415, 4.2063935223916191, -6.9166135608899282, 10.843153262137262, 5.3913075109409441, -10.744469667642237, -12.491640291445655, 14.141118162062066, 16.425476099516025, 9.8833761863476042, 2.8719064151687883, 14.982021915112442, 1.3588165304065343, -11.657839635726177, 11.066314862965077, -3.0565490195476204, 1.7820159270701772, -13.535333311782074, 4.4026933190218367, -11.097334550496313, -11.322820869044248, 16.418516996530371, 5.8239202459876136, 15.054905601216154, -9.3058742038490152, 8.48902767802557, -8.3853534273227748, -7.9255089736435176, -9.6156735881618811, 11.502594413898008, -6.0542015398269911, 7.1229229147355149, 0.31483632310264387, -11.482093481730768, -9.3225703551629309, 5.8001228713062831, -9.3515917458791051, 7.9778737065172969, 9.7095180444854847, -14.060064536791135, 4.9797253221020545, -6.9210799657794224, 6.6736460552213845, -7.7636429824024606, 10.233132490278882, 8.401747393605044, -10.861100567451366, 13.631509744686715, -15.723791754613185, -8.7931294115815923, -9.9520037489001609, -10.312792052906007, -8.0681893911111917, -15.411052087079765, 10.938779471602952, -8.751795633239853, -9.1302029882284419, -2.3357314769649777, -7.9130658335895596, 2.7508172894969509, -9.1666780515772324, 12.793063537524359, -13.39091818112591, 7.2827402370664842, -10.400778532411657, -1.90854156128735, -4.1272702472088971, 12.696932922959466, -4.0180403457213805, 10.828999052972396, 14.720617452742685, -8.3763729074389719, 3.955093172344033, 0.90932711822659873, -5.6696817865337819, -5.8822086115513805 +}; +const char HRA_112_112_detected_data[] = { +1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1 +}; + +} // FreeDV diff --git a/libfreedv/HRA_112_112.h b/libfreedv/HRA_112_112.h new file mode 100644 index 000000000..bc408a2ef --- /dev/null +++ b/libfreedv/HRA_112_112.h @@ -0,0 +1,23 @@ +/* + FILE....: HRA_112_112.h + + Static arrays for LDPC codec HRA_112_112, generated by ldpc_gen_c_h_file.m. +*/ + +#define HRA_112_112_NUMBERPARITYBITS 112 +#define HRA_112_112_MAX_ROW_WEIGHT 3 +#define HRA_112_112_CODELENGTH 224 +#define HRA_112_112_NUMBERROWSHCOLS 112 +#define HRA_112_112_MAX_COL_WEIGHT 3 +#define HRA_112_112_DEC_TYPE 0 +#define HRA_112_112_MAX_ITER 100 + +namespace FreeDV +{ + +extern const uint16_t HRA_112_112_H_rows[]; +extern const uint16_t HRA_112_112_H_cols[]; +extern const float HRA_112_112_input[]; +extern const char HRA_112_112_detected_data[]; + +} // FreeDV diff --git a/libfreedv/gp_interleaver.cpp b/libfreedv/gp_interleaver.cpp new file mode 100644 index 000000000..5a34b3b82 --- /dev/null +++ b/libfreedv/gp_interleaver.cpp @@ -0,0 +1,122 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: gp_interleaver.c + AUTHOR......: David Rowe + DATE CREATED: April 2018 + + Golden Prime Interleaver. My interpretation of "On the Analysis and + Design of Good Algebraic Interleavers", Xie et al,eq (5). + + See also octvae/gp_interleaver.m + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2018 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, as + published by the Free Software Foundation. This program is + distributed in the hope that it will be useful, but WITHOUT ANY + WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, see . +*/ + +#include +#include +#include "gp_interleaver.h" + +namespace FreeDV +{ + +/* + Choose b for Golden Prime Interleaver. b is chosen to be the + closest integer, which is relatively prime to N, to the Golden + section of N. + + Implemented with a LUT in C for convenience, Octave version + has a more complete implementation. +*/ + +int b_table[] = { + 112,71, + 224,139, + 448,277, + 672,419, + 896,557, + 1120,701, + 1344,839, + 1568,971, + 1792,1109, + 2016,1249, + 2240,1399, + 2464,1523, + 2688,1663, + 2912,1801, + 3136,1949, + 3360,2081, + 3584,2213 +}; + +int choose_interleaver_b(int Nbits) +{ + unsigned int i; + + for(i=0; i. + */ + +#include +#include +#include +#include +#include +#include + +#include "interldpc.h" +#include "codec2_ofdm.h" +#include "mpdecode_core.h" +#include "gp_interleaver.h" +#include "HRA_112_112.h" + +namespace FreeDV +{ + +/* CRC type function, used to compare QPSK vectors when debugging */ + +COMP test_acc(COMP v[], int n) { + COMP acc = {0.0f, 0.0f}; + int i; + + for (i = 0; i < n; i++) { + acc.real += roundf(v[i].real); + acc.imag += roundf(v[i].imag); + } + + return acc; +} + +void printf_n(COMP v[], int n) { + int i; + + for (i = 0; i < n; i++) { + fprintf(stderr, "%d %10f %10f\n", i, round(v[i].real), round(v[i].imag)); + } +} + +void set_up_hra_112_112(struct LDPC *ldpc, struct OFDM_CONFIG *config) { + ldpc->max_iter = HRA_112_112_MAX_ITER; + ldpc->dec_type = 0; + ldpc->q_scale_factor = 1; + ldpc->r_scale_factor = 1; + ldpc->CodeLength = HRA_112_112_CODELENGTH; + ldpc->NumberParityBits = HRA_112_112_NUMBERPARITYBITS; + ldpc->NumberRowsHcols = HRA_112_112_NUMBERROWSHCOLS; + ldpc->max_row_weight = HRA_112_112_MAX_ROW_WEIGHT; + ldpc->max_col_weight = HRA_112_112_MAX_COL_WEIGHT; + ldpc->H_rows = (uint16_t *) HRA_112_112_H_rows; + ldpc->H_cols = (uint16_t *) HRA_112_112_H_cols; + + /* provided for convenience and to match Octave vaiable names */ + + ldpc->data_bits_per_frame = HRA_112_112_CODELENGTH - HRA_112_112_NUMBERPARITYBITS; + ldpc->coded_bits_per_frame = HRA_112_112_CODELENGTH; + ldpc->coded_syms_per_frame = ldpc->coded_bits_per_frame / config->bps; +} + +void ldpc_encode_frame(struct LDPC *ldpc, int codeword[], unsigned char tx_bits_char[]) { + unsigned char *pbits = new unsigned char[ldpc->NumberParityBits]; + int i, j; + + encode(ldpc, tx_bits_char, pbits); + + for (i = 0; i < ldpc->data_bits_per_frame; i++) { + codeword[i] = tx_bits_char[i]; + } + + for (j = 0; i < ldpc->coded_bits_per_frame; i++, j++) { + codeword[i] = pbits[j]; + } + + delete[] pbits; +} + +void qpsk_modulate_frame(COMP tx_symbols[], int codeword[], int n) { + int s, i; + int dibit[2]; + std::complex qpsk_symb; + + for (s = 0, i = 0; i < n; s += 2, i++) { + dibit[0] = codeword[s + 1] & 0x1; + dibit[1] = codeword[s] & 0x1; + qpsk_symb = qpsk_mod(dibit); + tx_symbols[i].real = std::real(qpsk_symb); + tx_symbols[i].imag = std::imag(qpsk_symb); + } +} + +void interleaver_sync_state_machine(struct OFDM *ofdm, + struct LDPC *ldpc, + struct OFDM_CONFIG *config, + COMP codeword_symbols_de[], + float codeword_amps_de[], + float EsNo, int interleave_frames, + int *iter, int *parityCheckCount, int *Nerrs_coded) +{ + (void) config; + int coded_syms_per_frame = ldpc->coded_syms_per_frame; + int coded_bits_per_frame = ldpc->coded_bits_per_frame; + int data_bits_per_frame = ldpc->data_bits_per_frame; + float *llr = new float[coded_bits_per_frame]; + uint8_t *out_char = new uint8_t[coded_bits_per_frame]; + State next_sync_state_interleaver; + + next_sync_state_interleaver = ofdm->sync_state_interleaver; + + if ((ofdm->sync_state_interleaver == search) && (ofdm->frame_count >= (interleave_frames - 1))) { + symbols_to_llrs(llr, codeword_symbols_de, codeword_amps_de, EsNo, ofdm->mean_amp, coded_syms_per_frame); + iter[0] = run_ldpc_decoder(ldpc, out_char, llr, parityCheckCount); + Nerrs_coded[0] = data_bits_per_frame - parityCheckCount[0]; + + if ((Nerrs_coded[0] == 0) || (interleave_frames == 1)) { + /* sucessful decode! */ + next_sync_state_interleaver = synced; + ofdm->frame_count_interleaver = interleave_frames; + } + } + + ofdm->sync_state_interleaver = next_sync_state_interleaver; + delete[] out_char; + delete[] llr; +} + +/* measure uncoded (raw) bit errors over interleaver frame, note we + don't include txt bits as this is done after we dissassemmble the + frame */ + +int count_uncoded_errors(struct LDPC *ldpc, struct OFDM_CONFIG *config, int Nerrs_raw[], int interleave_frames, COMP codeword_symbols_de[]) { + int i, j, Nerrs, Terrs; + + int coded_syms_per_frame = ldpc->coded_syms_per_frame; + int coded_bits_per_frame = ldpc->coded_bits_per_frame; + int data_bits_per_frame = ldpc->data_bits_per_frame; + int *rx_bits_raw = new int[coded_bits_per_frame]; + + /* generate test codeword from known payload data bits */ + + int *test_codeword = new int[coded_bits_per_frame]; + uint16_t *r = new uint16_t[data_bits_per_frame]; + uint8_t *tx_bits = new uint8_t[data_bits_per_frame]; + + ofdm_rand(r, data_bits_per_frame); + + for (i = 0; i < data_bits_per_frame; i++) { + tx_bits[i] = r[i] > 16384; + } + + ldpc_encode_frame(ldpc, test_codeword, tx_bits); + + Terrs = 0; + for (j = 0; j < interleave_frames; j++) { + for (i = 0; i < coded_syms_per_frame; i++) { + int bits[2]; + std::complex s = std::complex{codeword_symbols_de[j * coded_syms_per_frame + i].real, codeword_symbols_de[j * coded_syms_per_frame + i].imag}; + qpsk_demod(s, bits); + rx_bits_raw[config->bps * i] = bits[1]; + rx_bits_raw[config->bps * i + 1] = bits[0]; + } + + Nerrs = 0; + + for (i = 0; i < coded_bits_per_frame; i++) { + if (test_codeword[i] != rx_bits_raw[i]) { + Nerrs++; + } + } + + Nerrs_raw[j] = Nerrs; + Terrs += Nerrs; + } + + delete[] tx_bits; + delete[] r; + delete[] test_codeword; + delete[] rx_bits_raw; + + return Terrs; +} + +int count_errors(uint8_t tx_bits[], uint8_t rx_bits[], int n) { + int i; + int Nerrs = 0; + + for (i = 0; i < n; i++) { + if (tx_bits[i] != rx_bits[i]) { + Nerrs++; + } + } + + return Nerrs; +} + +/* + Given an array of tx_bits, LDPC encodes, interleaves, and OFDM + modulates. + + Note this could be refactored to save memory, e.g. for embedded + applications we could call ofdm_txframe on a frame by frame + basis + */ + +void ofdm_ldpc_interleave_tx(struct OFDM *ofdm, struct LDPC *ldpc, std::complex tx_sams[], uint8_t tx_bits[], uint8_t txt_bits[], int interleave_frames, struct OFDM_CONFIG *config) { + int coded_syms_per_frame = ldpc->coded_syms_per_frame; + int coded_bits_per_frame = ldpc->coded_bits_per_frame; + int data_bits_per_frame = ldpc->data_bits_per_frame; + int ofdm_bitsperframe = ofdm_get_bits_per_frame(); + int *codeword = new int[coded_bits_per_frame]; + COMP *coded_symbols = new COMP[interleave_frames * coded_syms_per_frame]; + COMP *coded_symbols_inter = new COMP[interleave_frames * coded_syms_per_frame]; + int Nsamperframe = ofdm_get_samples_per_frame(); + std::complex *tx_symbols = new std::complex[ofdm_bitsperframe / config->bps]; + int j; + + for (j = 0; j < interleave_frames; j++) { + ldpc_encode_frame(ldpc, codeword, &tx_bits[j * data_bits_per_frame]); + qpsk_modulate_frame(&coded_symbols[j * coded_syms_per_frame], codeword, coded_syms_per_frame); + } + + gp_interleave_comp(coded_symbols_inter, coded_symbols, interleave_frames * coded_syms_per_frame); + + for (j = 0; j < interleave_frames; j++) { + ofdm_assemble_modem_frame_symbols(tx_symbols, &coded_symbols_inter[j * coded_syms_per_frame], &txt_bits[config->txtbits * j]); + ofdm_txframe(ofdm, &tx_sams[j * Nsamperframe], tx_symbols); + } + + delete[] tx_symbols; + delete[] coded_symbols_inter; + delete[] coded_symbols; + delete[] codeword; +} + +} // FreeDV diff --git a/libfreedv/mpdecode_core.cpp b/libfreedv/mpdecode_core.cpp new file mode 100644 index 000000000..50d0d0277 --- /dev/null +++ b/libfreedv/mpdecode_core.cpp @@ -0,0 +1,726 @@ +/* + FILE...: mpdecode_core.c + AUTHOR.: Matthew C. Valenti, Rohit Iyer Seshadri, David Rowe + CREATED: Sep 2016 + + C-callable core functions moved from MpDecode.c, so they can be used for + Octave and C programs. +*/ + +#include +#include +#include +#include +#include +#include "mpdecode_core.h" +#ifndef USE_ORIGINAL_PHI0 +#include "phi0.h" +#endif + +#ifdef __EMBEDDED__ +#include "machdep.h" +#endif + +#define QPSK_CONSTELLATION_SIZE 4 +#define QPSK_BITS_PER_SYMBOL 2 + +namespace FreeDV +{ + +/* QPSK constellation for symbol likelihood calculations */ + +static COMP S_matrix[] = { + { 1.0f, 0.0f}, + { 0.0f, 1.0f}, + { 0.0f, -1.0f}, + {-1.0f, 0.0f} +}; + +// c_nodes will be an array of NumberParityBits of struct c_node +// Each c_node contains an array of c_sub_node elements +// This structure reduces the indexing caluclations in SumProduct() + +struct c_sub_node { // Order is important here to keep total size small. + uint16_t index; // Values from H_rows (except last 2 entries) + uint16_t socket; // The socket number at the v_node + float message; // modified during operation! +}; + +struct c_node { + int degree; // A count of elements in the following arrays + struct c_sub_node *subs; +}; + +// v_nodes will be an array of CodeLength of struct v_node + +struct v_sub_node { + uint16_t index; // the index of a c_node it is connected to + // Filled with values from H_cols (except last 2 entries) + uint16_t socket; // socket number at the c_node + float message; // Loaded with input data + // modified during operation! + uint8_t sign; // 1 if input is negative + // modified during operation! +}; + +struct v_node { + int degree; // A count of ??? + float initial_value; + struct v_sub_node *subs; +}; + +void encode(struct LDPC *ldpc, unsigned char ibits[], unsigned char pbits[]) { + unsigned int tmp, par, prev=0; + int i, p, ind; + uint16_t *H_rows = ldpc->H_rows; + + for (p=0; pNumberParityBits; p++) { + par = 0; + + for (i=0; imax_row_weight; i++) { + ind = H_rows[p + i*ldpc->NumberParityBits]; + par = par + ibits[ind-1]; + } + + tmp = par + prev; + + tmp &= 1; // only retain the lsb + prev = tmp; + pbits[p] = tmp; + } +} + +#ifdef USE_ORIGINAL_PHI0 +/* Phi function */ +static float phi0( + float x ) +{ + float z; + + if (x>10) + return( 0 ); + else if (x< 9.08e-5 ) + return( 10 ); + else if (x > 9) + return( 1.6881e-4 ); + /* return( 1.4970e-004 ); */ + else if (x > 8) + return( 4.5887e-4 ); + /* return( 4.0694e-004 ); */ + else if (x > 7) + return( 1.2473e-3 ); + /* return( 1.1062e-003 ); */ + else if (x > 6) + return( 3.3906e-3 ); + /* return( 3.0069e-003 ); */ + else if (x > 5) + return( 9.2168e-3 ); + /* return( 8.1736e-003 ); */ + else { + z = (float) exp(x); + return( (float) log( (z+1)/(z-1) ) ); + } +} +#endif + + +/* Values for linear approximation (DecoderType=5) */ + +#define AJIAN -0.24904163195436 +#define TJIAN 2.50681740420944 + +/* The linear-log-MAP algorithm */ + +static float max_star0( + float delta1, + float delta2 ) +{ + register float diff; + + diff = delta2 - delta1; + + if ( diff > TJIAN ) + return( delta2 ); + else if ( diff < -TJIAN ) + return( delta1 ); + else if ( diff > 0 ) + return( delta2 + AJIAN*(diff-TJIAN) ); + else + return( delta1 - AJIAN*(diff+TJIAN) ); +} + +void init_c_v_nodes(struct c_node *c_nodes, + int shift, + int NumberParityBits, + int max_row_weight, + uint16_t *H_rows, + int H1, + int CodeLength, + struct v_node *v_nodes, + int NumberRowsHcols, + uint16_t *H_cols, + int max_col_weight, + int dec_type, + float *input) +{ + int i, j, k, count, cnt, c_index, v_index; + + /* first determine the degree of each c-node */ + + if (shift ==0){ + for (i=0;i 0 ) { + count++; + } + } + c_nodes[i].degree = count; + if (H1){ + if (i==0){ + c_nodes[i].degree=count+1; + } + else{ + c_nodes[i].degree=count+2; + } + } + } + } + else{ + cnt=0; + for (i=0;i<(NumberParityBits/shift);i++) { + for (k=0;k 0 ) { + count++; + } + } + c_nodes[cnt].degree = count; + if ((i==0)||(i==(NumberParityBits/shift)-1)){ + c_nodes[cnt].degree=count+1; + } + else{ + c_nodes[cnt].degree=count+2; + } + cnt++; + } + } + } + + if (H1){ + + if (shift ==0){ + for (i=0;i0){ + cnt=0; + for (i=0;i<(NumberParityBits/shift);i++){ + + for (k =0;k 0 ) { + count++; + } + } + v_nodes[i].degree = count; + } + + for(i=CodeLength-NumberParityBits+shift;i 0 ) { + count++; + } + } + v_nodes[i].degree = count; + } + } + + if (shift>0){ + v_nodes[CodeLength-1].degree =v_nodes[CodeLength-1].degree+1; + } + + + /* set up v_nodes */ + + for (i=0;i=CodeLength-NumberParityBits+shift)){ + v_nodes[i].subs[j].index=i-(CodeLength-NumberParityBits+shift)+count; + if (shift ==0){ + count=count+1; + } + else{ + count=count+shift; + } + } else { + v_nodes[i].subs[j].index = (H_cols[i+j*NumberRowsHcols] - 1); + } + + /* search the connected c-node for the proper message value */ + for (c_index=0;c_indexindex ].subs[ cp->socket ]; + phi_sum += vp->message; + sign ^= vp->sign; + } + + if (sign==0) ssum++; + + for (i=0;iindex ].subs[ cp->socket ]; + if ( sign ^ vp->sign ) { + cp->message = -phi0( phi_sum - vp->message ); // *r_scale_factor; + } else + cp->message = phi0( phi_sum - vp->message ); // *r_scale_factor; + } + } + + /* update q */ + for (i=0;iindex ].subs[ vp->socket ].message; + } + + /* make hard decision */ + if (Qi < 0) { + DecodedBits[i] = 1; + } + + /* now subtract to get the extrinsic information */ + for (j=0;jindex ].subs[ vp->socket ].message; + + vp->message = phi0( fabs( temp_sum ) ); // *q_scale_factor; + if (temp_sum > 0) + vp->sign = 0; + else + vp->sign = 1; + } + } + + /* count data bit errors, assuming that it is systematic */ + for (i=0;imax_iter; + dec_type = ldpc->dec_type; + q_scale_factor = ldpc->q_scale_factor; + r_scale_factor = ldpc->r_scale_factor; + + CodeLength = ldpc->CodeLength; /* length of entire codeword */ + NumberParityBits = ldpc->NumberParityBits; + NumberRowsHcols = ldpc->NumberRowsHcols; + + char *DecodedBits = (char*) calloc( CodeLength, sizeof( char ) ); + assert(DecodedBits); + + /* derive some parameters */ + + shift = (NumberParityBits + NumberRowsHcols) - CodeLength; + if (NumberRowsHcols == CodeLength) { + H1=0; + shift=0; + } else { + H1=1; + } + + max_row_weight = ldpc->max_row_weight; + max_col_weight = ldpc->max_col_weight; + + /* initialize c-node and v-node structures */ + + c_nodes = (struct c_node*) calloc( NumberParityBits, sizeof( struct c_node ) ); + assert(c_nodes); + v_nodes = (struct v_node*) calloc( CodeLength, sizeof( struct v_node)); + assert(v_nodes); + + init_c_v_nodes(c_nodes, shift, NumberParityBits, max_row_weight, ldpc->H_rows, H1, CodeLength, + v_nodes, NumberRowsHcols, ldpc->H_cols, max_col_weight, dec_type, input); + + int DataLength = CodeLength - NumberParityBits; + int *data_int = (int*) calloc( DataLength, sizeof(int) ); + + /* need to clear these on each call */ + + for(i=0; i 0.0L) - (sd[i] < 0.0L); + x = (sd[i]/mean - sign); + sum += x; + sumsq += x*x; + } + estvar = (n * sumsq - sum * sum) / (n * (n - 1)); + //fprintf(stderr, "mean: %f var: %f\n", mean, estvar); + + estEsN0 = 1.0/(2.0L * estvar + 1E-3); + for(i=0; i> 1; + } + mask = 1 << (bps - 1); + + for (k=0;k> 1; + } + } + for (k=0;kmax_iter = %d\n", ldpc->max_iter); +fprintf(stderr, "ldpc->dec_type = %d\n", ldpc->dec_type); +fprintf(stderr, "ldpc->q_scale_factor = %d\n", ldpc->q_scale_factor); +fprintf(stderr, "ldpc->r_scale_factor = %d\n", ldpc->r_scale_factor); +fprintf(stderr, "ldpc->CodeLength = %d\n", ldpc->CodeLength); +fprintf(stderr, "ldpc->NumberParityBits = %d\n", ldpc->NumberParityBits); +fprintf(stderr, "ldpc->NumberRowsHcols = %d\n", ldpc->NumberRowsHcols); +fprintf(stderr, "ldpc->max_row_weight = %d\n", ldpc->max_row_weight); +fprintf(stderr, "ldpc->max_col_weight = %d\n", ldpc->max_col_weight); +fprintf(stderr, "ldpc->data_bits_per_frame = %d\n", ldpc->data_bits_per_frame); +fprintf(stderr, "ldpc->coded_bits_per_frame = %d\n", ldpc->coded_bits_per_frame); +fprintf(stderr, "ldpc->coded_syms_per_frame = %d\n", ldpc->coded_syms_per_frame); +} + +} // FreeDV + +/* vi:set ts=4 et sts=4: */ diff --git a/libfreedv/ofdm.cpp b/libfreedv/ofdm.cpp new file mode 100644 index 000000000..102c36b8b --- /dev/null +++ b/libfreedv/ofdm.cpp @@ -0,0 +1,1752 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: ofdm.c + AUTHORS.....: David Rowe & Steve Sampson + DATE CREATED: June 2017 + + A Library of functions that implement a QPSK OFDM modem, C port of + the Octave functions in ofdm_lib.m + +\*---------------------------------------------------------------------------*/ +/* + Copyright (C) 2017/2018 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, as + published by the Free Software Foundation. This program is + distributed in the hope that it will be useful, but WITHOUT ANY + WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public + License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, see . + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "codec2/comp.h" +#include "ofdm_internal.h" +#include "codec2_ofdm.h" +#include "freedv_filter.h" + +namespace FreeDV +{ + +/* Static Prototypes */ + +static std::complex vector_sum(std::complex *, int); +static void dft(struct OFDM *, std::complex *, std::complex *); +static void idft(struct OFDM *, std::complex *, std::complex *); +static void ofdm_demod_core(struct OFDM *ofdm, int *rx_bits); +static int ofdm_sync_search_core(struct OFDM *ofdm); + +/* Defines */ + +#define max( a, b ) ( ((a) > (b)) ? (a) : (b) ) +#define min( a, b ) ( ((a) < (b)) ? (a) : (b) ) + +/* + * QPSK Quadrant bit-pair values - Gray Coded + */ +static const std::complex constellation[] = { + std::complex{1.0f, 0.0f}, + std::complex{0.0f, 1.0f}, + std::complex{0.0f, - 1.0f}, + std::complex{-1.0f, + 0.0f} +}; + +/* + * These pilots are compatible with Octave version + */ +static const int8_t pilotvalues[] = { + -1,-1, 1, 1,-1,-1,-1, 1, + -1, 1,-1, 1, 1, 1, 1, 1, + 1, 1, 1,-1,-1, 1,-1, 1, + -1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1,-1, 1, 1, 1, 1, + 1,-1,-1,-1,-1,-1,-1, 1, + -1, 1,-1, 1,-1,-1, 1,-1, + 1, 1, 1, 1,-1, 1,-1, 1 +}; + +static std::complex *tx_uw_syms; +static int *uw_ind; +static int *uw_ind_sym; + +/* static variables */ + +static struct OFDM_CONFIG ofdm_config; + +static float ofdm_tx_centre; /* TX Center frequency */ +static float ofdm_rx_centre; /* RX Center frequency */ +static float ofdm_fs; /* Sample rate */ +static float ofdm_ts; /* Symbol cycle time */ +static float ofdm_rs; /* Symbol rate */ +static float ofdm_tcp; /* Cyclic prefix duration */ +static float ofdm_inv_m; /* 1/m */ +static float ofdm_tx_nlower; /* TX lowest carrier freq */ +static float ofdm_rx_nlower; /* RX lowest carrier freq */ +static float ofdm_doc; /* division of radian circle */ +static float ofdm_timing_mx_thresh; /* See 700D Part 4 Acquisition blog post and ofdm_dev.m routines for how this was set */ + +static int ofdm_nc; /* NS-1 data symbols between pilots */ +static int ofdm_ns; +static int ofdm_bps; /* Bits per symbol */ +static int ofdm_m; /* duration of each symbol in samples */ +static int ofdm_ncp; /* duration of CP in samples */ + +static int ofdm_ftwindowwidth; +static int ofdm_bitsperframe; +static int ofdm_rowsperframe; +static int ofdm_samplesperframe; +static int ofdm_max_samplesperframe; +static int ofdm_rxbuf; +static int ofdm_ntxtbits; /* reserve bits/frame for auxillary text information */ +static int ofdm_nuwbits; /* Unique word, used for positive indication of lock */ + +/* Functions -------------------------------------------------------------------*/ + +static float cnormf(std::complex val) { + float realf = val.real(); + float imagf = val.imag(); + + return realf * realf + imagf * imagf; +} + +/* Gray coded QPSK modulation function */ + +std::complex qpsk_mod(int *bits) { + return constellation[(bits[1] << 1) | bits[0]]; +} + +/* Gray coded QPSK demodulation function */ + +void qpsk_demod(std::complex symbol, int *bits) { + std::complex rotate = symbol * cmplx(ROT45); + + bits[0] = rotate.real() < 0.0f; + bits[1] = rotate.imag() < 0.0f; +} + +/* + * ------------ + * ofdm_create + * ------------ + * + * Returns OFDM data structure on success + * Return NULL on fail + * + * If you want the defaults, call this with config structure + * and the NC setting to 0. This will fill the structure with + * default values of the original OFDM modem. + */ + +struct OFDM *ofdm_create(const struct OFDM_CONFIG *config) { + struct OFDM *ofdm; + float tval; + int i, j; + + /* Check if called correctly */ + + if (config == NULL) { + return NULL; + } + + if (config->nc == 0) { + /* Fill in default values */ + + ofdm_nc = 17; /* Number of carriers */ + ofdm_ns = 8; /* Number of Symbol frames */ + ofdm_bps = 2; /* Bits per Symbol */ + ofdm_ts = 0.018f; + ofdm_tcp = .002f; /* Cyclic Prefix duration */ + ofdm_tx_centre = 1500.0f; /* TX Centre Audio Frequency */ + ofdm_rx_centre = 1500.0f; /* RX Centre Audio Frequency */ + ofdm_fs = 8000.0f; /* Sample Frequency */ + ofdm_ntxtbits = 4; + ofdm_ftwindowwidth = 11; + ofdm_timing_mx_thresh = 0.30f; + } else { + /* Use the users values */ + + ofdm_nc = config->nc; /* Number of carriers */ + ofdm_ns = config->ns; /* Number of Symbol frames */ + ofdm_bps = config->bps; /* Bits per Symbol */ + ofdm_ts = config->ts; + ofdm_tcp = config->tcp; /* Cyclic Prefix duration */ + ofdm_tx_centre = config->tx_centre; /* TX Centre Audio Frequency */ + ofdm_rx_centre = config->rx_centre; /* RX Centre Audio Frequency */ + ofdm_fs = config->fs; /* Sample Frequency */ + ofdm_ntxtbits = config->txtbits; + ofdm_ftwindowwidth = config->ftwindowwidth; + ofdm_timing_mx_thresh = config->ofdm_timing_mx_thresh; + } + + ofdm_rs = (1.0f / ofdm_ts); /* Modulation Symbol Rate */ + ofdm_m = (int) (ofdm_fs / ofdm_rs); /* 144 */ + ofdm_ncp = (int) (ofdm_tcp * ofdm_fs); /* 16 */ + ofdm_inv_m = (1.0f / (float) ofdm_m); + + /* Copy structure into global */ + + ofdm_config.tx_centre = ofdm_tx_centre; + ofdm_config.rx_centre = ofdm_rx_centre; + ofdm_config.fs = ofdm_fs; + ofdm_config.rs = ofdm_rs; + ofdm_config.ts = ofdm_ts; + ofdm_config.tcp = ofdm_tcp; + ofdm_config.ofdm_timing_mx_thresh = ofdm_timing_mx_thresh; + ofdm_config.nc = ofdm_nc; + ofdm_config.ns = ofdm_ns; + ofdm_config.bps = ofdm_bps; + ofdm_config.txtbits = ofdm_ntxtbits; + ofdm_config.ftwindowwidth = ofdm_ftwindowwidth; + + /* Calculate sizes from config param */ + + ofdm_bitsperframe = (ofdm_ns - 1) * (ofdm_nc * ofdm_bps); + ofdm_rowsperframe = ofdm_bitsperframe / (ofdm_nc * ofdm_bps); + ofdm_samplesperframe = ofdm_ns * (ofdm_m + ofdm_ncp); + ofdm_max_samplesperframe = ofdm_samplesperframe + (ofdm_m + ofdm_ncp) / 4; + ofdm_rxbuf = 3 * ofdm_samplesperframe + 3 * (ofdm_m + ofdm_ncp); + ofdm_nuwbits = (ofdm_ns - 1) * ofdm_bps - ofdm_ntxtbits; // 10 + + /* Were ready to start filling in the OFDM structure now */ + + if ((ofdm = (struct OFDM *) malloc(sizeof (struct OFDM))) == NULL) { + return NULL; + } + + ofdm->pilot_samples = (std::complex*) malloc(sizeof (std::complex) * (ofdm_m + ofdm_ncp)); + ofdm->rxbuf = (std::complex*) malloc(sizeof (std::complex) * ofdm_rxbuf); + ofdm->pilots = (std::complex*) malloc(sizeof (std::complex) * (ofdm_nc + 2)); + + /* + * rx_sym is a 2D array of variable size + * + * allocate rx_sym row storage. It is a pointer to a pointer + */ + + ofdm->rx_sym = (std::complex**) malloc(sizeof (std::complex) * (ofdm_ns + 3)); + + /* allocate rx_sym column storage */ + + int free_last_rx_sym = 0; + for (i = 0; i < (ofdm_ns + 3); i++) + { + ofdm->rx_sym[i] = (std::complex *) malloc(sizeof(std::complex) * (ofdm_nc + 2)); + + if (ofdm->rx_sym[i] == NULL) { + free_last_rx_sym = i; + } + + free_last_rx_sym = (ofdm_ns + 3); + } + + /* The rest of these are 1D arrays of variable size */ + + ofdm->rx_np = (std::complex*) malloc(sizeof (std::complex) * (ofdm_rowsperframe * ofdm_nc)); + ofdm->rx_amp = (float*) malloc(sizeof (float) * (ofdm_rowsperframe * ofdm_nc)); + ofdm->aphase_est_pilot_log = (float*) malloc(sizeof (float) * (ofdm_rowsperframe * ofdm_nc)); + ofdm->tx_uw = (uint8_t*) malloc(sizeof (uint8_t) * ofdm_nuwbits); + + for (i = 0; i < ofdm_nuwbits; i++) { + ofdm->tx_uw[i] = 0; + } + + /* Null pointers to unallocated buffers */ + ofdm->ofdm_tx_bpf = NULL; + + /* store complex BPSK pilot symbols */ + + assert(sizeof (pilotvalues) >= (ofdm_nc + 2) * sizeof (int8_t)); + + /* There are only 64 pilot values available */ + + for (i = 0; i < (ofdm_nc + 2); i++) + { + ofdm->pilots[i].real((float) pilotvalues[i]); + ofdm->pilots[i].imag(0.0f); + } + + /* carrier tables for up and down conversion */ + + ofdm_doc = (TAU / (ofdm_fs / ofdm_rs)); + tval = ofdm_rs * ((float) ofdm_nc / 2); + ofdm_tx_nlower = floorf((ofdm_tx_centre - tval) / ofdm_rs); + ofdm_rx_nlower = floorf((ofdm_rx_centre - tval) / ofdm_rs); + + for (i = 0; i < ofdm_rxbuf; i++) { + ofdm->rxbuf[i] = 0.0f; + } + + for (i = 0; i < (ofdm_ns + 3); i++) { + for (j = 0; j < (ofdm_nc + 2); j++) { + ofdm->rx_sym[i][j] = 0.0f; + } + } + + for (i = 0; i < ofdm_rowsperframe * ofdm_nc; i++) { + ofdm->rx_np[i] = 0.0f; + } + + for (i = 0; i < ofdm_rowsperframe; i++) { + for (j = 0; j < ofdm_nc; j++) { + ofdm->aphase_est_pilot_log[ofdm_nc * i + j] = 0.0f; + ofdm->rx_amp[ofdm_nc * i + j] = 0.0f; + } + } + + /* default settings of options and states */ + + ofdm->verbose = 0; + ofdm->timing_en = true; + ofdm->foff_est_en = true; + ofdm->phase_est_en = true; + + ofdm->foff_est_gain = 0.05f; + ofdm->foff_est_hz = 0.0f; + ofdm->sample_point = 0; + ofdm->timing_est = 0; + ofdm->timing_valid = 0; + ofdm->timing_mx = 0.0f; + ofdm->nin = ofdm_samplesperframe; + ofdm->mean_amp = 0.0f; + ofdm->foff_metric = 0.0f; + + /* + * UW symbol placement, designed to get no false syncs at any freq + * offset. Use ofdm_dev.m, debug_false_sync() to test. Note we need + * to pair the UW bits so they fit into symbols. The LDPC decoder + * works on symbols so we can't break up any symbols into UW/LDPC bits. + */ + + uw_ind = (int*) malloc(sizeof (int) * ofdm_nuwbits); + uw_ind_sym = (int*) malloc(sizeof (int) * (ofdm_nuwbits / 2)); + + /* + * The Unique Word is placed in different indexes based on + * the number of carriers requested. + */ + + for (i = 0, j = 0; i < (ofdm_nuwbits / 2); i++, j += 2) { + int val = floorf((i + 1) * (ofdm_nc + 1) / 2); + uw_ind_sym[i] = val; // symbol index + uw_ind[j] = (val * 2); // bit index 1 + uw_ind[j + 1] = (val * 2) + 1; // bit index 2 + } + + tx_uw_syms = (std::complex*) malloc(sizeof (std::complex) * (ofdm_nuwbits / 2)); + + for (i = 0; i < (ofdm_nuwbits / 2); i++) { + tx_uw_syms[i] = 1.0f; // qpsk_mod(0:0) + } + + /* sync state machine */ + + ofdm->sync_state = search; + ofdm->last_sync_state = search; + ofdm->sync_state_interleaver = search; + ofdm->last_sync_state_interleaver = search; + + ofdm->uw_errors = 0; + ofdm->sync_counter = 0; + ofdm->frame_count = 0; + ofdm->sync_start = false; + ofdm->sync_end = false; + ofdm->sync_mode = autosync; + ofdm->frame_count_interleaver = 0; + + /* create the OFDM waveform */ + + std::complex *temp = (std::complex*) malloc(sizeof (std::complex) * ofdm_m); + + idft(ofdm, temp, ofdm->pilots); + + /* + * pilot_samples is 160 samples, but timing and freq offset est + * were found by experiment to work better without a cyclic + * prefix, so we uses zeroes instead. + */ + + /* zero out Cyclic Prefix (CP) values */ + + for (i = 0; i < ofdm_ncp; i++) { + ofdm->pilot_samples[i] = 0.0f; + } + + /* Now copy the whole thing after the above */ + + for (i = ofdm_ncp, j = 0; j < ofdm_m; i++, j++) { + ofdm->pilot_samples[i] = temp[j]; + } + free(temp); + + /* calculate constant used to normalise timing correlation maximum */ + + float acc = 0.0f; + + for (i = 0; i < (ofdm_m + ofdm_ncp); i++) { + acc += cnormf(ofdm->pilot_samples[i]); + } + + ofdm->timing_norm = (ofdm_m + ofdm_ncp) * acc; + ofdm->clock_offset_counter = 0; + ofdm->sig_var = ofdm->noise_var = 1.0f; + ofdm->tx_bpf_en = false; + + return ofdm; /* Success */ + + //// Error return points with free call in the reverse order of allocation: + + free(tx_uw_syms); + free(uw_ind); + free(uw_ind_sym); + free(ofdm->tx_uw); + free(ofdm->aphase_est_pilot_log); + free(ofdm->rx_amp); + free(ofdm->rx_np); + + for (i = 0; i < free_last_rx_sym; i++) { + free(ofdm->rx_sym[i]); + } + + free(ofdm->rx_sym); + free(ofdm->pilots); + free(ofdm->rxbuf); + free(ofdm->pilot_samples); + free(ofdm); + + return(NULL); +} + +void allocate_tx_bpf(struct OFDM *ofdm) { + ofdm->ofdm_tx_bpf = (struct quisk_cfFilter*) malloc(sizeof(struct quisk_cfFilter)); + + /* Transmit bandpass filter; complex coefficients, center frequency */ + + quisk_filt_cfInit(ofdm->ofdm_tx_bpf, filtP550S750, sizeof (filtP550S750) / sizeof (float)); + quisk_cfTune(ofdm->ofdm_tx_bpf, ofdm_tx_centre / ofdm_fs); +} + +void deallocate_tx_bpf(struct OFDM *ofdm) { + quisk_filt_destroy(ofdm->ofdm_tx_bpf); + free(ofdm->ofdm_tx_bpf); + ofdm->ofdm_tx_bpf = NULL; +} + +void ofdm_destroy(struct OFDM *ofdm) { + int i; + + if (ofdm->ofdm_tx_bpf) + deallocate_tx_bpf(ofdm); + + free(ofdm->pilot_samples); + free(ofdm->rxbuf); + free(ofdm->pilots); + + for (i = 0; i < (ofdm_ns + 3); i++) { /* 2D array */ + free(ofdm->rx_sym[i]); + } + + free(ofdm->rx_sym); + free(ofdm->rx_np); + free(ofdm->rx_amp); + free(ofdm->aphase_est_pilot_log); + free(ofdm->tx_uw); + free(ofdm); +} + +/* convert frequency domain into time domain */ + +static void idft(struct OFDM *ofdm, std::complex *result, std::complex *vector) { + (void) ofdm; + int row, col; + + result[0] = 0.0f; + + for (col = 0; col < (ofdm_nc + 2); col++) { + result[0] += vector[col]; // cexp(0j) == 1 + } + + result[0] *= ofdm_inv_m; + + for (row = 1; row < ofdm_m; row++) { + float tval1 = ofdm_tx_nlower * ofdm_doc *row; + float tval2 = ofdm_doc * row; + + std::complex c = cmplx(tval1); + std::complex delta = cmplx(tval2); + + result[row] = 0.0f; + + for (col = 0; col < (ofdm_nc + 2); col++) { + result[row] += (vector[col] * c); + c *= delta; + } + + result[row] *= ofdm_inv_m; + } +} + +/* convert time domain into frequency domain */ + +static void dft(struct OFDM *ofdm, std::complex *result, std::complex *vector) { + (void) ofdm; + int row, col, nlower; + + for (col = 0; col < (ofdm_nc + 2); col++) { + result[col] = vector[0]; // conj(cexp(0j)) == 1 + } + + for (col = 0, nlower = ofdm_rx_nlower; col < (ofdm_nc + 2); col++, nlower++) { + float tval = nlower * ofdm_doc; + std::complex c = cmplxconj(tval); + std::complex delta = c; + + for (row = 1; row < ofdm_m; row++) { + result[col] += (vector[row] * c); + c *= delta; + } + } +} + +static std::complex vector_sum(std::complex *a, int num_elements) { + int i; + + std::complex sum = 0.0f; + + for (i = 0; i < num_elements; i++) { + sum += a[i]; + } + + return sum; +} + +/* + * Correlates the OFDM pilot symbol samples with a window of received + * samples to determine the most likely timing offset. Combines two + * frames pilots so we need at least Nsamperframe+M+Ncp samples in rx. + * + * Can be used for acquisition (coarse timing), and fine timing. + * + * Unlike Octave version use states to return a few values. + */ + +static int est_timing(struct OFDM *ofdm, std::complex *rx, int length) { + std::complex csam; + int Ncorr = length - (ofdm_samplesperframe + (ofdm_m + ofdm_ncp)); + int SFrame = ofdm_samplesperframe; + float *corr = new float[Ncorr]; + int i, j; + + float acc = 0.0f; + + for (i = 0; i < length; i++) { + acc += cnormf(rx[i]); + } + + float av_level = 2.0f * sqrtf(ofdm->timing_norm * acc / length) + 1E-12f; + + for (i = 0; i < Ncorr; i++) { + std::complex corr_st = 0.0f; + std::complex corr_en = 0.0f; + + for (j = 0; j < (ofdm_m + ofdm_ncp); j++) { + csam = std::conj(ofdm->pilot_samples[j]); + + corr_st = corr_st + (rx[i + j ] * csam); + corr_en = corr_en + (rx[i + j + SFrame] * csam); + } + + corr[i] = (std::abs(corr_st) + std::abs(corr_en)) / av_level; + } + + /* find the max magnitude and its index */ + + float timing_mx = 0.0f; + int timing_est = 0; + + for (i = 0; i < Ncorr; i++) { + if (corr[i] > timing_mx) { + timing_mx = corr[i]; + timing_est = i; + } + } + + ofdm->timing_mx = timing_mx; + ofdm->timing_valid = (timing_mx > ofdm_timing_mx_thresh); /* bool but used as external int */ + + if (ofdm->verbose > 1) { + fprintf(stderr, " av_level: %f max: %f timing_est: %d timing_valid: %d\n", (double) av_level, (double) ofdm->timing_mx, timing_est, ofdm->timing_valid); + } + + delete[] corr; + + return timing_est; +} + +/* + * Determines frequency offset at current timing estimate, used for + * coarse freq offset estimation during acquisition. + * + * Freq offset is based on an averaged statistic that was found to be + * necessary to generate good quality estimates. + * + * Keep calling it when in trial or actual sync to keep statistic + * updated, in case we lose sync. + */ + +static float est_freq_offset(struct OFDM *ofdm, std::complex *rx, int timing_est) { + std::complex csam1, csam2; + float foff_est; + int j, k; + + /* + Freq offset can be considered as change in phase over two halves + of pilot symbols. We average this statistic over this and next + frames pilots. + */ + + std::complex p1, p2, p3, p4; + p1 = p2 = p3 = p4 = 0.0f; + + /* calculate phase of pilots at half symbol intervals */ + + for (j = 0, k = (ofdm_m + ofdm_ncp) / 2; j < (ofdm_m + ofdm_ncp) / 2; j++, k++) { + csam1 = std::conj(ofdm->pilot_samples[j]); + csam2 = std::conj(ofdm->pilot_samples[k]); + + /* pilot at start of frame */ + + p1 = p1 + (rx[timing_est + j] * csam1); + p2 = p2 + (rx[timing_est + k] * csam2); + + /* pilot at end of frame */ + + p3 = p3 + (rx[timing_est + j + ofdm_samplesperframe] * csam1); + p4 = p4 + (rx[timing_est + k + ofdm_samplesperframe] * csam2); + } + + /* Calculate sample rate of phase samples, we are sampling phase + of pilot at half a symbol intervals */ + + float Fs1 = ofdm_fs / ((ofdm_m + ofdm_ncp) / 2); + + /* + * subtract phase of adjacent samples, rate of change of phase is + * frequency est. We combine samples from either end of frame to + * improve estimate. Small real 1E-12 term to prevent instability + * with 0 inputs. + */ + + ofdm->foff_metric = 0.9f * ofdm->foff_metric + 0.1f * (std::conj(p1) * p2 + std::conj(p3) * p4); + foff_est = Fs1 * std::arg(ofdm->foff_metric + 1E-12f) / TAU; + + if (ofdm->verbose > 1) { + fprintf(stderr, " foff_metric: %f %f foff_est: %f\n", std::real(ofdm->foff_metric), std::imag(ofdm->foff_metric), (double) foff_est); + } + + return foff_est; +} + +/* + * ---------------------------------------------- + * ofdm_txframe - modulates one frame of symbols + * ---------------------------------------------- + */ + +void ofdm_txframe(struct OFDM *ofdm, std::complex *tx, std::complex *tx_sym_lin) { + int sx = ofdm_ns; + int sy = ofdm_nc + 2; + std::complex *aframe = new std::complex[sx*sy]; + std::complex *asymbol = new std::complex[ofdm_m]; + std::complex *asymbol_cp = new std::complex[ofdm_m + ofdm_ncp]; + int i, j, k, m; + + /* initialize aframe to complex zero */ + + for (i = 0; i < ofdm_ns; i++) { + for (j = 0; j < (ofdm_nc + 2); j++) { + aframe[i*sy+j] = 0.0f; + } + } + + /* copy in a row of complex pilots to first row */ + + for (i = 0; i < (ofdm_nc + 2); i++) { + aframe[0*sy+i] = ofdm->pilots[i]; + } + + /* + * Place symbols in multi-carrier frame with pilots + * This will place boundary values of complex zero around data + */ + + for (i = 1; i <= ofdm_rowsperframe; i++) { + + /* copy in the Nc complex values with [0 Nc 0] or (Nc + 2) total */ + + for (j = 1; j < (ofdm_nc + 1); j++) { + aframe[i*sy+j] = tx_sym_lin[((i - 1) * ofdm_nc) + (j - 1)]; + } + } + + /* OFDM up-convert symbol by symbol so we can add CP */ + + for (i = 0, m = 0; i < ofdm_ns; i++, m += (ofdm_m + ofdm_ncp)) { + idft(ofdm, asymbol, &aframe[i*sy]); + + /* Copy the last Ncp samples to the front */ + + for (j = (ofdm_m - ofdm_ncp), k = 0; j < ofdm_m; j++, k++) { + asymbol_cp[k] = asymbol[j]; + } + + /* Now copy the all samples for this row after it */ + + for (j = ofdm_ncp, k = 0; k < ofdm_m; j++, k++) { + asymbol_cp[j] = asymbol[k]; + } + + /* Now move row to the tx output */ + + for (j = 0; j < (ofdm_m + ofdm_ncp); j++) { + tx[m + j] = asymbol_cp[j]; + } + } + + /* optional Tx Band Pass Filter */ + + if (ofdm->tx_bpf_en == true) { + std::complex *tx_filt = new std::complex[ofdm_samplesperframe]; + + quisk_ccfFilter(tx, tx_filt, ofdm_samplesperframe, ofdm->ofdm_tx_bpf); + memcpy(tx, tx_filt, ofdm_samplesperframe * sizeof (std::complex)); + delete[] tx_filt; + } + + delete[] asymbol_cp; + delete[] asymbol; + delete[] aframe; + } + +struct OFDM_CONFIG *ofdm_get_config_param() { + return &ofdm_config; +} + +int ofdm_get_nin(struct OFDM *ofdm) { + return ofdm->nin; +} + +int ofdm_get_samples_per_frame() { + return ofdm_samplesperframe; +} + +int ofdm_get_max_samples_per_frame() { + return 2 * ofdm_max_samplesperframe; +} + +int ofdm_get_bits_per_frame() { + return ofdm_bitsperframe; +} + +void ofdm_set_verbose(struct OFDM *ofdm, int level) { + ofdm->verbose = level; +} + +void ofdm_set_timing_enable(struct OFDM *ofdm, bool val) { + ofdm->timing_en = val; + + if (ofdm->timing_en == false) { + /* manually set ideal timing instant */ + + ofdm->sample_point = (ofdm_ncp - 1); + } +} + +void ofdm_set_foff_est_enable(struct OFDM *ofdm, bool val) { + ofdm->foff_est_en = val; +} + +void ofdm_set_phase_est_enable(struct OFDM *ofdm, bool val) { + ofdm->phase_est_en = val; +} + +void ofdm_set_off_est_hz(struct OFDM *ofdm, float val) { + ofdm->foff_est_hz = val; +} + +void ofdm_set_tx_bpf(struct OFDM *ofdm, bool val) { + if (val == true) { + allocate_tx_bpf(ofdm); + ofdm->tx_bpf_en = true; + } else { + if (ofdm->ofdm_tx_bpf) + deallocate_tx_bpf(ofdm); + ofdm->tx_bpf_en = false; + } +} + +/* + * -------------------------------------- + * ofdm_mod - modulates one frame of bits + * -------------------------------------- + */ + +void ofdm_mod(struct OFDM *ofdm, COMP *result, const int *tx_bits) { + int length = ofdm_bitsperframe / ofdm_bps; + std::complex *tx = (std::complex *) &result[0]; // complex has same memory layout + std::complex *tx_sym_lin = new std::complex[length]; + int dibit[2]; + int s, i; + + if (ofdm_bps == 1) { + /* Here we will have Nbitsperframe / 1 */ + + for (s = 0; s < length; s++) { + tx_sym_lin[s] = (float) (2 * tx_bits[s] - 1); + } + } else if (ofdm_bps == 2) { + /* Here we will have Nbitsperframe / 2 */ + + for (s = 0, i = 0; i < length; s += 2, i++) { + dibit[0] = tx_bits[s + 1] & 0x1; + dibit[1] = tx_bits[s ] & 0x1; + tx_sym_lin[i] = qpsk_mod(dibit); + } + } + + ofdm_txframe(ofdm, tx, tx_sym_lin); + delete[] tx_sym_lin; + } + +/* + * ---------------------------------------------------------------------------------- + * ofdm_sync_search - attempts to find coarse sync parameters for modem initial sync + * ---------------------------------------------------------------------------------- + */ + +/* + * This is a wrapper to maintain the older functionality + * with an array of COMPs as input + */ +int ofdm_sync_search(struct OFDM *ofdm, COMP *rxbuf_in) { + std::complex *rx = (std::complex *) &rxbuf_in[0]; // complex has same memory layout + int i, j; + + /* + * insert latest input samples into rxbuf + * so it is primed for when we have to + * call ofdm_demod() + */ + for (i = 0, j = ofdm->nin; i < (ofdm_rxbuf - ofdm->nin); i++, j++) { + ofdm->rxbuf[i] = ofdm->rxbuf[j]; + } + + /* insert latest input samples onto tail of rxbuf */ + + for (i = (ofdm_rxbuf - ofdm->nin), j = 0; i < ofdm_rxbuf; i++, j++) { + ofdm->rxbuf[i] = rx[j]; + } + + return(ofdm_sync_search_core(ofdm)); +} + +/* + * This is a wrapper with a new interface to reduce memory allocated. + * This works with ofdm_demod and freedv_api. + */ +int ofdm_sync_search_shorts(struct OFDM *ofdm, short *rxbuf_in, float gain) { + int i, j; + + /* shift the buffer left based on nin */ + + for (i = 0, j = ofdm->nin; i < (ofdm_rxbuf - ofdm->nin); i++, j++) { + ofdm->rxbuf[i] = ofdm->rxbuf[j]; + } + + /* insert latest input samples onto tail of rxbuf */ + + for (i = (ofdm_rxbuf - ofdm->nin), j = 0; i < ofdm_rxbuf; i++, j++) { + ofdm->rxbuf[i] = ((float)rxbuf_in[j] * gain); + } + + return(ofdm_sync_search_core(ofdm)); +} + +/* + * This is the rest of the function which expects that the data is + * already in ofdm->rxbuf + */ +static int ofdm_sync_search_core(struct OFDM *ofdm) { + /* Attempt coarse timing estimate (i.e. detect start of frame) */ + + int st = ofdm_m + ofdm_ncp + ofdm_samplesperframe; + int en = st + 2 * ofdm_samplesperframe; + int ct_est = est_timing(ofdm, &ofdm->rxbuf[st], (en - st)); + + ofdm->coarse_foff_est_hz = est_freq_offset(ofdm, &ofdm->rxbuf[st], ct_est); + + if (ofdm->verbose) { + fprintf(stderr, " ct_est: %4d foff_est: %4.1f timing_valid: %d timing_mx: %5.4f\n", + ct_est, (double) ofdm->coarse_foff_est_hz, ofdm->timing_valid, (double) ofdm->timing_mx); + } + + if (ofdm->timing_valid) { + /* potential candidate found .... */ + + /* calculate number of samples we need on next buffer to get into sync */ + + ofdm->nin = ofdm_samplesperframe + ct_est; + + /* reset modem states */ + + ofdm->sample_point = ofdm->timing_est = 0; + ofdm->foff_est_hz = ofdm->coarse_foff_est_hz; + } else { + ofdm->nin = ofdm_samplesperframe; + } + + return ofdm->timing_valid; +} + +/* + * ------------------------------------------ + * ofdm_demod - Demodulates one frame of bits + * ------------------------------------------ + */ + +/* This is a wrapper to maintain the older functionality with an + * array of COMPs as input */ +void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { + std::complex *rx = (std::complex *) &rxbuf_in[0]; // complex has same memory layout + int i, j; + + /* shift the buffer left based on nin */ + for (i = 0, j = ofdm->nin; i < (ofdm_rxbuf - ofdm->nin); i++, j++) { + ofdm->rxbuf[i] = ofdm->rxbuf[j]; + } + + /* insert latest input samples onto tail of rxbuf */ + for (i = (ofdm_rxbuf - ofdm->nin), j = 0; i < ofdm_rxbuf; i++, j++) { + ofdm->rxbuf[i] = rx[j]; + } + + ofdm_demod_core(ofdm, rx_bits); +} + +/* This is a wrapper with a new interface to reduce memory allocated. + * This works with ofdm_demod and freedv_api. */ +void ofdm_demod_shorts(struct OFDM *ofdm, int *rx_bits, short *rxbuf_in, float gain) { + int i, j; + + /* shift the buffer left based on nin */ + for (i = 0, j = ofdm->nin; i < (ofdm_rxbuf - ofdm->nin); i++, j++) { + ofdm->rxbuf[i] = ofdm->rxbuf[j]; + } + + /* insert latest input samples onto tail of rxbuf */ + for (i = (ofdm_rxbuf - ofdm->nin), j = 0; i < ofdm_rxbuf; i++, j++) { + ofdm->rxbuf[i] = ((float)rxbuf_in[j] * gain); + } + + ofdm_demod_core(ofdm, rx_bits); +} + +/* + * This is the rest of the function which expects that the data is + * already in ofdm->rxbuf + */ +static void ofdm_demod_core(struct OFDM *ofdm, int *rx_bits) { + std::complex aphase_est_pilot_rect; + float *aphase_est_pilot = new float[ofdm_nc + 2]; + float *aamp_est_pilot = new float[ofdm_nc + 2]; + float freq_err_hz; + int i, j, k, rr, st, en, ft_est; + int prev_timing_est = ofdm->timing_est; + + /* + * get user and calculated freq offset + */ + + float woff_est = TAU * ofdm->foff_est_hz / ofdm_fs; + + /* update timing estimate -------------------------------------------------- */ + + if (ofdm->timing_en == true) { + /* update timing at start of every frame */ + + st = ((ofdm_m + ofdm_ncp) + ofdm_samplesperframe) - floorf(ofdm_ftwindowwidth / 2) + ofdm->timing_est; + en = st + ofdm_samplesperframe - 1 + (ofdm_m + ofdm_ncp) + ofdm_ftwindowwidth; + + std::complex *work = new std::complex[(en - st)]; + + /* + * Adjust for the frequency error by shifting the phase + * using a conjugate multiply + */ + + for (i = st, j = 0; i < en; i++, j++) { + float tval = woff_est * i; + work[j] = ofdm->rxbuf[i] * cmplxconj(tval); + } + + ft_est = est_timing(ofdm, work, (en - st)); + ofdm->timing_est += (ft_est - ceilf(ofdm_ftwindowwidth / 2)); + + /* + * keep the freq est statistic updated in case we lose sync, + * note we supply it with uncorrected rxbuf, note + * ofdm->coarse_fest_off_hz is unused in normal operation, + * but stored for use in tofdm.c + */ + + ofdm->coarse_foff_est_hz = est_freq_offset(ofdm, &ofdm->rxbuf[st], ft_est); + + /* first frame in trial sync will have a better freq offset est - lets use it */ + + if (ofdm->frame_count == 0) { + ofdm->foff_est_hz = ofdm->coarse_foff_est_hz; + woff_est = TAU * ofdm->foff_est_hz / ofdm_fs; + } + + if (ofdm->verbose > 1) { + fprintf(stderr, " ft_est: %2d timing_est: %2d sample_point: %2d\n", ft_est, ofdm->timing_est, ofdm->sample_point); + } + + /* Black magic to keep sample_point inside cyclic prefix. Or something like that. */ + + ofdm->sample_point = max(ofdm->timing_est + (ofdm_ncp / 4), ofdm->sample_point); + ofdm->sample_point = min(ofdm->timing_est + ofdm_ncp, ofdm->sample_point); + delete[] work; + } + + /* + * Convert the time-domain samples to the frequency-domain using the rx_sym + * data matrix. This will be Nc+2 carriers of 11 symbols. + * + * You will notice there are Nc+2 BPSK symbols for each pilot symbol, and that there + * are Nc QPSK symbols for each data symbol. + * + * XXXXXXXXXXXXXXXXX <-- Timing Slip + * PPPPPPPPPPPPPPPPPPP <-- Previous Frames Pilot + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD Ignore these past data symbols + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * PPPPPPPPPPPPPPPPPPP <-- This Frames Pilot + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD These are the current data symbols to be decoded + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * PPPPPPPPPPPPPPPPPPP <-- Next Frames Pilot + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD Ignore these next data symbols + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * DDDDDDDDDDDDDDDDD + * PPPPPPPPPPPPPPPPPPP <-- Future Frames Pilot + * XXXXXXXXXXXXXXXXX <-- Timing Slip + * + * So this algorithm will have seven data symbols and four pilot symbols to process. + * The average of the four pilot symbols is our phase estimation. + */ + + for (i = 0; i < (ofdm_ns + 3); i++) { + for (j = 0; j < (ofdm_nc + 2); j++) { + ofdm->rx_sym[i][j] = 0.0f; + } + } + + /* + * "Previous" pilot symbol is one modem frame above. + */ + + st = (ofdm_m + ofdm_ncp) + 1 + ofdm->sample_point; + en = st + ofdm_m; + + std::complex *work = new std::complex[ofdm_m]; + + /* down-convert at current timing instant---------------------------------- */ + + for (j = st, k = 0; j < en; j++, k++) { + float tval = woff_est * j; + work[k] = ofdm->rxbuf[j] * cmplxconj(tval); + } + + /* + * Each symbol is of course (ofdm_m + ofdm_ncp) samples long and + * becomes Nc+2 carriers after DFT. + * + * We put this carrier pilot symbol at the top of our matrix: + * + * 1 .................. Nc+2 + * + * +----------------------+ + * | Previous Pilot | rx_sym[0] + * +----------------------+ + * | | + * + */ + + dft(ofdm, ofdm->rx_sym[0], work); + + /* + * "This" pilot comes after the extra symbol allotted at the top, and after + * the "previous" pilot and previous data symbols (let's call it, the previous + * modem frame). + * + * So we will now be starting at "this" pilot symbol, and continuing to the + * "next" pilot symbol. + * + * In this routine we also process the current data symbols. + */ + + for (rr = 0; rr < (ofdm_ns + 1); rr++) { + st = (ofdm_m + ofdm_ncp) + ofdm_samplesperframe + (rr * (ofdm_m + ofdm_ncp)) + 1 + ofdm->sample_point; + en = st + ofdm_m; + + /* down-convert at current timing instant---------------------------------- */ + + for (j = st, k = 0; j < en; j++, k++) { + float tval = woff_est * j; + work[k] = ofdm->rxbuf[j] * cmplxconj(tval); + } + + /* + * We put these Nc+2 carrier symbols into our matrix after the previous pilot: + * + * 1 .................. Nc+2 + * + * | Previous Pilot | rx_sym[0] + * +----------------------+ + * | This Pilot | rx_sym[1] + * +----------------------+ + * | Data | rx_sym[2] + * +----------------------+ + * | Data | rx_sym[3] + * +----------------------+ + * | Data | rx_sym[4] + * +----------------------+ + * | Data | rx_sym[5] + * +----------------------+ + * | Data | rx_sym[6] + * +----------------------+ + * | Data | rx_sym[7] + * +----------------------+ + * | Data | rx_sym[8] + * +----------------------+ + * | Next Pilot | rx_sym[9] + * +----------------------+ + * | | rx_sym[10] + */ + + dft(ofdm, ofdm->rx_sym[rr + 1], work); + } + + /* + * OK, now we want to process to the "future" pilot symbol. This is after + * the "next" modem frame. + * + * We are ignoring the data symbols between the "next" pilot and "future" pilot. + * We only want the "future" pilot symbol, to perform the averaging of all pilots. + */ + + st = (ofdm_m + ofdm_ncp) + (3 * ofdm_samplesperframe) + 1 + ofdm->sample_point; + en = st + ofdm_m; + + /* down-convert at current timing instant---------------------------------- */ + + for (j = st, k = 0; j < en; j++, k++) { + float tval = woff_est * j; + work[k] = ofdm->rxbuf[j] * cmplxconj(tval); + } + + /* + * We put the future pilot after all the previous symbols in the matrix: + * + * 1 .................. Nc+2 + * + * | | rx_sym[9] + * +----------------------+ + * | Future Pilot | rx_sym[10] + * +----------------------+ + */ + + dft(ofdm, ofdm->rx_sym[ofdm_ns + 2], work); + + /* + * We are finished now with the DFT and down conversion + * From here on down we are in frequency domain + */ + + /* est freq err based on all carriers ------------------------------------ */ + + if (ofdm->foff_est_en == true) { + /* + * sym[1] is (this) pilot symbol, sym[9] is (next) pilot symbol. + * + * By subtracting the two averages of these pilots, we find the frequency + * by the change in phase over time. + */ + + std::complex freq_err_rect = + std::conj(vector_sum(ofdm->rx_sym[1], ofdm_nc + 2)) * + vector_sum(ofdm->rx_sym[ofdm_ns + 1], ofdm_nc + 2); + + /* prevent instability in atan(im/re) when real part near 0 */ + + freq_err_rect += 1E-6f; + + freq_err_hz = std::arg(freq_err_rect) * ofdm_rs / (TAU * ofdm_ns); + ofdm->foff_est_hz += (ofdm->foff_est_gain * freq_err_hz); + } + + /* OK - now estimate and correct pilot phase ---------------------------------- */ + + for (i = 0; i < (ofdm_nc + 2); i++) { + aphase_est_pilot[i] = 10.0f; + aamp_est_pilot[i] = 0.0f; + } + + /* + * Basically we divide the Nc+2 pilots into groups of 3 + * + * Then average the phase surrounding each of the data symbols. + */ + + for (i = 1; i < (ofdm_nc + 1); i++) { + std::complex symbol[3]; + + for (j = (i - 1), k = 0; j < (i + 2); j++, k++) { + symbol[k] = ofdm->rx_sym[1][j] * std::conj(ofdm->pilots[j]); /* this pilot conjugate */ + } + + aphase_est_pilot_rect = vector_sum(symbol, 3); + + for (j = (i - 1), k = 0; j < (i + 2); j++, k++) { + symbol[k] = ofdm->rx_sym[ofdm_ns + 1][j] * std::conj(ofdm->pilots[j]); /* next pilot conjugate */ + } + + aphase_est_pilot_rect = aphase_est_pilot_rect + vector_sum(symbol, 3); + + /* use next step of pilots in past and future */ + + for (j = (i - 1), k = 0; j < (i + 2); j++, k++) { + symbol[k] = ofdm->rx_sym[0][j] * std::conj(ofdm->pilots[j]); /* previous pilot */ + } + + aphase_est_pilot_rect = aphase_est_pilot_rect + vector_sum(symbol, 3); + + for (j = (i - 1), k = 0; j < (i + 2); j++, k++) { + symbol[k] = ofdm->rx_sym[ofdm_ns + 2][j] * std::conj(ofdm->pilots[j]); /* last pilot */ + } + + aphase_est_pilot_rect = aphase_est_pilot_rect + vector_sum(symbol, 3); + aphase_est_pilot[i] = std::arg(aphase_est_pilot_rect); + + /* amplitude is estimated over 12 pilots */ + + aamp_est_pilot[i] = std::abs(aphase_est_pilot_rect / 12.0f); + } + + /* + * correct phase offset using phase estimate, and demodulate + * bits, separate loop as it runs across cols (carriers) to get + * frame bit ordering correct + */ + + std::complex rx_corr; + int abit[2]; + int bit_index = 0; + float sum_amp = 0.0f; + + for (rr = 0; rr < ofdm_rowsperframe; rr++) { + /* + * Note the i starts with the second carrier, ends with Nc+1. + * so we ignore the first and last carriers. + * + * Also note we are using sym[2..8] or the seven data symbols. + */ + + for (i = 1; i < (ofdm_nc + 1); i++) { + if (ofdm->phase_est_en == true) { + rx_corr = ofdm->rx_sym[rr + 2][i] * cmplxconj(aphase_est_pilot[i]); + } else { + rx_corr = ofdm->rx_sym[rr + 2][i]; + } + + /* + * Output complex data symbols after phase correction; + * rx_np means the pilot symbols have been removed + */ + + ofdm->rx_np[(rr * ofdm_nc) + (i - 1)] = rx_corr; + + /* + * Note even though amp ests are the same for each col, + * the FEC decoder likes to have one amplitude per symbol + * so convenient to log them all + */ + + ofdm->rx_amp[(rr * ofdm_nc) + (i - 1)] = aamp_est_pilot[i]; + sum_amp += aamp_est_pilot[i]; + + /* + * Note like amps in this implementation phase ests are the + * same for each col, but we log them for each symbol anyway + */ + + ofdm->aphase_est_pilot_log[(rr * ofdm_nc) + (i - 1)] = aphase_est_pilot[i]; + + if (ofdm_bps == 1) { + rx_bits[bit_index++] = std::real(rx_corr) > 0.0f; + } else if (ofdm_bps == 2) { + /* + * Only one final task, decode what quadrant the phase + * is in, and return the dibits + */ + qpsk_demod(rx_corr, abit); + rx_bits[bit_index++] = abit[1]; + rx_bits[bit_index++] = abit[0]; + } + } + } + + /* update mean amplitude estimate for LDPC decoder scaling */ + + ofdm->mean_amp = 0.9f * ofdm->mean_amp + 0.1f * sum_amp / (ofdm_rowsperframe * ofdm_nc); + + /* Adjust nin to take care of sample clock offset */ + + ofdm->nin = ofdm_samplesperframe; + + if (ofdm->timing_en == true) { + ofdm->clock_offset_counter += prev_timing_est - ofdm->timing_est; + + int thresh = (ofdm_m + ofdm_ncp) / 8; + int tshift = (ofdm_m + ofdm_ncp) / 4; + + if (ofdm->timing_est > thresh) { + ofdm->nin = ofdm_samplesperframe + tshift; + ofdm->timing_est -= tshift; + ofdm->sample_point -= tshift; + } else if (ofdm->timing_est < -thresh) { + ofdm->nin = ofdm_samplesperframe - tshift; + ofdm->timing_est += tshift; + ofdm->sample_point += tshift; + } + } + + /* + * estimate signal and noise power, see ofdm_lib.m, + * cohpsk.m for more info + */ + + std::complex *rx_np = ofdm->rx_np; + + float sig_var = 0.0f; + + for (i = 0; i < (ofdm_rowsperframe * ofdm_nc); i++) { + sig_var += cnormf(rx_np[i]); + } + + sig_var /= (ofdm_rowsperframe * ofdm_nc); + float sig_rms = sqrtf(sig_var); + + float sum_x = 0.0f; + float sum_xx = 0.0f; + int n = 0; + + for (i = 0; i < (ofdm_rowsperframe * ofdm_nc); i++) { + std::complex s = rx_np[i]; + + if (fabsf(std::real(s)) > sig_rms) { + sum_x += std::imag(s); + sum_xx += std::imag(s) * std::imag(s); + n++; + } + } + + /* + * with large interfering carriers this alg can break down - in + * that case set a benign value for noise_var that will produce a + * sensible (probably low) SNR est + */ + + float noise_var = 1.0f; + + if (n > 1) { + noise_var = (n * sum_xx - sum_x * sum_x) / (n * (n - 1)); + } + + ofdm->noise_var = 2.0f * noise_var; + ofdm->sig_var = sig_var; + + delete[] work; + delete[] aamp_est_pilot; + delete[] aphase_est_pilot; +} + +/* iterate state machine ------------------------------------*/ + +void ofdm_sync_state_machine(struct OFDM *ofdm, uint8_t *rx_uw) { + int i; + + State next_state = ofdm->sync_state; + + ofdm->sync_start = false; + ofdm->sync_end = false; + + if (ofdm->sync_state == search) { + if (ofdm->timing_valid) { + ofdm->frame_count = 0; + ofdm->sync_counter = 0; + ofdm->sync_start = true; + ofdm->clock_offset_counter = 0; + next_state = trial; + } + } + + if ((ofdm->sync_state == synced) || (ofdm->sync_state == trial)) { + ofdm->frame_count++; + ofdm->frame_count_interleaver++; + + /* + * freq offset est may be too far out, and has aliases every 1/Ts, so + * we use a Unique Word to get a really solid indication of sync. + */ + + ofdm->uw_errors = 0; + + for (i = 0; i < ofdm_nuwbits; i++) { + ofdm->uw_errors += ofdm->tx_uw[i] ^ rx_uw[i]; + } + + /* + * during trial sync we don't tolerate errors so much, we look + * for 3 consecutive frames with low error rate to confirm sync + */ + + if (ofdm->sync_state == trial) { + if (ofdm->uw_errors > 2) { + /* if we exceed thresh stay in trial sync */ + + ofdm->sync_counter++; + ofdm->frame_count = 0; + } + + if (ofdm->sync_counter == 2) { + /* if we get two bad frames drop sync and start again */ + + next_state = search; + ofdm->sync_state_interleaver = search; + } + + if (ofdm->frame_count == 4) { + /* three good frames, sync is OK! */ + + next_state = synced; + } + } + + /* once we have synced up we tolerate a higher error rate to wait out fades */ + + if (ofdm->sync_state == synced) { + if (ofdm->uw_errors > 2) { + ofdm->sync_counter++; + } else { + ofdm->sync_counter = 0; + } + + if ((ofdm->sync_mode == autosync) && (ofdm->sync_counter == 12)) { + /* run of consecutive bad frames ... drop sync */ + + next_state = search; + ofdm->sync_state_interleaver = search; + } + } + } + + ofdm->last_sync_state = ofdm->sync_state; + ofdm->last_sync_state_interleaver = ofdm->sync_state_interleaver; + ofdm->sync_state = next_state; +} + +/*---------------------------------------------------------------------------* \ + + FUNCTIONS...: ofdm_set_sync + AUTHOR......: David Rowe + DATE CREATED: May 2018 + + Operator control of sync state machine. This mode is required to + acquire sync up at very low SNRS. This is difficult to implement, + for example we may get a false sync, or the state machine may fall + out of sync by mistake during a long fade. + + So with this API call we allow some operator assistance. + + Ensure this is called in the same thread as ofdm_sync_state_machine(). + +\*---------------------------------------------------------------------------*/ + +void ofdm_set_sync(struct OFDM *ofdm, Sync sync_cmd) { + assert(ofdm != NULL); + + switch (sync_cmd) { + case unsync: + /* + * force manual unsync, in case operator detects false sync, + * which will cause sync state machine to have another go at sync + */ + ofdm->sync_state = search; + ofdm->sync_state_interleaver = search; + break; + case autosync: + /* normal operating mode - sync state machine decides when to unsync */ + + ofdm->sync_mode = autosync; + break; + case manualsync: + /* + * allow sync state machine to sync, but not to unsync, the + * operator will decide that manually + */ + ofdm->sync_mode = manualsync; + break; + default: + assert(0); + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: ofdm_get_demod_stats() + AUTHOR......: David Rowe + DATE CREATED: May 2018 + + Fills stats structure with a bunch of demod information. + +\*---------------------------------------------------------------------------*/ + +void ofdm_get_demod_stats(struct OFDM *ofdm, struct MODEM_STATS *stats) { + int c, r; + + stats->Nc = ofdm_nc; + assert(stats->Nc <= MODEM_STATS_NC_MAX); + + float snr_est = 10.0f * log10f((0.1f + (ofdm->sig_var / ofdm->noise_var)) * ofdm_nc * ofdm_rs / 3000.0f); + float total = ofdm->frame_count * ofdm_samplesperframe; + + stats->snr_est = 0.9f * stats->snr_est + 0.1f * snr_est; + stats->sync = ((ofdm->sync_state == synced) || (ofdm->sync_state == trial)); + stats->foff = ofdm->foff_est_hz; + stats->rx_timing = ofdm->timing_est; + stats->clock_offset = 0; + + if (total != 0.0f) { + stats->clock_offset = ofdm->clock_offset_counter / total; + } + + stats->sync_metric = ofdm->timing_mx; + + assert(ofdm_rowsperframe < MODEM_STATS_NR_MAX); + stats->nr = ofdm_rowsperframe; + + for (c = 0; c < ofdm_nc; c++) { + for (r = 0; r < ofdm_rowsperframe; r++) { + std::complex rot = ofdm->rx_np[r * c] * cmplx(ROT45); + + stats->rx_symbols[r][c].real = std::real(rot); + stats->rx_symbols[r][c].imag = std::imag(rot); + } + } +} + +/* Assemble modem frame of bits from UW, payload bits, and txt bits */ + +void ofdm_assemble_modem_frame(struct OFDM *ofdm, uint8_t modem_frame[], uint8_t payload_bits[], uint8_t txt_bits[]) { + int b, t; + + int p = 0; + int u = 0; + + for (b = 0; b < (ofdm_bitsperframe - ofdm_ntxtbits); b++) { + if ((u < ofdm_nuwbits) && (b == uw_ind[u])) { + modem_frame[b] = ofdm->tx_uw[u++]; + } else { + modem_frame[b] = payload_bits[p++]; + } + } + + assert(u == ofdm_nuwbits); + assert(p == (ofdm_bitsperframe - ofdm_nuwbits - ofdm_ntxtbits)); + + for (t = 0; b < ofdm_bitsperframe; b++, t++) { + modem_frame[b] = txt_bits[t]; + } + + assert(t == ofdm_ntxtbits); +} + +/* Assemble modem frame from UW, payload symbols, and txt bits */ + +void ofdm_assemble_modem_frame_symbols(std::complex modem_frame[], COMP payload_syms[], uint8_t txt_bits[]) { + std::complex *payload = (std::complex *) &payload_syms[0]; // complex has same memory layout + int Nsymsperframe = ofdm_bitsperframe / ofdm_bps; + int Nuwsyms = ofdm_nuwbits / ofdm_bps; + int Ntxtsyms = ofdm_ntxtbits / ofdm_bps; + int dibit[2]; + int s, t; + + int p = 0; + int u = 0; + + for (s = 0; s < Nsymsperframe - Ntxtsyms; s++) { + if ((u < Nuwsyms) && (s == uw_ind_sym[u])) { + modem_frame[s] = tx_uw_syms[u++]; + } else { + modem_frame[s] = payload[p++]; + } + } + + assert(u == Nuwsyms); + assert(p == (Nsymsperframe - Nuwsyms - Ntxtsyms)); + + for (t = 0; s < Nsymsperframe; s++, t += ofdm_bps) { + dibit[0] = txt_bits[t + 1] & 0x1; + dibit[1] = txt_bits[t] & 0x1; + modem_frame[s] = qpsk_mod(dibit); + } + + assert(t == ofdm_ntxtbits); +} + +void ofdm_disassemble_modem_frame(struct OFDM *ofdm, uint8_t rx_uw[], + COMP codeword_syms[], + float codeword_amps[], + short txt_bits[]) { + std::complex *codeword = (std::complex *) &codeword_syms[0]; // complex has same memory layout + int Nsymsperframe = ofdm_bitsperframe / ofdm_bps; + int Nuwsyms = ofdm_nuwbits / ofdm_bps; + int Ntxtsyms = ofdm_ntxtbits / ofdm_bps; + int dibit[2]; + int s, t; + + int p = 0; + int u = 0; + + for (s = 0; s < (Nsymsperframe - Ntxtsyms); s++) { + if ((u < Nuwsyms) && (s == uw_ind_sym[u])) { + qpsk_demod(ofdm->rx_np[s], dibit); + + rx_uw[ofdm_bps * u ] = dibit[1]; + rx_uw[ofdm_bps * u + 1] = dibit[0]; + u++; + } else { + codeword[p] = ofdm->rx_np[s]; + codeword_amps[p] = ofdm->rx_amp[s]; + p++; + } + } + + assert(u == Nuwsyms); + assert(p == (Nsymsperframe - Nuwsyms - Ntxtsyms)); + + for (t = 0; s < Nsymsperframe; s++, t += ofdm_bps) { + qpsk_demod(ofdm->rx_np[s], dibit); + + txt_bits[t ] = dibit[1]; + txt_bits[t + 1] = dibit[0]; + } + + assert(t == ofdm_ntxtbits); +} + + +/* + * Pseudo-random number generator that we can implement in C with + * identical results to Octave. Returns an unsigned int between 0 + * and 32767. Used for generating test frames of various lengths. + */ + +void ofdm_rand(uint16_t r[], int n) { + uint64_t seed = 1; + int i; + + for (i = 0; i < n; i++) { + seed = (1103515245l * seed + 12345) % 32768; + r[i] = seed; + } +} + + +void ofdm_generate_payload_data_bits(uint8_t payload_data_bits[], int data_bits_per_frame) { + uint16_t *r = new uint16_t[data_bits_per_frame]; + int i; + + /* construct payload data bits */ + + ofdm_rand(r, data_bits_per_frame); + + for(i=0; i 16384; + } + + delete[] r; +} + +void ofdm_print_info(struct OFDM *ofdm) { + const char *syncmode[] = { + "unsync", + "autosync", + "manualsync" + }; + + fprintf(stderr, "ofdm->foff_est_gain = %g\n", (double)ofdm->foff_est_gain); + fprintf(stderr, "ofdm->foff_est_hz = %g\n", (double)ofdm->foff_est_hz); + fprintf(stderr, "ofdm->timing_mx = %g\n", (double)ofdm->timing_mx); + fprintf(stderr, "ofdm->coarse_foff_est_hz = %g\n", (double)ofdm->coarse_foff_est_hz); + fprintf(stderr, "ofdm->timing_norm = %g\n", (double)ofdm->timing_norm); + fprintf(stderr, "ofdm->sig_var = %g\n", (double)ofdm->sig_var); + fprintf(stderr, "ofdm->noise_var = %g\n", (double)ofdm->noise_var); + fprintf(stderr, "ofdm->mean_amp = %g\n", (double)ofdm->mean_amp); + fprintf(stderr, "ofdm->clock_offset_counter = %d\n", ofdm->clock_offset_counter); + fprintf(stderr, "ofdm->verbose = %d\n", ofdm->verbose); + fprintf(stderr, "ofdm->sample_point = %d\n", ofdm->sample_point); + fprintf(stderr, "ofdm->timing_est = %d\n", ofdm->timing_est); + fprintf(stderr, "ofdm->timing_valid = %d\n", ofdm->timing_valid); + fprintf(stderr, "ofdm->nin = %d\n", ofdm->nin); + fprintf(stderr, "ofdm->uw_errors = %d\n", ofdm->uw_errors); + fprintf(stderr, "ofdm->sync_counter = %d\n", ofdm->sync_counter); + fprintf(stderr, "ofdm->frame_count = %d\n", ofdm->frame_count); + fprintf(stderr, "ofdm->sync_start = %s\n", ofdm->sync_start ? "true" : "false"); + fprintf(stderr, "ofdm->sync_end = %s\n", ofdm->sync_end ? "true" : "false"); + fprintf(stderr, "ofdm->sync_mode = %s\n", syncmode[ofdm->sync_mode]); + fprintf(stderr, "ofdm->frame_count_interleaver = %d\n", ofdm->frame_count_interleaver); + fprintf(stderr, "ofdm->timing_en = %s\n", ofdm->timing_en ? "true" : "false"); + fprintf(stderr, "ofdm->foff_est_en = %s\n", ofdm->foff_est_en ? "true" : "false"); + fprintf(stderr, "ofdm->phase_est_en = %s\n", ofdm->phase_est_en ? "true" : "false"); + fprintf(stderr, "ofdm->tx_bpf_en = %s\n", ofdm->tx_bpf_en ? "true" : "false"); +}; + +} // FreeDV diff --git a/libfreedv/ofdm_internal.h b/libfreedv/ofdm_internal.h index 6dcd3f823..6c26d3796 100644 --- a/libfreedv/ofdm_internal.h +++ b/libfreedv/ofdm_internal.h @@ -43,8 +43,8 @@ #define TAU (2.0f * M_PI) #define ROT45 (M_PI / 4.0f) -#define cmplx(value) (COSF(value) + SINF(value) * I) -#define cmplxconj(value) (COSF(value) + SINF(value) * -I) +#define cmplx(value) (std::complex{cos(value), sin(value)}) +#define cmplxconj(value) (std::complex{cos(value), -sin(value)}) namespace FreeDV { diff --git a/libfreedv/phi0.cpp b/libfreedv/phi0.cpp new file mode 100644 index 000000000..b5e18da4a --- /dev/null +++ b/libfreedv/phi0.cpp @@ -0,0 +1,223 @@ + +// phi0.c +// +// An approximation of the function +// +// This file is generated by the gen_phi0 scritps +// Any changes should be made to that file, not this one + +#include + +#define SI16(f) ((int32_t)(f * (1<<16))) + +namespace FreeDV +{ + +float phi0( float xf ) { + + int32_t x = SI16(xf); + + if (x >= SI16(10.0f)) return(0.0f); + else { + if (x >= SI16(5.0f)) { + int i = 19 - (x >> 15); + switch (i) { + case 0: return(0.000116589f); // (9.5) + case 1: return(0.000192223f); // (9.0) + case 2: return(0.000316923f); // (8.5) + case 3: return(0.000522517f); // (8.0) + case 4: return(0.000861485f); // (7.5) + case 5: return(0.001420349f); // (7.0) + case 6: return(0.002341760f); // (6.5) + case 7: return(0.003860913f); // (6.0) + case 8: return(0.006365583f); // (5.5) + case 9: return(0.010495133f); // (5.0) + } + } + else { + if (x >= SI16(1.0f)) { + int i = 79 - (x >> 12); + switch (i) { + case 0: return(0.013903889f); // (4.9375) + case 1: return(0.014800644f); // (4.8750) + case 2: return(0.015755242f); // (4.8125) + case 3: return(0.016771414f); // (4.7500) + case 4: return(0.017853133f); // (4.6875) + case 5: return(0.019004629f); // (4.6250) + case 6: return(0.020230403f); // (4.5625) + case 7: return(0.021535250f); // (4.5000) + case 8: return(0.022924272f); // (4.4375) + case 9: return(0.024402903f); // (4.3750) + case 10: return(0.025976926f); // (4.3125) + case 11: return(0.027652501f); // (4.2500) + case 12: return(0.029436184f); // (4.1875) + case 13: return(0.031334956f); // (4.1250) + case 14: return(0.033356250f); // (4.0625) + case 15: return(0.035507982f); // (4.0000) + case 16: return(0.037798579f); // (3.9375) + case 17: return(0.040237016f); // (3.8750) + case 18: return(0.042832850f); // (3.8125) + case 19: return(0.045596260f); // (3.7500) + case 20: return(0.048538086f); // (3.6875) + case 21: return(0.051669874f); // (3.6250) + case 22: return(0.055003924f); // (3.5625) + case 23: return(0.058553339f); // (3.5000) + case 24: return(0.062332076f); // (3.4375) + case 25: return(0.066355011f); // (3.3750) + case 26: return(0.070637993f); // (3.3125) + case 27: return(0.075197917f); // (3.2500) + case 28: return(0.080052790f); // (3.1875) + case 29: return(0.085221814f); // (3.1250) + case 30: return(0.090725463f); // (3.0625) + case 31: return(0.096585578f); // (3.0000) + case 32: return(0.102825462f); // (2.9375) + case 33: return(0.109469985f); // (2.8750) + case 34: return(0.116545700f); // (2.8125) + case 35: return(0.124080967f); // (2.7500) + case 36: return(0.132106091f); // (2.6875) + case 37: return(0.140653466f); // (2.6250) + case 38: return(0.149757747f); // (2.5625) + case 39: return(0.159456024f); // (2.5000) + case 40: return(0.169788027f); // (2.4375) + case 41: return(0.180796343f); // (2.3750) + case 42: return(0.192526667f); // (2.3125) + case 43: return(0.205028078f); // (2.2500) + case 44: return(0.218353351f); // (2.1875) + case 45: return(0.232559308f); // (2.1250) + case 46: return(0.247707218f); // (2.0625) + case 47: return(0.263863255f); // (2.0000) + case 48: return(0.281099022f); // (1.9375) + case 49: return(0.299492155f); // (1.8750) + case 50: return(0.319127030f); // (1.8125) + case 51: return(0.340095582f); // (1.7500) + case 52: return(0.362498271f); // (1.6875) + case 53: return(0.386445235f); // (1.6250) + case 54: return(0.412057648f); // (1.5625) + case 55: return(0.439469363f); // (1.5000) + case 56: return(0.468828902f); // (1.4375) + case 57: return(0.500301872f); // (1.3750) + case 58: return(0.534073947f); // (1.3125) + case 59: return(0.570354566f); // (1.2500) + case 60: return(0.609381573f); // (1.1875) + case 61: return(0.651427083f); // (1.1250) + case 62: return(0.696805010f); // (1.0625) + case 63: return(0.745880827f); // (1.0000) + } + } + else { + if (x > SI16(0.007812f)) { + if (x > SI16(0.088388f)) { + if (x > SI16(0.250000f)) { + if (x > SI16(0.500000f)) { + if (x > SI16(0.707107f)) { + return(0.922449644f); + } else { + return(1.241248638f); + } + } else { + if (x > SI16(0.353553f)) { + return(1.573515241f); + } else { + return(1.912825912f); + } + } + } else { + if (x > SI16(0.125000f)) { + if (x > SI16(0.176777f)) { + return(2.255740095f); + } else { + return(2.600476919f); + } + } else { + return(2.946130351f); + } + } + } else { + if (x > SI16(0.022097f)) { + if (x > SI16(0.044194f)) { + if (x > SI16(0.062500f)) { + return(3.292243417f); + } else { + return(3.638586634f); + } + } else { + if (x > SI16(0.031250f)) { + return(3.985045009f); + } else { + return(4.331560985f); + } + } + } else { + if (x > SI16(0.011049f)) { + if (x > SI16(0.015625f)) { + return(4.678105767f); + } else { + return(5.024664952f); + } + } else { + return(5.371231340f); + } + } + } + } else { + if (x > SI16(0.000691f)) { + if (x > SI16(0.001953f)) { + if (x > SI16(0.003906f)) { + if (x > SI16(0.005524f)) { + return(5.717801329f); + } else { + return(6.064373119f); + } + } else { + if (x > SI16(0.002762f)) { + return(6.410945809f); + } else { + return(6.757518949f); + } + } + } else { + if (x > SI16(0.000977f)) { + if (x > SI16(0.001381f)) { + return(7.104092314f); + } else { + return(7.450665792f); + } + } else { + return(7.797239326f); + } + } + } else { + if (x > SI16(0.000173f)) { + if (x > SI16(0.000345f)) { + if (x > SI16(0.000488f)) { + return(8.143812888f); + } else { + return(8.490386464f); + } + } else { + if (x > SI16(0.000244f)) { + return(8.836960047f); + } else { + return(9.183533634f); + } + } + } else { + if (x > SI16(0.000086f)) { + if (x > SI16(0.000122f)) { + return(9.530107222f); + } else { + return(9.876680812f); + } + } else { + return(10.000000000f); + } + } + } + } + } + } + } + return(10.0f); +} + +} // FreeDV diff --git a/libfreedv/phi0.h b/libfreedv/phi0.h new file mode 100644 index 000000000..e1e67d859 --- /dev/null +++ b/libfreedv/phi0.h @@ -0,0 +1,11 @@ +// phi0.h +#ifndef PHI0_H +#define PHI0_H + +namespace FreeDV { + +extern float phi0( float xf ); + +} // FreeDV + +#endif