diff options
Diffstat (limited to 'src/ch.c')
| -rw-r--r-- | src/ch.c | 479 |
1 files changed, 479 insertions, 0 deletions
diff --git a/src/ch.c b/src/ch.c new file mode 100644 index 0000000..57411bd --- /dev/null +++ b/src/ch.c @@ -0,0 +1,479 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: ch.c + AUTHOR......: David Rowe + DATE CREATED: May 2015 + + Channel simulation program for testing command line versions of modems. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2015-2022 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 <string.h> +#include <math.h> +#include <errno.h> +#include <getopt.h> + +#include "freedv_api.h" +#include "codec2_cohpsk.h" +#include "comp_prim.h" +#include "ht_coeff.h" +#include "ssbfilt_coeff.h" + +#include "debug_alloc.h" + +#define BUF_N 160 +#define MPG_DELAY_MS 0.5 +#define MPP_DELAY_MS 2.0 +#define MPD_DELAY_MS 4.0 + +/* see instructions below for how to generate these files */ + +#define DEFAULT_FADING_DIR "unittest" +#define MPG_FADING_FILE_NAME "slow_fading_samples.float" +#define MPP_FADING_FILE_NAME "fast_fading_samples.float" +#define MPD_FADING_FILE_NAME "faster_fading_samples.float" + +// Gaussian from uniform: +float gaussian(void) { + double x = (double)rand() / RAND_MAX; + double y = (double)rand() / RAND_MAX; + double z = sqrt(-2 * log(x)) * cos(2 * M_PI * y); + return sqrt(1./2.) * z; +} + +// complex noise sample +COMP noise(void) { + COMP n = {gaussian(),gaussian()}; + return n; +} + +int main(int argc, char *argv[]) +{ + FILE *fin, *ffading, *fout; + char *fading_dir; + float NodB, foff_hz; + int fading_en, nhfdelay; + + short buf[BUF_N]; + float htbuf[HT_N+BUF_N]; + COMP ch_in[BUF_N]; + COMP ch_fdm[BUF_N]; + COMP ssbfiltbuf[SSBFILT_N+BUF_N]; + COMP ssbfiltout[BUF_N]; + + COMP phase_ch; + float No, variance; + COMP scaled_noise; + float hf_gain; + COMP *ch_fdm_delay = NULL, aspread, aspread_2ms, delayed, direct; + float tx_pwr, tx_pwr_fade, noise_pwr, user_multipath_delay; + int frames, i, j, k, Fs, ret, nclipped, noutclipped, ssbfilt_en, complex_out, ctest; + float sam, peak, clip, papr, CNo, snr3k, gain; + + if (argc < 3) { + helpmsg: + fprintf(stderr, "Command line channel simulation tool.\n" + "\n" + "usage: %s InputRealModemRawFile OutputRealModemRawFile [Options]\n" + "\n" + " real int16 input -> Gain -> Hilbert Transform -> clipper -> freq shift ->\n" + " Multipath -> AWGN noise -> SSB filter -> real int16 output\n" + "\n" + "[--clip int16] Hilbert clipper (clip complex signal magnitude, default 32767)\n" + "[--complexout] Optional int16 IQ complex output (default real int16)\n" + "[--ctest] Check PAPR is around 0dB, used to support ctests\n" + "[--freq FoffHz] Frequency offset (default 0Hz)\n" + "[--fading_dir Path] path to multipath fading files (default 'unittest')\n" + "[--Fs SampleRateHz] Sample rate of simulation (default 8000 Hz)\n" + "[--gain G] Linear gain (default 1.0)\n" + "[--mpg] Multipath good 0.1Hz Doppler, 0.5ms delay\n" + "[--mpp] Multipath poor 1.0Hz Doppler, 2.0ms delay\n" + "[--mpd] Multipath disturbed 2.0Hz Doppler, 4.0ms delay\n" + "[--ssbfilt 0|1] SSB bandwidth filter (default 1 on)\n" + "[--mulipath_delay ms] Optionally adjust multipath delay\n" + "[--No dBHz] AWGN Noise density dB/Hz (default -100)" + "\n" + , argv[0]); + exit(1); + } + + if (strcmp(argv[1], "-") == 0) fin = stdin; + else if ( (fin = fopen(argv[1],"rb")) == NULL ) { + fprintf(stderr, "ch: Error opening input modem raw 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, "ch: Error opening output modem raw file: %s: %s.\n", + argv[2], strerror(errno)); + exit(1); + } + + NodB = -100; + Fs = 8000; foff_hz = 0.0; fading_en = 0; ctest = 0; + clip =32767; gain = 1.0; + ssbfilt_en = 1; complex_out = 0; + fading_dir = strdup(DEFAULT_FADING_DIR); user_multipath_delay = -1.0; + + int o = 0; + int opt_idx = 0; + while( o != -1 ){ + static struct option long_opts[] = { + {"complexout", no_argument, 0, 'o'}, + {"ctest", no_argument, 0, 't'}, + {"clip", required_argument, 0, 'c'}, + {"fading_dir", required_argument, 0, 'u'}, + {"freq", required_argument, 0, 'f'}, + {"Fs", required_argument, 0, 'r'}, + {"gain", required_argument, 0, 'g'}, + {"ssbfilt", required_argument, 0, 's'}, + {"help", no_argument, 0, 'h'}, + {"mpg", no_argument, 0, 'i'}, + {"mpp", no_argument, 0, 'p'}, + {"mpd", no_argument, 0, 'd'}, + {"multipath_delay", required_argument, 0, 'm'}, + {"No", required_argument, 0, 'n'}, + {0, 0, 0, 0} + }; + + o = getopt_long(argc,argv,"c:df:g:im:n:opr:s:tu:h",long_opts,&opt_idx); + + switch(o) { + case 'c': + clip = atof(optarg); + break; + case 'd': + fading_en = 3; + break; + case 'f': + foff_hz = atof(optarg); + break; + case 'g': + gain = atof(optarg); + break; + case 'i': + fading_en = 1; + break; + case 'm': + user_multipath_delay = atof(optarg); + break; + case 'n': + NodB = atof(optarg); + break; + case 'o': + complex_out = 1; + break; + case 'p': + fading_en = 2; + break; + case 'r': + Fs = atoi(optarg); + break; + case 's': + ssbfilt_en = atoi(optarg); + break; + case 't': + ctest = 1; + break; + case 'u': + fading_dir = strdup(optarg); + break; + case 'h': + case '?': + goto helpmsg; + break; + } + } + + phase_ch.real = 1.0; phase_ch.imag = 0.0; + + /* N = var = NoFs */ + + // arbitrary noise scaling, to maintain backwards compatibility with many tests. TODO make the No + // units more sensible, and fix all the tests that depend on this scaling + No = pow(10.0, NodB/10.0)*1000*1000; + variance = Fs*No; + + tx_pwr = tx_pwr_fade = noise_pwr = 0.0; + noutclipped = 0; nclipped = 0; + peak = 0.0; + + /* init HF fading model */ + + ffading = NULL; + nhfdelay = 0; + if (fading_en) { + char fname[256]; + + if (fading_en == 1) { + sprintf(fname, "%s/%s", fading_dir, MPG_FADING_FILE_NAME); + ffading = fopen(fname, "rb"); + if (ffading == NULL) { + cant_load_fading_file: + fprintf(stderr, "-----------------------------------------------------\n"); + fprintf(stderr, "ch ERROR: Can't find fading file: %s\n", fname); + fprintf(stderr, "\nAdjust path --fading_dir or use GNU Octave to generate:\n\n"); + gen_fading_file: + fprintf(stderr, "$ octave --no-gui\n"); + fprintf(stderr, "octave:24> pkg load signal\n"); + fprintf(stderr, "octave:24> time_secs=60\n"); + fprintf(stderr, "octave:25> ch_fading(\"faster_fading_samples.float\", 8000, 2.0, 8000*time_secs)\n"); + fprintf(stderr, "octave:26> ch_fading(\"fast_fading_samples.float\", 8000, 1.0, 8000*time_secs)\n"); + fprintf(stderr, "octave:27> ch_fading(\"slow_fading_samples.float\", 8000, 0.1, 8000*time_secs)\n"); + fprintf(stderr, "-----------------------------------------------------\n"); + exit(1); + } + nhfdelay = floor(MPG_DELAY_MS*Fs/1000); + } + + if (fading_en == 2) { + sprintf(fname, "%s/%s", fading_dir, MPP_FADING_FILE_NAME); + ffading = fopen(fname, "rb"); + if (ffading == NULL) goto cant_load_fading_file; + nhfdelay = floor(MPP_DELAY_MS*Fs/1000); + } + + if (fading_en == 3) { + sprintf(fname, "%s/%s", fading_dir, MPD_FADING_FILE_NAME); + ffading = fopen(fname, "rb"); + if (ffading == NULL) goto cant_load_fading_file; + nhfdelay = floor(MPD_DELAY_MS*Fs/1000); + } + + ch_fdm_delay = (COMP*)MALLOC((nhfdelay+COHPSK_NOM_SAMPLES_PER_FRAME)*sizeof(COMP)); + assert(ch_fdm_delay != NULL); + for(i=0; i<nhfdelay+COHPSK_NOM_SAMPLES_PER_FRAME; i++) { + ch_fdm_delay[i].real = 0.0; + ch_fdm_delay[i].imag = 0.0; + } + + /* optionally override delay from command line */ + if (user_multipath_delay >= 0.0) nhfdelay = floor(user_multipath_delay*Fs/1000); + + /* first values in file are HF gains */ + + for (i=0; i<4; i++) + ret = fread(&hf_gain, sizeof(float), 1, ffading); + //fprintf(stderr, "hf_gain: %f\n", hf_gain); + } + + assert(HT_N == sizeof(ht_coeff)/sizeof(COMP)); + for(i=0; i<HT_N; i++) { + htbuf[i] = 0.0; + } + for(i=0; i<SSBFILT_N; i++) { + ssbfiltbuf[i].real = 0.0; ssbfiltbuf[i].imag = 0.0; + } + COMP lo_phase = {1.0,0.0}; + COMP lo_freq; + lo_freq.real = cos(2.0*M_PI*SSBFILT_CENTRE/Fs); + lo_freq.imag = sin(2.0*M_PI*SSBFILT_CENTRE/Fs); + + fprintf(stderr, "ch: Fs: %d NodB: %4.2f foff: %4.2f Hz fading: %d nhfdelay: %d clip: %4.2f ssbfilt: %d complexout: %d\n", + Fs, NodB, foff_hz, fading_en, nhfdelay, clip, ssbfilt_en, complex_out); + + /* --------------------------------------------------------*\ + Main Loop + \*---------------------------------------------------------*/ + + frames = 0; + while(fread(buf, sizeof(short), BUF_N, fin) == BUF_N) { + frames++; + + /* Hilbert Transform to produce complex signal so we can do + single sided freq shifts, HF channel models, and analog + compression. Allows us to use real signal I/O. + + As the real and imag filters both have unity gain, ch_in[] + has twice the power of the real input signal buf[]. + */ + + for(i=0, j=HT_N; i<BUF_N; i++,j++) { + + htbuf[j] = (float)buf[i]*gain; + + /* FIR filter with HT to get imag, just delay to get real */ + + ch_in[i].real = 0.0; + ch_in[i].imag = 0.0; + for(k=0; k<HT_N; k++) { + ch_in[i].real += htbuf[j-k]*ht_coeff[k].real; + ch_in[i].imag += htbuf[j-k]*ht_coeff[k].imag; + } + //printf("%d %f %f\n", i, ch_in[i].real, ch_in[i].imag); + } + assert(j <= (BUF_N+HT_N)); + + /* update HT memory */ + for(i=0; i<HT_N; i++) + htbuf[i] = htbuf[i+BUF_N]; + + /* --------------------------------------------------------*\ + Clipping mag of complex signal + \*---------------------------------------------------------*/ + + for(i=0; i<BUF_N; i++) { + float mag = sqrt(ch_in[i].real*ch_in[i].real + ch_in[i].imag*ch_in[i].imag); + float angle = atan2(ch_in[i].imag, ch_in[i].real); + if (mag > clip) { + mag = clip; + nclipped++; + } + tx_pwr += mag*mag; + /* we get a bit of overshoot in peak measurements if HT filter hasn't been primed */ + if (frames*BUF_N > HT_N) + if (mag > peak) { + peak = mag; + //fprintf(stderr, "%d %f\n",frames, mag); + } + ch_in[i].real = mag*cos(angle); + ch_in[i].imag = mag*sin(angle); + } + + /* --------------------------------------------------------*\ + Channel + \*---------------------------------------------------------*/ + + fdmdv_freq_shift_coh(ch_fdm, ch_in, foff_hz, Fs, &phase_ch, BUF_N); + + /* optional HF fading -------------------------------------*/ + + if (fading_en) { + + /* update delayed signal buffer */ + + for(i=0; i<nhfdelay; i++) + ch_fdm_delay[i] = ch_fdm_delay[i+BUF_N]; + for(j=0; j<BUF_N; i++, j++) + ch_fdm_delay[i] = ch_fdm[j]; + + /* combine direct and delayed paths, both multiplied by + "spreading" (Doppler) functions */ + + for(i=0; i<BUF_N; i++) { + ret = fread(&aspread, sizeof(COMP), 1, ffading); + if (ret == 0) { + fprintf(stderr, "ch: Fading file finished - simulation stopping. You may need more samples:\n"); + goto gen_fading_file; + } + ret = fread(&aspread_2ms, sizeof(COMP), 1, ffading); + if (ret == 0) { + fprintf(stderr, "ch: Fading file finished - simulation stopping. You may need more samples:\n"); + goto gen_fading_file; + } + //printf("%f %f %f %f\n", aspread.real, aspread.imag, aspread_2ms.real, aspread_2ms.imag); + + direct = cmult(aspread, ch_fdm[i]); + delayed = cmult(aspread_2ms, ch_fdm_delay[i]); + ch_fdm[i] = fcmult(hf_gain, cadd(direct, delayed)); + } + } + + /* Measure power after fading model to make sure average pwr + is the same as AWGN channels. We only output the real + signal, which is half the power. */ + + for(i=0; i<BUF_N; i++) { + tx_pwr_fade += pow(ch_fdm[i].real, 2.0); + } + + /* AWGN noise ------------------------------------------*/ + + for(i=0; i<BUF_N; i++) { + COMP n = noise(); + scaled_noise = fcmult(sqrt(variance), n); + ch_fdm[i] = cadd(ch_fdm[i], scaled_noise); + noise_pwr += pow(scaled_noise.real, 2.0) + pow(scaled_noise.imag, 2.0); + } + + /* FIR filter to simulate (a rather flat) SSB filter. We + filter the complex signal by shifting it down to DC and + using real coefficients. */ + + for(i=0, j=SSBFILT_N; i<BUF_N; i++,j++) { + if (ssbfilt_en) { + ssbfiltbuf[j] = cmult(ch_fdm[i], cconj(lo_phase)); + ssbfiltout[i].real = 0.0; ssbfiltout[i].imag = 0.0; + for(k=0; k<SSBFILT_N; k++) { + ssbfiltout[i].real += ssbfiltbuf[j-k].real*ssbfilt_coeff[k]; + ssbfiltout[i].imag += ssbfiltbuf[j-k].imag*ssbfilt_coeff[k]; + } + ssbfiltout[i] = cmult(ssbfiltout[i], lo_phase); + lo_phase = cmult(lo_phase, lo_freq); + } + else { + ssbfiltout[i] = ch_fdm[i]; + } + } + + /* update SSB filter memory */ + for(i=0; i<SSBFILT_N; i++) + ssbfiltbuf[i] = ssbfiltbuf[i+BUF_N]; + + int nout = (complex_out+1)*BUF_N; + short bufout[nout], *pout=bufout; + for(i=0; i<BUF_N; i++) { + sam = ssbfiltout[i].real; + if (sam > 32767.0) { noutclipped++; sam = 32767.0; } + if (sam < -32767.0) { noutclipped++; sam = -32767.0; } + *pout++ = sam; + if (complex_out) { + sam = ssbfiltout[i].imag; + if (sam > 32767.0) { noutclipped++; sam = 32767.0; } + if (sam < -32767.0) { noutclipped++; sam = -32767.0; } + *pout++ = sam; + } + } + + fwrite(bufout, sizeof(short), nout, fout); + + /* if this is in a pipeline, we probably don't want the usual + buffering to occur */ + + if (fout == stdout) fflush(stdout); + } + + fclose(fin); + fclose(fout); + + int nsamples = frames*BUF_N; + papr = 10*log10(peak*peak/(tx_pwr/nsamples)); + CNo = 10*log10(tx_pwr/(noise_pwr/(Fs))); + snr3k = CNo - 10*log10(3000); + float outclipped_percent = noutclipped*100.0/nsamples; + fprintf(stderr, "ch: SNR3k(dB): %8.2f C/No....: %8.2f\n", snr3k, CNo); + fprintf(stderr, "ch: peak.....: %8.2f RMS.....: %8.2f CPAPR.....: %5.2f \n", peak, sqrt(tx_pwr/nsamples), papr); + fprintf(stderr, "ch: Nsamples.: %8d clipped.: %8.2f%% OutClipped: %5.2f%%\n", + nsamples, nclipped*100.0/nsamples, outclipped_percent); + if (outclipped_percent > 0.1) fprintf(stderr, "ch: WARNING output clipping\n"); + + if (ffading != NULL) fclose(ffading); + if (ch_fdm_delay != NULL) FREE(ch_fdm_delay); + if (ctest) { + /* special ctest mode: check CPAPR is around 0dB */ + if (fabs(papr) < 0.7) return 0; else return 1; + } + else return 0; +} |