VE7UWU/sdrangel
1
0
Fork 0
mirror of https://github.com/f4exb/sdrangel.git synced 2024-11-29 11:18:56 -05:00
Code Issues Projects Releases Wiki Activity

libfreedv: added 700D related stuff

Browse Source
This commit is contained in:
f4exb 2019-03-05 03:12:41 +01:00
parent 96836e7ff6
commit 0951b96b93
10 changed files with 3154 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
libfreedv/CMakeLists.txt
View File

@ -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

27
libfreedv/HRA_112_112.cpp Normal file
View File

@ -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 <stdint.h>
#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

23
libfreedv/HRA_112_112.h Normal file
View File

@ -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

122
libfreedv/gp_interleaver.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdio.h>
#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<sizeof(b_table)/(2*sizeof(int)); i+=2) {
if (b_table[i] == Nbits) {
return b_table[i+1];
}
}
/* if we get it means a Nbits we dont have in our table so choke */
assert(0);
return 0;
}
void gp_interleave_comp(COMP interleaved_frame[], COMP frame[], int Nbits) {
int b = choose_interleaver_b(Nbits);
int i,j;
for (i=0; i<Nbits; i++) {
j = (b*i) % Nbits;
interleaved_frame[j] = frame[i];
}
}
void gp_deinterleave_comp(COMP frame[], COMP interleaved_frame[], int Nbits) {
int b = choose_interleaver_b(Nbits);
int i,j;
for (i=0; i<Nbits; i++) {
j = (b*i) % Nbits;
frame[i] = interleaved_frame[j];
}
}
void gp_interleave_float(float interleaved_frame[], float frame[], int Nbits) {
int b = choose_interleaver_b(Nbits);
int i,j;
for (i=0; i<Nbits; i++) {
j = (b*i) % Nbits;
interleaved_frame[j] = frame[i];
}
}
void gp_deinterleave_float(float frame[], float interleaved_frame[], int Nbits) {
int b = choose_interleaver_b(Nbits);
int i,j;
for (i=0; i<Nbits; i++) {
j = (b*i) % Nbits;
frame[i] = interleaved_frame[j];
}
}
} // FreeDV

260
libfreedv/interldpc.cpp Normal file
View File

@ -0,0 +1,260 @@
/*---------------------------------------------------------------------------*\
FILE........: interldpc.c
AUTHOR......: David Rowe
DATE CREATED: April 2018
Helper functions for interleaved LDPC waveforms.
\*---------------------------------------------------------------------------*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#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<float> 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<float> s = std::complex<float>{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<float> 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<float> *tx_symbols = new std::complex<float>[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

726
libfreedv/mpdecode_core.cpp Normal file
View File

@ -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 <math.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <assert.h>
#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 <degree> 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; p<ldpc->NumberParityBits; p++) {
par = 0;
for (i=0; i<ldpc->max_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<NumberParityBits;i++) {
count = 0;
for (j=0;j<max_row_weight;j++) {
if ( H_rows[i+j*NumberParityBits] > 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<shift;k++){
count = 0;
for (j=0;j<max_row_weight;j++) {
if ( H_rows[cnt+j*NumberParityBits] > 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;i<NumberParityBits;i++) {
// Allocate sub nodes
c_nodes[i].subs = (struct c_sub_node*) calloc(c_nodes[i].degree, sizeof(struct c_sub_node));
assert(c_nodes[i].subs);
// Populate sub nodes
for (j=0;j<c_nodes[i].degree-2;j++) {
c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
}
j=c_nodes[i].degree-2;
if (i==0){
c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
}
else {
c_nodes[i].subs[j].index = (CodeLength-NumberParityBits)+i-1;
}
j=c_nodes[i].degree-1;
c_nodes[i].subs[j].index = (CodeLength-NumberParityBits)+i;
}
}
if (shift >0){
cnt=0;
for (i=0;i<(NumberParityBits/shift);i++){
for (k =0;k<shift;k++){
// Allocate sub nodes
c_nodes[cnt].subs = (struct c_sub_node*) calloc(c_nodes[cnt].degree, sizeof(struct c_sub_node));
assert(c_nodes[cnt].subs);
// Populate sub nodes
for (j=0;j<c_nodes[cnt].degree-2;j++) {
c_nodes[cnt].subs[j].index = (H_rows[cnt+j*NumberParityBits] - 1);
}
j=c_nodes[cnt].degree-2;
if ((i ==0)||(i==(NumberParityBits/shift-1))){
c_nodes[cnt].subs[j].index = (H_rows[cnt+j*NumberParityBits] - 1);
}
else{
c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i);
}
j=c_nodes[cnt].degree-1;
c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i+1);
if (i== (NumberParityBits/shift-1))
{
c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i);
}
cnt++;
}
}
}
} else {
for (i=0;i<NumberParityBits;i++) {
// Allocate sub nodes
c_nodes[i].subs = (struct c_sub_node*) calloc(c_nodes[i].degree, sizeof(struct c_sub_node));
assert(c_nodes[i].subs);
// Populate sub nodes
for (j=0;j<c_nodes[i].degree;j++){
c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
}
}
}
/* determine degree of each v-node */
for(i=0;i<(CodeLength-NumberParityBits+shift);i++){
count=0;
for (j=0;j<max_col_weight;j++) {
if ( H_cols[i+j*NumberRowsHcols] > 0 ) {
count++;
}
}
v_nodes[i].degree = count;
}
for(i=CodeLength-NumberParityBits+shift;i<CodeLength;i++){
count=0;
if (H1){
if(i!=CodeLength-1){
v_nodes[i].degree=2;
} else{
v_nodes[i].degree=1;
}
} else{
for (j=0;j<max_col_weight;j++) {
if ( H_cols[i+j*NumberRowsHcols] > 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;i++) {
// Allocate sub nodes
v_nodes[i].subs = (struct v_sub_node*) calloc(v_nodes[i].degree, sizeof(struct v_sub_node));
assert(v_nodes[i].subs);
// Populate sub nodes
/* index tells which c-nodes this v-node is connected to */
v_nodes[i].initial_value = input[i];
count=0;
for (j=0;j<v_nodes[i].degree;j++) {
if ((H1)&& (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_index<c_nodes[ v_nodes[i].subs[j].index ].degree;c_index++)
if ( c_nodes[ v_nodes[i].subs[j].index ].subs[c_index].index == i ) {
v_nodes[i].subs[j].socket = c_index;
break;
}
/* initialize v-node with received LLR */
if ( dec_type == 1)
v_nodes[i].subs[j].message = fabs(input[i]);
else
v_nodes[i].subs[j].message = phi0( fabs(input[i]) );
if (input[i] < 0)
v_nodes[i].subs[j].sign = 1;
}
}
/* now finish setting up the c_nodes */
for (i=0;i<NumberParityBits;i++) {
/* index tells which v-nodes this c-node is connected to */
for (j=0;j<c_nodes[i].degree;j++) {
/* search the connected v-node for the proper message value */
for (v_index=0;v_index<v_nodes[ c_nodes[i].subs[j].index ].degree;v_index++)
if (v_nodes[ c_nodes[i].subs[j].index ].subs[v_index].index == i ) {
c_nodes[i].subs[j].socket = v_index;
break;
}
}
}
}
///////////////////////////////////////
/* function for doing the MP decoding */
// Returns the iteration count
int SumProduct( int *parityCheckCount,
char DecodedBits[],
struct c_node c_nodes[],
struct v_node v_nodes[],
int CodeLength,
int NumberParityBits,
int max_iter,
float r_scale_factor,
float q_scale_factor,
int data[] )
{
(void) r_scale_factor;
(void) q_scale_factor;
int result;
int bitErrors;
int i,j, iter;
float phi_sum;
int sign;
float temp_sum;
float Qi;
int ssum;
result = max_iter;
for (iter=0;iter<max_iter;iter++) {
for(i=0; i<CodeLength; i++) DecodedBits[i] = 0; // Clear each pass!
bitErrors = 0;
/* update r */
ssum = 0;
for (j=0;j<NumberParityBits;j++) {
sign = v_nodes[ c_nodes[j].subs[0].index ].subs[ c_nodes[j].subs[0].socket ].sign;
phi_sum = v_nodes[ c_nodes[j].subs[0].index ].subs[ c_nodes[j].subs[0].socket ].message;
for (i=1;i<c_nodes[j].degree;i++) {
// Compiler should optomize this but write the best we can to start from.
struct c_sub_node *cp = &c_nodes[j].subs[i];
struct v_sub_node *vp = &v_nodes[ cp->index ].subs[ cp->socket ];
phi_sum += vp->message;
sign ^= vp->sign;
}
if (sign==0) ssum++;
for (i=0;i<c_nodes[j].degree;i++) {
struct c_sub_node *cp = &c_nodes[j].subs[i];
struct v_sub_node *vp = &v_nodes[ cp->index ].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;i<CodeLength;i++) {
/* first compute the LLR */
Qi = v_nodes[i].initial_value;
for (j=0;j<v_nodes[i].degree;j++) {
struct v_sub_node *vp = &v_nodes[i].subs[j];
Qi += c_nodes[ vp->index ].subs[ vp->socket ].message;
}
/* make hard decision */
if (Qi < 0) {
DecodedBits[i] = 1;
}
/* now subtract to get the extrinsic information */
for (j=0;j<v_nodes[i].degree;j++) {
struct v_sub_node *vp = &v_nodes[i].subs[j];
temp_sum = Qi - c_nodes[ vp->index ].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;i<CodeLength-NumberParityBits;i++)
if ( DecodedBits[i] != data[i] )
bitErrors++;
/* Halt if zero errors */
if (bitErrors == 0) {
result = iter + 1;
break;
}
// count the number of PC satisfied and exit if all OK
*parityCheckCount = ssum;
if (ssum==NumberParityBits) {
result = iter + 1;
break;
}
}
return(result);
}
/* Convenience function to call LDPC decoder from C programs */
int run_ldpc_decoder(struct LDPC *ldpc, uint8_t out_char[], float input[], int *parityCheckCount) {
int max_iter, dec_type;
float q_scale_factor, r_scale_factor;
int max_row_weight, max_col_weight;
int CodeLength, NumberParityBits, NumberRowsHcols, shift, H1;
int i;
struct c_node *c_nodes;
struct v_node *v_nodes;
/* default values */
max_iter = ldpc->max_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<CodeLength; i++) DecodedBits[i] = 0;
/* Call function to do the actual decoding */
int iter = SumProduct( parityCheckCount, DecodedBits, c_nodes, v_nodes,
CodeLength, NumberParityBits, max_iter,
r_scale_factor, q_scale_factor, data_int );
for (i=0; i<CodeLength; i++) out_char[i] = DecodedBits[i];
/* Clean up memory */
free(DecodedBits);
free( data_int );
for (i=0;i<NumberParityBits;i++) {
free( c_nodes[i].subs );
}
free( c_nodes );
for (i=0;i<CodeLength;i++) {
free( v_nodes[i].subs);
}
free( v_nodes );
return iter;
}
void sd_to_llr(float llr[], double sd[], int n) {
double sum, mean, sign, sumsq, estvar, estEsN0, x;
int i;
/* convert SD samples to LLRs -------------------------------*/
sum = 0.0;
for(i=0; i<n; i++)
sum += fabs(sd[i]);
mean = sum/n;
/* find variance from +/-1 symbol position */
sum = sumsq = 0.0;
for(i=0; i<n; i++) {
sign = (sd[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<n; i++)
llr[i] = 4.0L * estEsN0 * sd[i];
}
/*
Determine symbol likelihood from received QPSK symbols.
Notes:
1) We assume fading[] is real, it is also possible to compute
with complex fading, see CML library Demod2D.c source code.
2) Using floats instead of doubles, for stm32.
Testing shows good BERs with floats.
*/
void Demod2D(float symbol_likelihood[], /* output, M*number_symbols */
COMP r[], /* received QPSK symbols, number_symbols */
COMP S_matrix[], /* constellation of size M */
float EsNo,
float fading[], /* real fading values, number_symbols */
float mean_amp,
int number_symbols)
{
int M=QPSK_CONSTELLATION_SIZE;
int i,j;
float tempsr, tempsi, Er, Ei;
/* determine output */
for (i=0;i<number_symbols;i++) { /* go through each received symbol */
for (j=0;j<M;j++) { /* each postulated symbol */
tempsr = fading[i]*S_matrix[j].real/mean_amp;
tempsi = fading[i]*S_matrix[j].imag/mean_amp;
Er = r[i].real/mean_amp - tempsr;
Ei = r[i].imag/mean_amp - tempsi;
symbol_likelihood[i*M+j] = -EsNo*(Er*Er+Ei*Ei);
//printf("symbol_likelihood[%d][%d] = %f\n", i,j,symbol_likelihood[i*M+j]);
}
//exit(0);
}
}
void Somap(float bit_likelihood[], /* number_bits, bps*number_symbols */
float symbol_likelihood[], /* M*number_symbols */
int number_symbols)
{
int M=QPSK_CONSTELLATION_SIZE, bps = QPSK_BITS_PER_SYMBOL;
int n,i,j,k,mask;
float *num = new float[bps];
float *den = new float[bps];
float metric;
for (n=0; n<number_symbols; n++) { /* loop over symbols */
for (k=0;k<bps;k++) {
/* initialize */
num[k] = -1000000;
den[k] = -1000000;
}
for (i=0;i<M;i++) {
metric = symbol_likelihood[n*M+i]; /* channel metric for this symbol */
mask = 1 << (bps - 1);
for (j=0;j<bps;j++) {
mask = mask >> 1;
}
mask = 1 << (bps - 1);
for (k=0;k<bps;k++) { /* loop over bits */
if (mask&i) {
/* this bit is a one */
num[k] = max_star0( num[k], metric );
} else {
/* this bit is a zero */
den[k] = max_star0( den[k], metric );
}
mask = mask >> 1;
}
}
for (k=0;k<bps;k++) {
bit_likelihood[bps*n+k] = num[k] - den[k];
}
}
delete[] den;
delete[] num;
}
void symbols_to_llrs(float llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, float mean_amp, int nsyms) {
int i;
float *symbol_likelihood = new float[nsyms*QPSK_CONSTELLATION_SIZE];
float *bit_likelihood = new float[nsyms*QPSK_BITS_PER_SYMBOL];
Demod2D(symbol_likelihood, rx_qpsk_symbols, S_matrix, EsNo, rx_amps, mean_amp, nsyms);
Somap(bit_likelihood, symbol_likelihood, nsyms);
for(i=0; i<nsyms*QPSK_BITS_PER_SYMBOL; i++) {
llr[i] = -bit_likelihood[i];
}
delete[] bit_likelihood;
delete[] symbol_likelihood;
}
void ldpc_print_info(struct LDPC *ldpc) {
fprintf(stderr, "ldpc->max_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: */

1752
libfreedv/ofdm.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

4
libfreedv/ofdm_internal.h
View File

@ -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<float>{cos(value), sin(value)})
#define cmplxconj(value) (std::complex<float>{cos(value), -sin(value)})
namespace FreeDV
{

223
libfreedv/phi0.cpp Normal file
View File

@ -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 <stdint.h>
#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

11
libfreedv/phi0.h Normal file
View File

@ -0,0 +1,11 @@
// phi0.h
#ifndef PHI0_H
#define PHI0_H
namespace FreeDV {
extern float phi0( float xf );
} // FreeDV
#endif