From: drowe67 Date: Sat, 2 Aug 2014 02:00:09 +0000 (+0000) Subject: improved fdmdv channel simulation function and wrote cmd line util to test it. getti... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=554de3e35938cc4a117570b259e49ed36e33f287;p=freetel-svn-tracking.git improved fdmdv channel simulation function and wrote cmd line util to test it. getting same BER results from Octave and C demods which is a good sign. Some of the fdmdv_demod_c plots a bit messed up but BER results correct git-svn-id: https://svn.code.sf.net/p/freetel/code@1780 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/fdmdv_demod.m b/codec2-dev/octave/fdmdv_demod.m index 65f0acaa..6cc7043f 100644 --- a/codec2-dev/octave/fdmdv_demod.m +++ b/codec2-dev/octave/fdmdv_demod.m @@ -19,7 +19,7 @@ function fdmdv_demod(rawfilename, nbits, NumCarriers, errorpatternfilename, symb frames = nbits/(Nc*Nb); prev_rx_symbols = ones(Nc+1,1); - foff_phase = 1; + foff_phase_rect = 1; % BER stats @@ -40,6 +40,7 @@ function fdmdv_demod(rawfilename, nbits, NumCarriers, errorpatternfilename, symb rx_symbols_log = []; rx_timing_log = []; + foff_coarse_log = []; foff_log = []; rx_fdm_log = []; snr_est_log = []; @@ -106,22 +107,22 @@ function fdmdv_demod(rawfilename, nbits, NumCarriers, errorpatternfilename, symb if sync == 0 foff = foff_coarse; end - foff_log = [ foff_log foff ]; + foff_coarse_log = [foff_coarse_log foff_coarse]; + foff_rect = exp(j*2*pi*foff/Fs); for i=1:nin - foff_phase *= foff_rect'; - rx_fdm(i) = rx_fdm(i)*foff_phase; + foff_phase_rect *= foff_rect'; + rx_fdm_fcorr(i) = rx_fdm(i)*foff_phase_rect; end % baseband processing - rx_baseband = fdm_downconvert(rx_fdm, nin); - rx_filt = rx_filter(rx_baseband, nin); - + rx_fdm_filter = rxdec_filter(rx_fdm_fcorr, nin); + rx_filt = down_convert_and_rx_filter(rx_fdm_filter, nin, M/Q); [rx_symbols rx_timing] = rx_est_timing(rx_filt, nin); - rx_timing_log = [rx_timing_log rx_timing]; + nin = M; if rx_timing > 2*M/P nin += M/P; @@ -130,11 +131,7 @@ function fdmdv_demod(rawfilename, nbits, NumCarriers, errorpatternfilename, symb nin -= M/P; end - if strcmp(modulation,'dqpsk') - rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols./abs(prev_rx_symbols))*exp(j*pi/4)]; - else - rx_symbols_log = [rx_symbols_log rx_symbols]; - endif + rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols./abs(prev_rx_symbols))*exp(j*pi/4)]; [rx_bits sync_bit f_err pd] = psk_to_bits(prev_rx_symbols, rx_symbols, modulation); % optionally save output symbols @@ -167,11 +164,12 @@ function fdmdv_demod(rawfilename, nbits, NumCarriers, errorpatternfilename, symb % update some states + prev_rx_symbols = rx_symbols; [sig_est noise_est] = snr_update(sig_est, noise_est, pd); snr_est = calc_snr(sig_est, noise_est); snr_est_log = [snr_est_log snr_est]; foff -= 0.5*f_err; - prev_rx_symbols = rx_symbols; + foff_log = [foff_log foff]; % freq est state machine diff --git a/codec2-dev/octave/fdmdv_demod_c.m b/codec2-dev/octave/fdmdv_demod_c.m index 98b571b3..b5c75334 100644 --- a/codec2-dev/octave/fdmdv_demod_c.m +++ b/codec2-dev/octave/fdmdv_demod_c.m @@ -9,8 +9,8 @@ % Version 2 % -function fdmdv_demod_c(dumpfilename, bits) - NumCarriers = 16; +function fdmdv_demod_c(dumpfilename, bits, NumCarriers) + fdmdv; % include modem code frames = bits/(Nc*Nb); @@ -69,6 +69,9 @@ function fdmdv_demod_c(dumpfilename, bits) test_frame_sync_log = [test_frame_sync_log test_frame_sync_state]; end + ber = total_bit_errors / total_bits; + printf("%d bits %d errors BER: %1.4f\n",total_bits, total_bit_errors, ber); + % --------------------------------------------------------------------- % Plots % --------------------------------------------------------------------- diff --git a/codec2-dev/octave/tfdmdv.m b/codec2-dev/octave/tfdmdv.m index bba5a7ab..0b9bc846 100644 --- a/codec2-dev/octave/tfdmdv.m +++ b/codec2-dev/octave/tfdmdv.m @@ -104,7 +104,6 @@ for f=1:frames fbb_phase_rx /= mag; [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, nin); - [foff_coarse S1 S2] = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin); %sync = 0; % when debugging good idea to uncomment this to "open loop" @@ -147,7 +146,7 @@ for f=1:frames end nin_log = [nin_log nin]; - [rx_bits sync_bit foff_fine pd] = psk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk'); + [rx_bits sync_bit f_err pd] = psk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk'); phase_difference_log = [phase_difference_log pd]; [sig_est noise_est] = snr_update(sig_est, noise_est, pd); @@ -158,7 +157,7 @@ for f=1:frames rx_bits_log = [rx_bits_log rx_bits]; foff_fine_log = [foff_fine_log foff_fine]; sync_bit_log = [sync_bit_log sync_bit]; - foff -= 0.5*foff_fine; + foff -= 0.5*f_err; foff_log = [foff_log foff]; % freq est state machine diff --git a/codec2-dev/src/CMakeLists.txt b/codec2-dev/src/CMakeLists.txt index 97487369..dca5aa0c 100644 --- a/codec2-dev/src/CMakeLists.txt +++ b/codec2-dev/src/CMakeLists.txt @@ -226,6 +226,9 @@ target_link_libraries(fdmdv_demod ${CMAKE_REQUIRED_LIBRARIES}) add_executable(fdmdv_put_test_bits fdmdv_put_test_bits.c fdmdv.c kiss_fft.c) target_link_libraries(fdmdv_put_test_bits ${CMAKE_REQUIRED_LIBRARIES}) +add_executable(fdmdv_channel fdmdv_channel.c fdmdv.c kiss_fft.c) +target_link_libraries(fdmdv_channel ${CMAKE_REQUIRED_LIBRARIES}) + add_executable(fdmdv_interleave fdmdv_interleave.c) target_link_libraries(fdmdv_interleave ${CMAKE_REQUIRED_LIBRARIES}) diff --git a/codec2-dev/src/fdmdv.c b/codec2-dev/src/fdmdv.c index d16b3184..8b7448ca 100644 --- a/codec2-dev/src/fdmdv.c +++ b/codec2-dev/src/fdmdv.c @@ -1750,6 +1750,37 @@ void fdmdv_dump_osc_mags(struct FDMDV *f) fprintf(stderr, "\n\n"); } + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: randn() + AUTHOR......: David Rowe + DATE CREATED: 2 August 2014 + + Simple approximation to normal (gaussian) random number generator + with 0 mean and unit variance. + +\*---------------------------------------------------------------------------*/ + +#define RANDN_IT 12 /* This magic number of iterations gives us a + unit variance. I think beacuse var = + (b-a)^2/12 for one uniform random variable, so + for a sum of n random variables it's + n(b-a)^2/12, or for b=1, a = 0, n=12, we get + var = 12(1-0)^2/12 = 1 */ + +static float randn() { + int i; + float rn = 0.0; + + for(i=0; isig_pwr_av = 0.9*f->sig_pwr_av + 0.1*(sig_pwr/nin); + sig_pwr /= nin; + + f->sig_pwr_av = 0.9*f->sig_pwr_av + 0.1*sig_pwr; /* det noise to meet target SNR */ target_snr_linear = powf(10.0, target_snr/10.0); - noise_pwr = f->sig_pwr_av/target_snr_linear; - noise_pwr_4k = 0.75*noise_pwr; /* this is the equivalent power scaled to a 4000 Hz BW */ - - /* for convenience we are using a uniform random number generator */ - - var_uniform = 1.0/12.0; - noise_gain = sqrtf(0.5*noise_pwr_4k/var_uniform); + noise_pwr = f->sig_pwr_av/target_snr_linear; /* noise pwr in a 3000 Hz BW */ + noise_pwr_1Hz = noise_pwr/3000.0; /* noise pwr in a 1 Hz bandwidth */ + noise_pwr_4000Hz = noise_pwr_1Hz*4000.0; /* noise pwr in a 4000 Hz BW, which + due to fs=8000 Hz in our simulation noise BW */ + noise_gain = sqrtf(0.5*noise_pwr_4000Hz); /* split noise pwr between real and imag sides */ + for(i=0; isig_pwr_av: %f target_snr_linear: %f noise_pwr_4k: %f noise_gain: %f\n", - // f->sig_pwr_av, target_snr_linear, noise_pwr_4k, noise_gain); + /* + fprintf(stderr, "sig_pwr: %f f->sig_pwr_av: %e target_snr_linear: %f noise_pwr_4000Hz: %e noise_gain: %e\n", + sig_pwr, f->sig_pwr_av, target_snr_linear, noise_pwr_4000Hz, noise_gain); + */ } diff --git a/codec2-dev/src/fdmdv_channel.c b/codec2-dev/src/fdmdv_channel.c new file mode 100644 index 00000000..e4d87cbe --- /dev/null +++ b/codec2-dev/src/fdmdv_channel.c @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: fdmdv_channel.c + AUTHOR......: David Rowe + DATE CREATED: 2 August 2014 + + Given an input raw file (8kHz, 16 bit shorts) of FDMDV modem + samples, adds channel impairments and outputs to another raw file. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2014 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, 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 "codec2_fdmdv.h" + +int main(int argc, char *argv[]) +{ + FILE *fin, *fout; + short rx_fdm_buf[FDMDV_NOM_SAMPLES_PER_FRAME]; + COMP rx_fdm[FDMDV_NOM_SAMPLES_PER_FRAME]; + struct FDMDV *fdmdv; + float snrdB, sam; + int i; + + if (argc < 3) { + printf("usage: %s InputModemRawFile OutputModemRawFile SNRdB\n", argv[0]); + printf("e.g %s test_in.raw test_out.raw 4\n", argv[0]); + exit(1); + } + + if (strcmp(argv[1], "-") == 0) fin = stdin; + else if ( (fin = fopen(argv[1],"rb")) == NULL ) { + fprintf(stderr, "Error opening input modem sample file: %s: %s.\n", + argv[1], strerror(errno)); + exit(1); + } + + if (strcmp(argv[2], "-") == 0) fout = stdout; + else if ( (fout = fopen(argv[2],"wb")) == NULL ) { + fprintf(stderr, "Error opening output modem sample file: %s: %s.\n", + argv[2], strerror(errno)); + exit(1); + } + + snrdB = atof(argv[3]); + fdmdv = fdmdv_create(FDMDV_NC); + + while(fread(rx_fdm_buf, sizeof(short), FDMDV_NOM_SAMPLES_PER_FRAME, fin) == FDMDV_NOM_SAMPLES_PER_FRAME) { + + for(i=0; i 32767.0) sam = 32767.0; + if (sam < -32767.0) sam = -32767.0; + rx_fdm_buf[i] = sam; + } + + fwrite(rx_fdm_buf, sizeof(short), FDMDV_NOM_SAMPLES_PER_FRAME, fout); + + /* if this is in a pipeline, we probably don't want the usual + buffering to occur */ + + if (fout == stdout) fflush(stdout); + if (fin == stdin) fflush(stdin); + } + + fclose(fin); + fclose(fout); + fdmdv_destroy(fdmdv); + + return 0; +} + diff --git a/codec2-dev/src/globals.c b/codec2-dev/src/globals.c deleted file mode 100644 index f2182f79..00000000 --- a/codec2-dev/src/globals.c +++ /dev/null @@ -1,49 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - FILE........: globals.c - AUTHOR......: David Rowe - DATE CREATED: 11/5/94 - - Globals for sinusoidal speech coder. - -\*---------------------------------------------------------------------------*/ - -/* - Copyright (C) 2009 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 "sine.h" /* global defines for coder */ - -/* Globals used in encoder and decoder */ - -int frames; /* number of frames processed so far */ -float Sn[M]; /* float input speech samples */ -MODEL model; /* model parameters for the current frame */ -int Nw; /* number of samples in analysis window */ -float sig; /* energy of current frame */ - -/* Globals used in encoder */ - -float w[M]; /* time domain hamming window */ -COMP W[FFT_ENC]; /* DFT of w[] */ -COMP Sw[FFT_ENC]; /* DFT of current frame */ - -/* Globals used in decoder */ - -COMP Sw_[FFT_ENC]; /* DFT of all voiced synthesised signal */ -float Sn_[AW_DEC]; /* synthesised speech */ -float Pn[AW_DEC]; /* time domain Parzen (trapezoidal) window */ - diff --git a/codec2-dev/src/globals.h b/codec2-dev/src/globals.h deleted file mode 100644 index cef72034..00000000 --- a/codec2-dev/src/globals.h +++ /dev/null @@ -1,47 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - FILE........: globals.h - AUTHOR......: David Rowe - DATE CREATED: 1/11/94 - - Globals for sinusoidal speech coder. - -\*---------------------------------------------------------------------------*/ - -/* - Copyright (C) 2009 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 . -*/ - -/* Globals used in encoder and decoder */ - -extern int frames; /* number of frames processed so far */ -extern float Sn[]; /* float input speech samples */ -extern MODEL model; /* model parameters for the current frame */ -extern int Nw; /* number of samples in analysis window */ -extern float sig; /* energy of current frame */ - -/* Globals used in encoder */ - -extern float w[]; /* time domain hamming window */ -extern COMP W[]; /* frequency domain hamming window */ -extern COMP Sw[]; /* DFT of current frame */ -extern COMP Sw_[]; /* DFT of all voiced synthesised signal */ - -/* Globals used in decoder */ - -extern float Sn_[]; /* output synthesised speech samples */ -extern float Pn[]; /* time domain Parzen (trapezoidal) window */ - diff --git a/codec2-dev/src/glottal.c b/codec2-dev/src/glottal.c deleted file mode 100644 index 8ac3ff4a..00000000 --- a/codec2-dev/src/glottal.c +++ /dev/null @@ -1,257 +0,0 @@ -const float glottal[]={ - 0.000000, - -0.057687, - -0.115338, - -0.172917, - -0.230385, - -0.287707, - -0.344845, - -0.401762, - -0.458419, - -0.514781, - -0.570809, - -0.626467, - -0.681721, - -0.736537, - -0.790884, - -0.844733, - -0.898057, - -0.950834, - -1.003044, - -1.054670, - -1.105700, - -1.156124, - -1.205936, - -1.255132, - -1.303711, - -1.351675, - -1.399026, - -1.445769, - -1.491908, - -1.537448, - -1.582393, - -1.626747, - -1.670514, - -1.713693, - -1.756285, - -1.798288, - -1.839697, - -1.880507, - -1.920712, - -1.960302, - -1.999269, - -2.037603, - -2.075295, - -2.112335, - -2.148716, - -2.184430, - -2.219472, - -2.253839, - -2.287531, - -2.320550, - -2.352900, - -2.384588, - -2.415623, - -2.446019, - -2.475788, - -2.504946, - -2.533512, - -2.561501, - -2.588934, - -2.615827, - -2.642198, - -2.668064, - -2.693439, - -2.718337, - -2.742767, - -2.766738, - -2.790256, - -2.813322, - -2.835936, - -2.858094, - -2.879790, - -2.901016, - -2.921759, - -2.942008, - -2.961747, - -2.980961, - -2.999632, - -3.017745, - -3.035282, - -3.052228, - -3.068567, - -3.084285, - -3.099371, - -3.113813, - -3.127605, - -3.140738, - 3.129975, - 3.118167, - 3.107022, - 3.096537, - 3.086709, - 3.077531, - 3.068996, - 3.061096, - 3.053821, - 3.047159, - 3.041102, - 3.035636, - 3.030753, - 3.026441, - 3.022690, - 3.019491, - 3.016836, - 3.014718, - 3.013132, - 3.012072, - 3.011535, - 3.011521, - 3.012028, - 3.013057, - 3.014612, - 3.016695, - 3.019310, - 3.022463, - 3.026160, - 3.030407, - 3.035212, - 3.040580, - 3.046520, - 3.053038, - 3.060141, - 3.067836, - 3.076128, - 3.085023, - 3.094525, - 3.104639, - 3.115367, - 3.126712, - 3.138674, - -3.131930, - -3.118731, - -3.104915, - -3.090485, - -3.075444, - -3.059795, - -3.043543, - -3.026695, - -3.009254, - -2.991229, - -2.972625, - -2.953449, - -2.933710, - -2.913414, - -2.892567, - -2.871176, - -2.849248, - -2.826787, - -2.803798, - -2.780284, - -2.756247, - -2.731689, - -2.706609, - -2.681005, - -2.654875, - -2.628213, - -2.601015, - -2.573272, - -2.544977, - -2.516121, - -2.486694, - -2.456686, - -2.426084, - -2.394879, - -2.363060, - -2.330616, - -2.297538, - -2.263816, - -2.229444, - -2.194416, - -2.158727, - -2.122375, - -2.085359, - -2.047682, - -2.009347, - -1.970361, - -1.930732, - -1.890470, - -1.849587, - -1.808098, - -1.766017, - -1.723360, - -1.680145, - -1.636388, - -1.592105, - -1.547313, - -1.502025, - -1.456256, - -1.410016, - -1.363314, - -1.316157, - -1.268547, - -1.220486, - -1.171971, - -1.122997, - -1.073555, - -1.023636, - -0.973227, - -0.922312, - -0.870875, - -0.818899, - -0.766366, - -0.713257, - -0.659554, - -0.605242, - -0.550303, - -0.494723, - -0.438492, - -0.381598, - -0.324036, - -0.265800, - -0.206889, - -0.147303, - -0.087046, - -0.026121, - 0.035463, - 0.097698, - 0.160576, - 0.224087, - 0.288221, - 0.352969, - 0.418323, - 0.484276, - 0.550822, - 0.617958, - 0.685681, - 0.753991, - 0.822889, - 0.892378, - 0.962462, - 1.033144, - 1.104430, - 1.176325, - 1.248833, - 1.321956, - 1.395696, - 1.470051, - 1.545019, - 1.620593, - 1.696763, - 1.773516, - 1.850837, - 1.928705, - 2.007097, - 2.085987, - 2.165347, - 2.245145, - 2.325347, - 2.405919, - 2.486824, - 2.568025, - 2.649485, - 2.731167, - 2.813033, - 2.895045, - 2.977167, - 3.059362};