From ea45dca9d999d99c5a8f101342c234b6e4af0331 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Mon, 8 Jun 2015 11:03:11 +0000 Subject: [PATCH] soft decision support, steps towards error masking, Eb/No estimation in octave and C git-svn-id: https://svn.code.sf.net/p/freetel/code@2176 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/cohpsk.m | 52 ++++++++++---------- codec2-dev/octave/cohpsk_demod_plot.m | 41 +++++++++++++--- codec2-dev/octave/tcohpsk.m | 34 ++++++++++---- codec2-dev/src/c2dec.c | 8 ++-- codec2-dev/src/c2enc.c | 34 +++++++------- codec2-dev/src/codec2.c | 52 ++++++++++++++------ codec2-dev/src/codec2.h | 38 +++++---------- codec2-dev/src/cohpsk.c | 68 +++++++++++++++------------ codec2-dev/src/cohpsk_internal.h | 4 ++ codec2-dev/src/cohpsk_mod.c | 7 ++- codec2-dev/unittest/tcohpsk.c | 13 +++-- codec2-dev/unittest/test_cohpsk_ch.c | 5 +- 12 files changed, 216 insertions(+), 140 deletions(-) diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index 06f73855..fdeb9357 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -225,7 +225,7 @@ end % Symbol rate processing for rx side (demodulator) ------------------------------------------------------- -function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_to_bits(cohpsk, ct_symb_buf) +function [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms cohpsk] = qpsk_symbols_to_bits(cohpsk, ct_symb_buf) framesize = cohpsk.framesize; Nsymb = cohpsk.Nsymb; Nsymbrow = cohpsk.Nsymbrow; @@ -269,7 +269,7 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_ % now correct phase of data symbols rx_symb = zeros(Nsymbrow, Nc); - rx_symb_linear = zeros(1, Nsymbrow*Nc); + rx_symb_linear = zeros(1, Nsymbrow*Nc*Nd); rx_bits = zeros(1, framesize); for c=1:Nc*Nd for r=1:Nsymbrow @@ -279,6 +279,7 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_ rx_symb(r,c) = ct_symb_buf(2+r,c); end i = (c-1)*Nsymbrow + r; + rx_symb_linear(i) = rx_symb(r,c); %printf("phi_ %d %d %f %f\n", r,c,real(exp(-j*phi_(r,c))), imag(exp(-j*phi_(r,c)))); end end @@ -292,7 +293,6 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_ for d=1:Nd-1 div_symb += rx_symb(r,c + Nc*d); end - rx_symb_linear(i) = div_symb; rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(div_symb); end end @@ -302,34 +302,32 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_ % position. However this is complicated by fading, which means the % amplitude of the symbols is constantly changing. - % Now the scatter diagram in a fading channel is a X shape. The - % noise can be resolved into two components at right angles to - % each other. The component along the the "thickness" of the arms - % is proportional to the noise power and not affected by fading. + % Now the scatter diagram in a fading channel is a X or cross + % shape. The noise can be resolved into two components at right + % angles to each other. The component along the the "thickness" + % of the arms is proportional to the noise power and not affected + % by fading. We only use points further along the real axis than + % the mean amplitude so we keep out of the central nosiey blob - v = zeros(1,Nsymb); - for i=1:Nsymb + sig_rms = mean(abs(rx_symb_linear)); + + sum_x = 0; + sum_xx = 0; + n = 0; + for i=1:Nsymb*Nd s = rx_symb_linear(i); - if abs(real(s)) > abs(imag(s)) - v(i) = imag(s); - else - v(i) = real(s); + if abs(real(s)) > sig_rms + sum_x += imag(s); + sum_xx += imag(s)*imag(s); + n++; end - %printf("s: %f %f v: %f\n", real(s), imag(s), v(i)); end - - % Note we are only measuring variance in one axis, as other axis is obscured by fading. We assume - % that total noise power is shared between both axis so multiply by sqrt(2) to get an estimate of - % total noise pwr. Small constant prevents divide by zero errors on start up. - - No_ = var(v)*sqrt(2) + 1E-6; - - % Estimate signal power - - Es_ = mean(amp_ .^ 2); - - EsNo_ = Es_/No_; - %printf("Es_: %f No_: %f Es/No: %f Es/No dB: %f\n", Es_, No_, Es_/No_, 10*log10(EsNo_)); + + noise_var = 0; + if n > 1 + noise_var = (n*sum_xx - sum_x*sum_x)/(n*(n-1)); + end + noise_rms = sqrt(noise_var); endfunction diff --git a/codec2-dev/octave/cohpsk_demod_plot.m b/codec2-dev/octave/cohpsk_demod_plot.m index 2d4ab30d..44376bac 100644 --- a/codec2-dev/octave/cohpsk_demod_plot.m +++ b/codec2-dev/octave/cohpsk_demod_plot.m @@ -4,9 +4,9 @@ % Plot Octave outputs from cohpsk_demod, c2dec, to visualise whats going on % when errors hit the system -% $ ./c2enc 700 ../../raw/ve9qrp_10s.raw - | ./cohpsk_mod - - | ./cohpsk_ch - - -60 50 1 1 | ./cohpsk_demod - - | ./c2dec 700 - - --dump ve9qrp | play -t raw -r 8000 -s -2 - -q +% $ ./c2enc 700 ../../raw/ve9qrp_10s.raw - | ./cohpsk_mod - - | ./cohpsk_ch - - -60 50 1 1 | ./cohpsk_demod - - cohpsk_demod.txt | ./c2dec 700 - - --dump ve9qrp | play -t raw -r 8000 -s -2 - -q -% ./c2enc 700 ../../raw/ve9qrp_10s.raw - | ./cohpsk_mod - - | ./cohpsk_ch - - -30 50 1 1 | ./cohpsk_demod - - | ./c2dec 700 - - --dump ve9qrp_snr3 | play -t raw -r 8000 -s -2 - -q +% ./c2enc 700 ../../raw/ve9qrp_10s.raw - | ./cohpsk_mod - - | ./cohpsk_ch - - -30 50 1 1 | ./cohpsk_demod - - cohpsk_demod.txt | ./c2dec 700 - - --dump ve9qrp_snr3 | play -t raw -r 8000 -s -2 - -q graphics_toolkit ("gnuplot"); @@ -110,9 +110,10 @@ hold off figure(8) clf; -f1=fft(ve9qrp_ak_(st:en,:)',128); -f2=fft(ve9qrp_snr3_ak_(1:en-st+1,:)',128); -d = (20*log10(abs(f1)) - 20*log10(abs(f2))); +f1=1./fft(ve9qrp_ak_(st:en,:)',128); +f2=1./fft(ve9qrp_snr3_ak_(1:en-st+1,:)',128); +%d = (20*log10(abs(f1)) - 20*log10(abs(f2))); +d = 20*log10(abs(f2)); sdsq = mean(d.^2); plot(sdsq) title('spectral distortion clean and channel SNR=3dB') @@ -122,9 +123,37 @@ clf; y = 1:en-st+1; x = 1:40; %mesh(y,x,-20*log10(abs(f2(1:40,:)))); -mesh(y,x,d(1:40,:)); +mesh(y,x,d(x,:)); grid title('Synthesis filter difference between clean and channel SNR=3dB'); xlabel('Time (codec frames)') ylabel('Frequency 0 to 2500Hz'); zlabel('Difference (dB)'); + +% map soft decn information to LSPs + +mel1 = ve9qrp_snr3_softdec(:,10:12); +mel2 = ve9qrp_snr3_softdec(:,13:14); +mel3 = ve9qrp_snr3_softdec(:,15:18); +mel4 = ve9qrp_snr3_softdec(:,19:21); +mel5 = ve9qrp_snr3_softdec(:,22:24); +mel6 = ve9qrp_snr3_softdec(:,25:26); +softdec_mel = [sum(mel1'.^2); sum(mel2'.^2); sum(mel3'.^2); sum(mel4'.^2); sum(mel5'.^2); sum(mel6'.^2)]; + +figure(10) +clf; +y = 1:en-st+1; +x = 1:6; +%mesh(y,x,-20*log10(abs(f2(1:40,:)))); +mesh(y, x, softdec_mel(:,y)); +grid +xlabel('Codec frame') +ylabel('LSP') +zlabel('Power') +%axis([1 (en-st+1) 1 6 -10 5]) + +% plot symbol energy against SD + +figure(11) +semilogx(mean(softdec_mel(:,1:en-st+1)), sdsq,'+') +grid diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 93eadb94..51e1e7a5 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -41,7 +41,7 @@ % + Only a small change in fading perf with filter on/off % + however other filters may have other effects, should test this, % e.g. scatter plots, some sort of BER metric? -% [ ] pilot based EsNo estimation +% [X] EsNo estimation % [ ] filter reqd with compression? % + make sure not too much noise passed into noise floor % [X] different diversity combination @@ -51,8 +51,6 @@ % + ssb filter % + compression % + make sure it's flat with many errors -% [ ] linear tracking of ampl ests -% + seem to only change every graphics_toolkit ("gnuplot"); more off; @@ -66,8 +64,8 @@ randn('state',1); % select which test ---------------------------------------------------------- -%test = 'compare to c'; -test = 'awgn'; +test = 'compare to c'; +%test = 'awgn'; %test = 'fading'; % some parameters that can be over ridden, e.g. to disable parts of modem @@ -94,10 +92,10 @@ end % should be BER around 0.015 to 0.02 if strcmp(test, 'awgn') - frames = 10; + frames = 100; foff = 0; dfoff = -0/Fs; - EsNodB = 80; + EsNodB = 8; fading_en = 0; hf_delay_ms = 2; compare_with_c = 0; @@ -107,10 +105,10 @@ end % Similar to AWGN - should be BER around 0.015 to 0.02 if strcmp(test, 'fading'); - frames = 200; + frames = 100; foff = -10.5; dfoff = 0.0/Fs; - EsNodB = 12; + EsNodB = 9; fading_en = 1; hf_delay_ms = 2; compare_with_c = 0; @@ -229,6 +227,9 @@ rx_timing_log = []; ratio_log = []; foff_log = []; f_est_log = []; +sig_rms_log = []; +noise_rms_log = []; +noise_rms_filt_log = []; % Channel modeling and BER measurement ---------------------------------------- @@ -259,6 +260,8 @@ ch_fdm_delay = zeros(1, acohpsk.Nsymbrowpilot*M + nhfdelay); h = freqz(b,a,(600:2600)/(Fs/(2*pi))); filt_gain = (2600-600)/sum(abs(h) .^ 2); % ensures power after filter == before filter +noise_rms_filt = 0; + % main loop -------------------------------------------------------------------- % run mod and channel as aseparate loop so we can resample to simulate sample rate differences @@ -479,7 +482,7 @@ for f=1:frames; % if we are in sync complete demodulation with symbol rate processing if (next_sync == 1) || (sync == 1) - [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_] = qpsk_symbols_to_bits(acohpsk, acohpsk.ct_symb_ff_buf); + [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms] = qpsk_symbols_to_bits(acohpsk, acohpsk.ct_symb_ff_buf); rx_symb_log = [rx_symb_log; rx_symb]; rx_amp_log = [rx_amp_log; amp_]; rx_phi_log = [rx_phi_log; phi_]; @@ -487,6 +490,10 @@ for f=1:frames; tx_bits_prev_log = [tx_bits_prev_log prev_tx_bits2]; ratio_log = [ratio_log acohpsk.ratio]; ct_symb_ff_log = [ct_symb_ff_log; acohpsk.ct_symb_ff_buf(1:acohpsk.Nsymbrowpilot,:)]; + sig_rms_log = [sig_rms_log sig_rms]; + noise_rms_log = [noise_rms_log noise_rms]; + noise_rms_filt = 0.9*noise_rms_filt + 0.1*noise_rms; + noise_rms_filt_log = [noise_rms_filt_log noise_rms_filt]; % BER stats @@ -594,6 +601,8 @@ if compare_with_c check(rx_timing_log, rx_timing_log_c, 'rx_timing',0.001); check(rx_bits_log, rx_bits_log_c, 'rx_bits'); check(f_est_log, f_est_log_c, 'f_est'); + check(sig_rms_log, sig_rms_log_c, 'sig_rms'); + check(noise_rms_log, noise_rms_log_c, 'noise_rms'); else @@ -653,8 +662,12 @@ else figure(6) clf; + subplot(211) stem(nerr_log) title('Bit Errors'); + subplot(212) + plot(noise_rms_filt_log,'r', sig_rms_log,'g'); + title('Est rms signal and noise') figure(7); clf; @@ -682,6 +695,7 @@ else plot(error_positions_hist) title('histogram of bit errors') + end diff --git a/codec2-dev/src/c2dec.c b/codec2-dev/src/c2dec.c index cb20bd9d..8d186382 100644 --- a/codec2-dev/src/c2dec.c +++ b/codec2-dev/src/c2dec.c @@ -116,7 +116,7 @@ int main(int argc, char *argv[]) ber = 0.0; burst_length = burst_period = 0.0; burst_timer = 0.0; - dump = natural = 0; + dump = natural = softdec = 0; codec2 = codec2_create(mode); nsam = codec2_samples_per_frame(codec2); @@ -124,7 +124,7 @@ int main(int argc, char *argv[]) buf = (short*)malloc(nsam*sizeof(short)); nbyte = (nbit + 7) / 8; bits = (unsigned char*)malloc(nbyte*sizeof(char)); - softdec_bits = (int*)malloc(nbit*sizeof(float)); + softdec_bits = (float*)malloc(nbit*sizeof(float)); frames = bit_errors = bits_proc = 0; nstart_bit = 0; nend_bit = nbit-1; @@ -174,6 +174,7 @@ int main(int argc, char *argv[]) codec2_set_natural_or_gray(codec2, !natural); //printf("%d %d\n", nstart_bit, nend_bit); + //fprintf(stderr, "softdec: %d natural: %d\n", softdec, natural); if (softdec) ret = (fread(softdec_bits, sizeof(float), nbit, fin) == (size_t)nbit); else @@ -263,7 +264,8 @@ int main(int argc, char *argv[]) byte++; } } - } + codec2_set_softdec(codec2, softdec_bits); + } codec2_decode_ber(codec2, buf, bits, ber_est); fwrite(buf, sizeof(short), nsam, fout); diff --git a/codec2-dev/src/c2enc.c b/codec2-dev/src/c2enc.c index 5f5fddab..5a000b86 100644 --- a/codec2-dev/src/c2enc.c +++ b/codec2-dev/src/c2enc.c @@ -41,12 +41,12 @@ int main(int argc, char *argv[]) FILE *fout; short *buf; unsigned char *bits; - int nsam, nbit, nbyte, gray, unpacked; - int *unpacked_bits; + int nsam, nbit, nbyte, gray, softdec; + float *unpacked_bits; int bit, byte,i; if (argc < 4) { - printf("usage: c2enc 3200|2400|1600|1400|1300|1200|700 InputRawspeechFile OutputBitFile [--natural] [--unpacked]\n"); + printf("usage: c2enc 3200|2400|1600|1400|1300|1200|700 InputRawspeechFile OutputBitFile [--natural] [--softdec]\n"); printf("e.g c2enc 1400 ../raw/hts1a.raw hts1a.c2\n"); printf("e.g c2enc 1300 ../raw/hts1a.raw hts1a.c2 --natural\n"); exit(1); @@ -92,38 +92,38 @@ int main(int argc, char *argv[]) nbyte = (nbit + 7) / 8; bits = (unsigned char*)malloc(nbyte*sizeof(char)); - unpacked_bits = (int*)malloc(nbit*sizeof(int)); + unpacked_bits = (float*)malloc(nbit*sizeof(float)); - if (argc == 5) { - if (strcmp(argv[4], "--natural") == 0) + gray = 1; + softdec = 0; + for (i=4; i> bit) & 0x1; + unpacked_bits[i] = 1.0 - 2.0*((bits[byte] >> bit) & 0x1); bit--; if (bit < 0) { bit = 7; byte++; } } - fwrite(unpacked_bits, sizeof(int), nbit, fout); + fwrite(unpacked_bits, sizeof(float), nbit, fout); } else fwrite(bits, sizeof(char), nbyte, fout); diff --git a/codec2-dev/src/codec2.c b/codec2-dev/src/codec2.c index ad1f986f..d1208641 100644 --- a/codec2-dev/src/codec2.c +++ b/codec2-dev/src/codec2.c @@ -92,7 +92,7 @@ static void ear_protection(float in_out[], int n); \*---------------------------------------------------------------------------*/ -struct CODEC2 * CODEC2_WIN32SUPPORT codec2_create(int mode) +struct CODEC2 * codec2_create(int mode) { struct CODEC2 *c2; int i,l; @@ -171,7 +171,7 @@ struct CODEC2 * CODEC2_WIN32SUPPORT codec2_create(int mode) \*---------------------------------------------------------------------------*/ -void CODEC2_WIN32SUPPORT codec2_destroy(struct CODEC2 *c2) +void codec2_destroy(struct CODEC2 *c2) { assert(c2 != NULL); free(c2->bpf_buf); @@ -191,7 +191,7 @@ void CODEC2_WIN32SUPPORT codec2_destroy(struct CODEC2 *c2) \*---------------------------------------------------------------------------*/ -int CODEC2_WIN32SUPPORT codec2_bits_per_frame(struct CODEC2 *c2) { +int codec2_bits_per_frame(struct CODEC2 *c2) { if (c2->mode == CODEC2_MODE_3200) return 64; if (c2->mode == CODEC2_MODE_2400) @@ -221,7 +221,7 @@ int CODEC2_WIN32SUPPORT codec2_bits_per_frame(struct CODEC2 *c2) { \*---------------------------------------------------------------------------*/ -int CODEC2_WIN32SUPPORT codec2_samples_per_frame(struct CODEC2 *c2) { +int codec2_samples_per_frame(struct CODEC2 *c2) { if (c2->mode == CODEC2_MODE_3200) return 160; if (c2->mode == CODEC2_MODE_2400) @@ -240,7 +240,7 @@ int CODEC2_WIN32SUPPORT codec2_samples_per_frame(struct CODEC2 *c2) { return 0; /* shouldnt get here */ } -void CODEC2_WIN32SUPPORT codec2_encode(struct CODEC2 *c2, unsigned char *bits, short speech[]) +void codec2_encode(struct CODEC2 *c2, unsigned char *bits, short speech[]) { assert(c2 != NULL); assert( @@ -269,12 +269,12 @@ void CODEC2_WIN32SUPPORT codec2_encode(struct CODEC2 *c2, unsigned char *bits, s codec2_encode_700(c2, bits, speech); } -void CODEC2_WIN32SUPPORT codec2_decode(struct CODEC2 *c2, short speech[], const unsigned char *bits) +void codec2_decode(struct CODEC2 *c2, short speech[], const unsigned char *bits) { codec2_decode_ber(c2, speech, bits, 0.0); } -void CODEC2_WIN32SUPPORT codec2_decode_ber(struct CODEC2 *c2, short speech[], const unsigned char *bits, float ber_est) +void codec2_decode_ber(struct CODEC2 *c2, short speech[], const unsigned char *bits, float ber_est) { assert(c2 != NULL); assert( @@ -1404,9 +1404,9 @@ void codec2_encode_700(struct CODEC2 *c2, unsigned char * bits, short speech[]) encode_mels_scalar(indexes, mel, LPC_ORD_LOW); for(i=0; igray); } - + pack_natural_or_gray(bits, &nbit, spare, 2, c2->gray); assert(nbit == (unsigned)codec2_bits_per_frame(c2)); @@ -1459,8 +1459,9 @@ void codec2_decode_700(struct CODEC2 *c2, short speech[], const unsigned char * e[3] = decode_energy(e_index, 3); for(i=0; igray); } + decode_mels_scalar(mel, indexes, LPC_ORD_LOW); for(i=0; isoftdec) { + float e = 0.0; + for(i=9; i<9+17; i++) + e += c2->softdec[i]*c2->softdec[i]; + e /= 6.0; + //fprintf(stderr, "e: %f\n", e); + //if (e < 0.3) + // bw_expand_lsps(&lsps[3][0], LPC_ORD_LOW, 150.0, 300.0); + } + #endif + /* interpolate ------------------------------------------------*/ /* LSPs, Wo, and energy are sampled every 40ms so we interpolate @@ -1644,7 +1660,7 @@ static void ear_protection(float in_out[], int n) { } } -void CODEC2_WIN32SUPPORT codec2_set_lpc_post_filter(struct CODEC2 *c2, int enable, int bass_boost, float beta, float gamma) +void codec2_set_lpc_post_filter(struct CODEC2 *c2, int enable, int bass_boost, float beta, float gamma) { assert((beta >= 0.0) && (beta <= 1.0)); assert((gamma >= 0.0) && (gamma <= 1.0)); @@ -1660,7 +1676,7 @@ void CODEC2_WIN32SUPPORT codec2_set_lpc_post_filter(struct CODEC2 *c2, int enabl Experimental method of sending voice/data frames for FreeDV. */ -int CODEC2_WIN32SUPPORT codec2_get_spare_bit_index(struct CODEC2 *c2) +int codec2_get_spare_bit_index(struct CODEC2 *c2) { assert(c2 != NULL); @@ -1684,7 +1700,7 @@ int CODEC2_WIN32SUPPORT codec2_get_spare_bit_index(struct CODEC2 *c2) for convenience. */ -int CODEC2_WIN32SUPPORT codec2_rebuild_spare_bit(struct CODEC2 *c2, int unpacked_bits[]) +int codec2_rebuild_spare_bit(struct CODEC2 *c2, int unpacked_bits[]) { int v1,v3; @@ -1732,9 +1748,15 @@ int CODEC2_WIN32SUPPORT codec2_rebuild_spare_bit(struct CODEC2 *c2, int unpacked return -1; } -void CODEC2_WIN32SUPPORT codec2_set_natural_or_gray(struct CODEC2 *c2, int gray) +void codec2_set_natural_or_gray(struct CODEC2 *c2, int gray) { assert(c2 != NULL); c2->gray = gray; } +void codec2_set_softdec(struct CODEC2 *c2, float *softdec) +{ + assert(c2 != NULL); + c2->softdec = softdec; +} + diff --git a/codec2-dev/src/codec2.h b/codec2-dev/src/codec2.h index 834d9a74..d3c6d27d 100644 --- a/codec2-dev/src/codec2.h +++ b/codec2-dev/src/codec2.h @@ -33,19 +33,6 @@ #ifndef __CODEC2__ #define __CODEC2__ -/* set up the calling convention for DLL function import/export for - WIN32 cross compiling */ - -#ifdef __CODEC2_WIN32__ -#ifdef __CODEC2_BUILDING_DLL__ -#define CODEC2_WIN32SUPPORT __declspec(dllexport) __stdcall -#else -#define CODEC2_WIN32SUPPORT __declspec(dllimport) __stdcall -#endif -#else -#define CODEC2_WIN32SUPPORT -#endif - #define CODEC2_MODE_3200 0 #define CODEC2_MODE_2400 1 #define CODEC2_MODE_1600 2 @@ -56,18 +43,19 @@ struct CODEC2; -struct CODEC2 * CODEC2_WIN32SUPPORT codec2_create(int mode); -void CODEC2_WIN32SUPPORT codec2_destroy(struct CODEC2 *codec2_state); -void CODEC2_WIN32SUPPORT codec2_encode(struct CODEC2 *codec2_state, unsigned char * bits, short speech_in[]); -void CODEC2_WIN32SUPPORT codec2_decode(struct CODEC2 *codec2_state, short speech_out[], const unsigned char *bits); -void CODEC2_WIN32SUPPORT codec2_decode_ber(struct CODEC2 *codec2_state, short speech_out[], const unsigned char *bits, float ber_est); -int CODEC2_WIN32SUPPORT codec2_samples_per_frame(struct CODEC2 *codec2_state); -int CODEC2_WIN32SUPPORT codec2_bits_per_frame(struct CODEC2 *codec2_state); - -void CODEC2_WIN32SUPPORT codec2_set_lpc_post_filter(struct CODEC2 *codec2_state, int enable, int bass_boost, float beta, float gamma); -int CODEC2_WIN32SUPPORT codec2_get_spare_bit_index(struct CODEC2 *codec2_state); -int CODEC2_WIN32SUPPORT codec2_rebuild_spare_bit(struct CODEC2 *codec2_state, int unpacked_bits[]); -void CODEC2_WIN32SUPPORT codec2_set_natural_or_gray(struct CODEC2 *codec2_state, int gray); +struct CODEC2 * codec2_create(int mode); +void codec2_destroy(struct CODEC2 *codec2_state); +void codec2_encode(struct CODEC2 *codec2_state, unsigned char * bits, short speech_in[]); +void codec2_decode(struct CODEC2 *codec2_state, short speech_out[], const unsigned char *bits); +void codec2_decode_ber(struct CODEC2 *codec2_state, short speech_out[], const unsigned char *bits, float ber_est); +int codec2_samples_per_frame(struct CODEC2 *codec2_state); +int codec2_bits_per_frame(struct CODEC2 *codec2_state); + +void codec2_set_lpc_post_filter(struct CODEC2 *codec2_state, int enable, int bass_boost, float beta, float gamma); +int codec2_get_spare_bit_index(struct CODEC2 *codec2_state); +int codec2_rebuild_spare_bit(struct CODEC2 *codec2_state, int unpacked_bits[]); +void codec2_set_natural_or_gray(struct CODEC2 *codec2_state, int gray); +void codec2_set_softdec(struct CODEC2 *c2, float *softdec); #endif diff --git a/codec2-dev/src/cohpsk.c b/codec2-dev/src/cohpsk.c index fbbc91f8..a87203ea 100644 --- a/codec2-dev/src/cohpsk.c +++ b/codec2-dev/src/cohpsk.c @@ -155,6 +155,8 @@ struct COHPSK *cohpsk_create(void) coh->fdmdv = fdmdv; + coh->sig_rms = coh->noise_rms = 0.0; + /* disable optional logging by default */ coh->rx_baseband_log = NULL; @@ -268,34 +270,19 @@ void bits_to_qpsk_symbols(COMP tx_symb[][COHPSK_NC*ND], int tx_bits[], int nbits void qpsk_symbols_to_bits(struct COHPSK *coh, float rx_bits[], COMP ct_symb_buf[][COHPSK_NC*ND]) { - int p, r, c, i, pc, d; + int p, r, c, i, pc, d, n; float x[NPILOTSFRAME+2], x1; COMP y[NPILOTSFRAME+2], yfit; + COMP rx_symb_linear[NSYMROW*COHPSK_NC*ND]; COMP m, b; COMP __attribute__((unused)) corr, rot, pi_on_4, phi_rect, div_symb; float mag, __attribute__((unused)) phi_, __attribute__((unused)) amp_; + float sum_x, sum_xx, noise_var; + COMP s; pi_on_4.real = cosf(M_PI/4); pi_on_4.imag = sinf(M_PI/4); for(c=0; cpilot2[p][c], ct_symb_buf[sampling_points[p]][c])); - mag += cabsolute(ct_symb_buf[sampling_points[p]][c]); - } - - phi_ = atan2f(corr.imag, corr.real); - amp_ = mag/(NPILOTSFRAME+2); - for(r=0; rphi_[r][c] = phi_; - coh->amp_[r][c] = amp_; - } -#else /* Set up lin reg model and interpolate phase. Works better than average for channels with quickly changing phase, like HF. */ @@ -304,17 +291,15 @@ void qpsk_symbols_to_bits(struct COHPSK *coh, float rx_bits[], COMP ct_symb_buf[ x[p] = sampling_points[p]; pc = c % COHPSK_NC; y[p] = fcmult(coh->pilot2[p][pc], ct_symb_buf[sampling_points[p]][c]); - //printf("%f %f\n", y[p].real, y[p].imag); } - //printf("\n"); + linreg(&m, &b, x, y, NPILOTSFRAME+2); for(r=0; rphi_[r][c] = atan2(yfit.imag, yfit.real); - //printf(" %f", coh->phi_[r][c]); } - //printf("\n"); + /* amplitude estimation */ mag = 0.0; @@ -325,22 +310,20 @@ void qpsk_symbols_to_bits(struct COHPSK *coh, float rx_bits[], COMP ct_symb_buf[ for(r=0; ramp_[r][c] = amp_; } -#endif } - //exit(0); + /* now correct phase of data symbols */ for(c=0; cphi_[r][c]); phi_rect.imag = -sinf(coh->phi_[r][c]); coh->rx_symb[r][c] = cmult(ct_symb_buf[NPILOTSFRAME + r][c], phi_rect); - //printf("%d %d %f %f\n", r,c, coh->rx_symb[r][c].real, coh->rx_symb[r][c].imag); - //printf("phi_ %d %d %f %f\n", r,c, ct_symb_buf[NPILOTSFRAME + r][c].real, ct_symb_buf[NPILOTSFRAME + r][c].imag); + i = c*NSYMROW + r; + rx_symb_linear[i] = coh->rx_symb[r][c]; } } - /* and finally optional diversity combination and make decn on bits */ + /* and finally optional diversity combination, note output is soft decm a "1" is < 0 */ for(c=0; csig_rms = mag/(NSYMROW*COHPSK_NC*ND); + + sum_x = 0; + sum_xx = 0; + n = 0; + for (i=0; i coh->sig_rms) { + sum_x += s.imag; + sum_xx += s.imag*s.imag; + n++; + } + } + + noise_var = 0; + if (n > 1) { + noise_var = (n*sum_xx - sum_x*sum_x)/(n*(n-1)); + } + coh->noise_rms = sqrtf(noise_var); + } diff --git a/codec2-dev/src/cohpsk_internal.h b/codec2-dev/src/cohpsk_internal.h index 5efc222a..6ffd25bb 100644 --- a/codec2-dev/src/cohpsk_internal.h +++ b/codec2-dev/src/cohpsk_internal.h @@ -65,8 +65,12 @@ struct COHPSK { int frame; float ratio; + float sig_rms; + float noise_rms; + struct FDMDV *fdmdv; + /* optional log variables used for testing Octave to C port */ COMP *rx_baseband_log; diff --git a/codec2-dev/src/cohpsk_mod.c b/codec2-dev/src/cohpsk_mod.c index 1cbe375a..5220f826 100644 --- a/codec2-dev/src/cohpsk_mod.c +++ b/codec2-dev/src/cohpsk_mod.c @@ -4,7 +4,7 @@ AUTHOR......: David Rowe DATE CREATED: April 5 2015 - Given an input file of bits (note one bit per int, not compressed), + Given an input file of bits (note one bit per float, soft decision format), outputs a raw file (8kHz, 16 bit shorts) of COHPSK modem samples ready to send over a HF radio channel. @@ -41,6 +41,7 @@ int main(int argc, char *argv[]) { FILE *fin, *fout; struct COHPSK *cohpsk; + float tx_bits_sd[COHPSK_BITS_PER_FRAME]; int tx_bits[COHPSK_BITS_PER_FRAME]; COMP tx_fdm[COHPSK_SAMPLES_PER_FRAME]; short tx_fdm_scaled[COHPSK_SAMPLES_PER_FRAME]; @@ -70,9 +71,11 @@ int main(int argc, char *argv[]) frames = 0; - while(fread(tx_bits, sizeof(int), COHPSK_BITS_PER_FRAME, fin) == COHPSK_BITS_PER_FRAME) { + while(fread(tx_bits_sd, sizeof(float), COHPSK_BITS_PER_FRAME, fin) == COHPSK_BITS_PER_FRAME) { frames++; + for(i=0; i 0.0; + rx_bits[i] = rx_bits_sd[i] < 0.0; ch_fdm_frame_log_index += tmp; @@ -241,7 +241,10 @@ int main(int argc, char *argv[]) } memcpy(&rx_bits_log[COHPSK_BITS_PER_FRAME*log_bits], rx_bits, sizeof(int)*COHPSK_BITS_PER_FRAME); log_bits++; - f_est_log[f_est_samples++] = coh->f_est; + f_est_log[f_est_samples] = coh->f_est; + sig_rms_log[f_est_samples] = coh->sig_rms; + noise_rms_log[f_est_samples] = coh->noise_rms; + f_est_samples++;; } assert(log_r <= NSYMROWPILOT*FRAMES); @@ -275,6 +278,8 @@ int main(int argc, char *argv[]) octave_save_complex(fout, "rx_symb_log_c", (COMP*)rx_symb_log, log_data_r, COHPSK_NC*ND, COHPSK_NC*ND); octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1, COHPSK_BITS_PER_FRAME*log_bits); octave_save_float(fout, "f_est_log_c", &f_est_log[1], 1, f_est_samples-1, f_est_samples-1); + octave_save_float(fout, "sig_rms_log_c", sig_rms_log, 1, f_est_samples, f_est_samples-1); + octave_save_float(fout, "noise_rms_log_c", noise_rms_log, 1, f_est_samples, f_est_samples); #ifdef XX #endif fclose(fout); diff --git a/codec2-dev/unittest/test_cohpsk_ch.c b/codec2-dev/unittest/test_cohpsk_ch.c index 32f40f82..8ba27868 100644 --- a/codec2-dev/unittest/test_cohpsk_ch.c +++ b/codec2-dev/unittest/test_cohpsk_ch.c @@ -55,7 +55,8 @@ tcohpsk first (any variant) to load the function into Octave, e.g.: octave:17> tcohpsk - octave:18> write_noise_file("../raw/fading_samples.float", 7500, 7500*60) + octave:18> write_noise_file("../raw/fast_fading_samples.float", 7500, 7500*60) + octave:19> write_noise_file("../raw/slow_fading_samples.float", 75000, 7500*60) */ #define FADING_FILE_NAME "../../raw/fading_samples.float" @@ -89,7 +90,7 @@ int main(int argc, char *argv[]) int ch_buf_n; float tx_pwr, rx_pwr, noise_pwr; int error_positions_hist[COHPSK_BITS_PER_FRAME]; - int log_data_r, c, j, tmp, ret; + int log_data_r, c, j, tmp; for(i=0; i