From aaf86486819b4d5071f6b3e45daa1562e7ce4e5b Mon Sep 17 00:00:00 2001 From: baobrien Date: Fri, 12 Feb 2016 03:17:01 +0000 Subject: [PATCH] Testing of fmfsk is a bit more fleshed out; fsk has been made more suitable for DV bitrate operaton git-svn-id: https://svn.code.sf.net/p/freetel/code@2697 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fmfsk.m | 14 +-- codec2-dev/octave/tfmfsk.m | 57 +++++------ codec2-dev/src/fmfsk.c | 2 +- codec2-dev/src/fsk.c | 153 +++++++++++++++++++++++++---- codec2-dev/src/fsk.h | 17 +++- codec2-dev/src/fsk_demod.c | 54 ++++++---- codec2-dev/src/fsk_get_test_bits.c | 39 ++++---- 7 files changed, 228 insertions(+), 108 deletions(-) diff --git a/codec2-dev/octave/fmfsk.m b/codec2-dev/octave/fmfsk.m index 54627b69..1d3cfabb 100644 --- a/codec2-dev/octave/fmfsk.m +++ b/codec2-dev/octave/fmfsk.m @@ -30,17 +30,13 @@ function states = fmfsk_init(Fs,Rb) assert(mod(Fs,Rb*2)==0); %Current fixed processing buffer size, in non-ME bits - nbit = 192; - + nbit = 96; states.Rb = Rb; states.Rs = Rb*2; % Manchester-encoded bitrate states.Fs = Fs; states.Ts = Fs/states.Rs; states.N = nbit*2*states.Ts; - %states.N - %states.N = floor(states.Rs*.080)*states.Ts - %states.N states.nin = states.N; % Samples in the next demod cycle states.nstash = states.Ts*2; % How many samples to stash away between proc cycles for timing adjust states.nmem = states.N+(4*states.Ts); @@ -143,7 +139,7 @@ function [rx_bits states] = fmfsk_demod(states,rx) apeven = 0; apodd = 0; - sample_offset = (Ts/2)+Ts+rx_timing-1 + sample_offset = (Ts/2)+Ts+rx_timing-1; symsamp = zeros(1,nsym); @@ -251,8 +247,6 @@ function fmfsk_run_sim(EbNodB,timing_offset=0,de=0,of=0,hpf=0) [b, a] = cheby1(4, 1, 300/Fs, 'high'); % 300Hz HPF to simulate FM radios tx_pmod = fmfsk_mod(states, tx_bits); - figure(10) - plot(tx_pmod); tx = analog_fm_mod(fm_states, tx_pmod); tx = tx(10:length(tx)); @@ -276,8 +270,10 @@ function fmfsk_run_sim(EbNodB,timing_offset=0,de=0,of=0,hpf=0) rx = filter(filt,1,rx); figure(4) + title("Spectrum of rx-ed signal after FM demod and FM radio channel"); plot(20*log10(abs(fft(rx)))) - figure(5) + figure(5) + title("Time domain of rx-ed signal after FM demod and FM radio channel"); plot(rx) %rx = real(rx); %b1 = fir2(100, [0 4000 5200 48000]/48000, [1 1 0.5 0.5]); diff --git a/codec2-dev/octave/tfmfsk.m b/codec2-dev/octave/tfmfsk.m index a60ba1e8..8af56c07 100644 --- a/codec2-dev/octave/tfmfsk.m +++ b/codec2-dev/octave/tfmfsk.m @@ -103,19 +103,13 @@ endfunction function test_stats = fmfsk_demod_xt(Fs,Rs,mod,tname,M=2) global tfsk_location; - size(mod) - %Name of executable containing the modulator modvecfilename = sprintf('fmfsk_demod_ut_modvec_%d',getpid()); bitvecfilename = sprintf('fmfsk_demod_ut_bitvec_%d',getpid()); tvecfilename = sprintf('fmfsk_demod_ut_tracevec_%d.txt',getpid()); - delete(bitvecfilename); - delete(modvecfilename); - delete(tvecfilename); - %command to be run by system to launch the demod - command = sprintf('%s D %d %d %s %s %s',tfsk_location,Fs,Rs,modvecfilename,bitvecfilename,tvecfilename) + command = sprintf('%s D %d %d %s %s %s',tfsk_location,Fs,Rs,modvecfilename,bitvecfilename,tvecfilename); %save modulated input into a file modvecfile = fopen(modvecfilename,'wb+'); @@ -303,7 +297,7 @@ function stats = tfmfsk_run_sim(EbNodB,timing_offset=0,de=0,of=0,hpf=0,df=0,M=2) end %Add frequency drift - fdrift = df/Fs + fdrift = df/Fs; fshift = 2*pi*fdrift*(1:length(tx)); fshift = exp(j*(fshift.^2)); tx = tx.*fshift; @@ -315,7 +309,6 @@ function stats = tfmfsk_run_sim(EbNodB,timing_offset=0,de=0,of=0,hpf=0,df=0,M=2) %High-pass filter to simulate the FM radios if hpf>0 - printf("high-pass filtering!\n") rx = filter(b,a,rx); end @@ -390,13 +383,13 @@ function pass = ebno_battery_test(timing_offset,drift,hpf,deemp,outfilt) ebnodbs = length(ebnodbrange); %Replication of other parameters for parcellfun - timingv = repmat(timing_offset,1,ebnodbs); - driftv = repmat(timing_offset,1,ebnodbs); - hpfv = repmat(timing_offset,1,ebnodbs); - deempv = repmat(timing_offset,1,ebnodbs); - outfv = repmat(timing_offset,1,ebnodbs); + timingv = repmat(timing_offset ,1,ebnodbs); + driftv = repmat(drift ,1,ebnodbs); + hpfv = repmat(hpf ,1,ebnodbs); + deempv = repmat(deemp ,1,ebnodbs); + outfv = repmat(outfilt ,1,ebnodbs); - statv = pararrayfun(floor(1.25*nproc()),@tfmfsk_run_sim,ebnodbrange,timingv,deempv,outfv,hpfv,driftv); + statv = pararrayfun(floor(.75*nproc()),@tfmfsk_run_sim,ebnodbrange,timingv,deempv,outfv,hpfv,driftv); %statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,mv); passv = zeros(1,length(statv)); @@ -429,7 +422,7 @@ function pass = test_drift_var(hpf,deemp,outfilt) endfunction function pass = test_fmfsk_battery() - pass = test_mod_horuscfg_randbits; + pass = test_mod_fdvbcfg_randbits; assert(pass) pass = pass && test_drift_var(1,1,1); assert(pass) @@ -439,24 +432,18 @@ function pass = test_fmfsk_battery() endfunction function plot_fsk_bers(M=2) - %Range of EbNodB over which to plot - ebnodbrange = (3:13); - - berc = ones(1,length(ebnodbrange)); - bero = ones(1,length(ebnodbrange)); - berinc = ones(1,length(ebnodbrange)); - ebnodbs = length(ebnodbrange) - mode = 2; - %Replication of other parameters for parcellfun - modev = repmat(mode,1,ebnodbs); - timingv = repmat(0,1,ebnodbs); - fadingv = repmat(0,1,ebnodbs); - dfv = repmat(0,1,ebnodbs); - dav = repmat(1,1,ebnodbs); - Mv = repmat(M,1,ebnodbs); + %Range of EbNodB over which to test + ebnodbrange = (8:14); + ebnodbs = length(ebnodbrange); + %Replication of other parameters for parcellfun + timingv = repmat(1 ,1,ebnodbs); + driftv = repmat(1 ,1,ebnodbs); + hpfv = repmat(1 ,1,ebnodbs); + deempv = repmat(1 ,1,ebnodbs); + outfv = repmat(1 ,1,ebnodbs); - statv = pararrayfun(floor(1.25*nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv); + statv = pararrayfun(nproc(),@tfmfsk_run_sim,ebnodbrange,timingv,deempv,outfv,hpfv,driftv); %statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv); for ii = (1:length(statv)) @@ -468,10 +455,10 @@ function plot_fsk_bers(M=2) clf; figure(M) - semilogy(ebnodbrange, berinc,sprintf('r;%dFSK non-coherent theory;',M)) + semilogy(ebnodbrange, berinc,sprintf('r;2FSK non-coherent theory;',M)) hold on; - semilogy(ebnodbrange, bero ,sprintf('g;Octave fsk horus %dFSK Demod;',M)) - semilogy(ebnodbrange, berc,sprintf('v;C fsk horus %dFSK Demod;',M)) + semilogy(ebnodbrange, bero ,sprintf('g;Octave ME-FM-FSK Demod;',M)) + semilogy(ebnodbrange, berc,sprintf('v;C ME-FM-FSK Demod;',M)) hold off; grid("minor"); axis([min(ebnodbrange) max(ebnodbrange) 1E-5 1]) diff --git a/codec2-dev/src/fmfsk.c b/codec2-dev/src/fmfsk.c index 3b04163d..06878557 100644 --- a/codec2-dev/src/fmfsk.c +++ b/codec2-dev/src/fmfsk.c @@ -38,7 +38,7 @@ #include "modem_probe.h" #include "comp_prim.h" -#define STD_PROC_BITS 192 +#define STD_PROC_BITS 96 /* diff --git a/codec2-dev/src/fsk.c b/codec2-dev/src/fsk.c index c887cf8b..8333f8ed 100644 --- a/codec2-dev/src/fsk.c +++ b/codec2-dev/src/fsk.c @@ -32,7 +32,7 @@ \*---------------------------------------------------------------------------*/ /* P oversampling rate constant -- should probably be init-time configurable */ -#define ct_P 8 +#define horus_P 8 /* Define this to enable EbNodB estimate */ /* This needs square roots, may take more cpu time than it's worth */ @@ -110,6 +110,129 @@ static inline COMP comp_normalize(COMP a){ } +/*---------------------------------------------------------------------------*\ + + FUNCTION....: fsk_create_hbr + AUTHOR......: Brady O'Brien + DATE CREATED: 11 February 2016 + + Create and initialize an instance of the FSK modem. Returns a pointer + to the modem state/config struct. One modem config struct may be used + for both mod and demod. returns NULL on failure. + +\*---------------------------------------------------------------------------*/ + +struct FSK * fsk_create_hbr(int Fs, int Rs,int P,int M, int tx_f1, int tx_fs) +{ + struct FSK *fsk; + int i; + int memold; + int Ndft = 0; + + /* Number of symbols in a processing frame */ + int nsyms = 192; + + /* Check configuration validity */ + assert(Fs > 0 ); + assert(Rs > 0 ); + assert(tx_f1 > 0); + assert(tx_fs > 0); + assert(P > 0); + /* Ts (Fs/Rs) must be an integer */ + assert( (Fs%Rs) == 0 ); + /* Ts/P (Fs/Rs/P) must be an integer */ + assert( ((Fs/Rs)%P) == 0 ); + assert( M==2 || M==4); + + fsk = (struct FSK*) malloc(sizeof(struct FSK)); + if(fsk == NULL) return NULL; + + + /* Set constant config parameters */ + fsk->Fs = Fs; + fsk->Rs = Rs; + fsk->Ts = Fs/Rs; + fsk->N = fsk->Ts*nsyms; + fsk->P = P; + fsk->Nsym = nsyms; + fsk->Nmem = fsk->N+(2*fsk->Ts); + fsk->f1_tx = tx_f1; + fsk->fs_tx = tx_fs; + fsk->nin = fsk->N; + fsk->mode = M==2 ? MODE_2FSK : MODE_4FSK; + fsk->Nbits = M==2 ? fsk->Nsym : fsk->Nsym*2; + + /* Find smallest 2^N value that fits Fs for efficient FFT */ + /* It would probably be better to use KISS-FFt's routine here */ + for(i=1; i; i<<=1) + if((fsk->N)&i) + Ndft = i; + + fsk->Ndft = Ndft; + + fsk->est_min = tx_f1-3*Rs; + if(fsk->est_min<0) fsk->est_min = 0; + + fsk->est_max = (Fs/2)-Rs; + + fsk->est_space = (Rs/2)-(Rs/10); + + /* Set up rx state */ + fsk->phi1_c = comp_exp_j(0); + fsk->phi2_c = comp_exp_j(0); + fsk->phi3_c = comp_exp_j(0); + fsk->phi4_c = comp_exp_j(0); + + fsk->phi1_c.real = 0; + fsk->phi1_c.imag = 1; + + memold = (4*fsk->Ts); + + fsk->nstash = memold; + fsk->samp_old = (float*) malloc(sizeof(float)*memold); + if(fsk->samp_old == NULL){ + free(fsk); + return NULL; + } + + for(i=0;isamp_old[i]=0; + + fsk->fft_cfg = kiss_fftr_alloc(fsk->Ndft,0,NULL,NULL); + if(fsk->fft_cfg == NULL){ + free(fsk->samp_old); + free(fsk); + return NULL; + } + + fsk->fft_est = (float*)malloc(sizeof(float)*fsk->Ndft/2); + if(fsk->fft_est == NULL){ + free(fsk->samp_old); + free(fsk->fft_cfg); + free(fsk); + return NULL; + } + + for(i=0;iNdft/2;i++)fsk->fft_est[i] = 0; + + fsk->norm_rx_timing = 0; + + /* Set up tx state */ + fsk->tx_phase_c = comp_exp_j(0); + + /* Set up demod stats */ + fsk->EbNodB = 0; + fsk->f1_est = 0; + fsk->f2_est = 0; + fsk->f3_est = 0; + fsk->f4_est = 0; + fsk->ppm = 0; + + return fsk; +} + +#define HORUS_MIN 800 +#define HORUS_MAX 2500 +#define HORUS_MIN_SPACING 200 /*---------------------------------------------------------------------------*\ @@ -135,31 +258,24 @@ struct FSK * fsk_create(int Fs, int Rs,int M, int tx_f1, int tx_fs) assert(Rs > 0 ); assert(tx_f1 > 0); assert(tx_fs > 0); - assert(ct_P > 0); + assert(horus_P > 0); /* Ts (Fs/Rs) must be an integer */ assert( (Fs%Rs) == 0 ); /* Ts/P (Fs/Rs/P) must be an integer */ - assert( ((Fs/Rs)%ct_P) == 0 ); + assert( ((Fs/Rs)%horus_P) == 0 ); assert( M==2 || M==4); fsk = (struct FSK*) malloc(sizeof(struct FSK)); if(fsk == NULL) return NULL; - - /* Find smallest 2^N value that fits Fs for efficient FFT */ - /* It would probably be better to use KISS-FFt's routine here */ - for(i=1; i; i<<=1) - if(Fs&i) - Ndft = i<<1; Ndft = 1024; - //Ndft = 4096; /* Set constant config parameters */ fsk->Fs = Fs; fsk->Rs = Rs; fsk->Ts = Fs/Rs; fsk->N = Fs; - fsk->P = ct_P; + fsk->P = horus_P; fsk->Nsym = fsk->N/fsk->Ts; fsk->Ndft = Ndft; fsk->Nmem = fsk->N+(2*fsk->Ts); @@ -168,6 +284,9 @@ struct FSK * fsk_create(int Fs, int Rs,int M, int tx_f1, int tx_fs) fsk->nin = fsk->N; fsk->mode = M==2 ? MODE_2FSK : MODE_4FSK; fsk->Nbits = M==2 ? fsk->Nsym : fsk->Nsym*2; + fsk->est_min = HORUS_MIN; + fsk->est_max = HORUS_MAX; + fsk->est_space = HORUS_MIN_SPACING; /* Set up rx state */ fsk->phi1_c = comp_exp_j(0); @@ -232,9 +351,6 @@ void fsk_destroy(struct FSK *fsk){ free(fsk); } -#define FEST_MIN 800 -#define FEST_MAX 2500 -#define FEST_MIN_SPACING 200 /* * Internal function to estimate the frequencies of the two tones within a block of samples. @@ -265,9 +381,9 @@ void fsk_demod_freq_est(struct FSK *fsk, float fsk_in[],float *freqs,int M){ kiss_fft_cpx *fftout = (kiss_fft_cpx*)alloca(sizeof(kiss_fft_cpx)*(Ndft/2)+1); fft_samps = Ndft; - f_min = (FEST_MIN*Ndft)/Fs; - f_max = (FEST_MAX*Ndft)/Fs; - f_zero = (FEST_MIN_SPACING*Ndft)/(2*Fs); + f_min = (fsk->est_min*Ndft)/Fs; + f_max = (fsk->est_max*Ndft)/Fs; + f_zero = (fsk->est_space*Ndft)/(2*Fs); int fft_loops = nin/Ndft; for(j=0; j6) - modem_probe_init("fsk2",argv[6]); + if(argc>7) + modem_probe_init("fsk2",argv[7]); /* set up FSK */ - fsk = fsk_create(Fs,Rs,M,1200,400); + if(!hbr) + fsk = fsk_create(Fs,Rs,M,1200,400); + else + fsk = fsk_create_hbr(Fs,Rs,P,M,1200,400); + if(fin==NULL || fout==NULL || fsk==NULL){ fprintf(stderr,"Couldn't open test vector files\n"); @@ -102,6 +117,7 @@ int main(int argc,char *argv[]){ } } + modem_probe_close(); cleanup: fclose(fin); diff --git a/codec2-dev/src/fsk_get_test_bits.c b/codec2-dev/src/fsk_get_test_bits.c index a303e8fb..d4f38bad 100644 --- a/codec2-dev/src/fsk_get_test_bits.c +++ b/codec2-dev/src/fsk_get_test_bits.c @@ -43,8 +43,6 @@ int main(int argc,char *argv[]){ FILE *fout; uint8_t *bitbuf; - /* Seed the RNG with some constant */ - srand(158324); if(argc != 3){ fprintf(stderr,"usage: %s OutputBitsOnePerByte FrameCount\n",argv[0]); @@ -69,29 +67,26 @@ int main(int argc,char *argv[]){ bitbuf = (uint8_t*)alloca(sizeof(uint8_t)*FSK_FRAME_SIZE); /* Write out sync frame and sequence */ - for(i=0; i