From: baobrien Date: Tue, 9 Feb 2016 06:31:27 +0000 (+0000) Subject: fmfsk.c now producing correct bits. found and fixed stupid error. octave unit test... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=5592851cb808bb0c0c49795ac9568a9b5c0ab975;p=freetel-svn-tracking.git fmfsk.c now producing correct bits. found and fixed stupid error. octave unit test driver started. git-svn-id: https://svn.code.sf.net/p/freetel/code@2692 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/fmfsk.m b/codec2-dev/octave/fmfsk.m index 7697fd8e..66fd8ab9 100644 --- a/codec2-dev/octave/fmfsk.m +++ b/codec2-dev/octave/fmfsk.m @@ -113,9 +113,8 @@ function [rx_bits states] = fmfsk_demod(states,rx) en = st+Ts-1; rx_filt(ii) = sum(ssamps(st:en)); end - length(rx_filt) - + states.rx_filt = rx_filt; % Fine timing estimation ------------------------------------------------------ % Estimate fine timing using line at Rs/2 that Manchester encoding provides @@ -125,6 +124,7 @@ function [rx_bits states] = fmfsk_demod(states,rx) x = (rx_filt .^ 2) * exp(-j*w*(0:Np-1))'; norm_rx_timing = angle(x)/(2*pi) - 0.42; rx_timing = round(norm_rx_timing*Ts); + exp(-j*w*(0:Np-1))(1) %If rx timing is too far out, ask for more or less sample the next time % around to even it all out @@ -137,16 +137,19 @@ function [rx_bits states] = fmfsk_demod(states,rx) end states.nin = next_nin; - + states.norm_rx_timing = norm_rx_timing; + norm_rx_timing %'Even' and 'Odd' manchester bitstream. % We'll figure out which to produce later rx_even = zeros(1,nbits); rx_odd = zeros(1,nbits); apeven = 0; apodd = 0; - + sample_offset = (Ts/2)+Ts+rx_timing-1; + symsamp = zeros(1,nsym); + % Figure out the bits of the 'even' and 'odd' ME streams % Also sample rx_filt offset by what fine timing determined along the way % Note: ii is a zero-indexed array pointer, for less mind-breaking c portage @@ -157,7 +160,7 @@ function [rx_bits states] = fmfsk_demod(states,rx) mdiff = lastv-currv; lastv = currv; mbit = mdiff>0; - + symsamp(ii+1) = currv; if mod(ii,2)==1 apeven += abs(mdiff); rx_even( floor(ii/2)+1 ) = mbit; @@ -166,6 +169,7 @@ function [rx_bits states] = fmfsk_demod(states,rx) rx_odd( floor(ii/2)+1 ) = mbit; end end + states.symsamp = symsamp; % Decide on the correct ME alignment if(apeven>apodd) @@ -174,7 +178,7 @@ function [rx_bits states] = fmfsk_demod(states,rx) rx_bits = rx_odd; end - states.lastint = rx_filt(sample_offset+((nsym-1)*Ts)+1); + states.lastint = lastv; endfunction % run_sim copypasted from fsk_horus.m @@ -311,14 +315,6 @@ function fmfsk_run_sim(EbNodB,timing_offset=0,de=0,of=0,hpf=0) rx_bits_buf(1:nbit) = rx_bits_buf(nbit+1:2*nbit); rx_bits_buf(nbit+1:2*nbit) = rx_bits; rx_bits_log = [rx_bits_log rx_bits]; - %rx_bits_sd_log = [rx_bits_sd_log states.rx_bits_sd]; - - %norm_rx_timing_log = [norm_rx_timing_log states.norm_rx_timing]; - %x_log = [x_log states.x]; - %timing_nl_log = [timing_nl_log states.timing_nl]; - %f_int_resample_log = [f_int_resample_log abs(states.f_int_resample(:,:))]; - %f_log = [f_log; states.f]; - %EbNodB_log = [EbNodB_log states.EbNodB]; end diff --git a/codec2-dev/octave/tfmfsk.m b/codec2-dev/octave/tfmfsk.m new file mode 100644 index 00000000..13c317f3 --- /dev/null +++ b/codec2-dev/octave/tfmfsk.m @@ -0,0 +1,517 @@ +% tfsk.m +% Author: Brady O'Brien 8 February 2016 + + + +% Copyright 2016 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 . + + +% Octave script to check c port of mancyfsk/fmfsk against the fmfsk.m +% +#{ + + FMFSK Modem automated test instructions: + + 1. Use cmake to build in debug mode to ensure unittest/tfsk is built: + + $ cd ~/codec2 + $ rm -Rf build_linux && mkdir build_linux + $ cd build_linux + $ cmake -DCMAKE_BUILD_TYPE=Debug .. + $ make + + 2 - Change tfsk_location below if required + 3 - Ensure Octave packages signal and parallel are installed + 4 - Start Octave and run tfsk.m. It will perform all tests automatically + +#} + +%tfsk executable path/file +global tfsk_location = '../build_linux/unittest/tfmfsk'; + + + +fmfsk +pkg load signal; +pkg load parallel; +graphics_toolkit('gnuplot'); + + +global mod_pass_fail_maxdiff = 1e-3/5000; + +function mod = fmfsk_mod_c(Fs,Rs,bits) + global tfsk_location; + %command to be run by system to launch the modulator + command = sprintf('%s M %d %d fsk_mod_ut_bitvec fsk_mod_ut_modvec fmfsk_mod_ut_log.txt',tfsk_location,Fs,Rs); + %save input bits into a file + bitvecfile = fopen('fsk_mod_ut_bitvec','wb+'); + fwrite(bitvecfile,bits,'uint8'); + fclose(bitvecfile); + + %run the modulator + system(command); + + modvecfile = fopen('fsk_mod_ut_modvec','rb'); + mod = fread(modvecfile,'single'); + fclose(modvecfile); + +endfunction + + +%Compare 2 vectors, fail if they are not close enough +function pass = vcompare(vc,voct,vname,tname,tol,pnum) + + %Get delta of vectors + dvec = abs(abs(vc)-abs(voct)); + + %Normalize difference + dvec = dvec ./ abs(max(abs(voct))+1e-8); + + maxdvec = abs(max(dvec)); + pass = maxdvec=states.nin + ninold = states.nin; + [bitbuf,states] = fmfsk_demod(states, modin(1:states.nin)); + modin=modin(ninold+1:length(modin)); + obits = [obits bitbuf]; + + o_norm_rx_timing = [o_norm_rx_timing states.norm_rx_timing]; + o_symsamp = [o_symsamp states.symsamp]; + o_rx_filt = [o_rx_filt states.rx_filt]; + + end + + close all + + % One part-per-thousand allowed on important parameters + pass = 1; + + figure(8) + plot((1:length(t_rx_filt)),t_rx_filt,(1:length(o_rx_filt)),o_rx_filt); + figure(9) + plot((1:length(t_symsamp)),t_symsamp,(1:length(o_symsamp)),o_symsamp); + + %vcompare(t_rx_filt,o_rx_filt,.001,8); + vcompare(t_norm_rx_timing,o_norm_rx_timing,'norm rx timing',.001,9); + vcompare(t_symsamp,o_symsamp,'symsamp',.001,10); + + assert(pass); + diffpass = sum(xor(obits,bits'))<4; + diffbits = sum(xor(obits,bits')); + + + + if diffpass==0 + printf('\n***bitcompare test failed test %s diff %d\n\n',tname,sum(xor(obits,bits'))) + figure(15) + plot(xor(obits,bits')) + title(sprintf('Bitcompare failure test %s',tname)) + end + + pass = pass && diffpass; + + + test_stats.pass = pass; + test_stats.diff = sum(xor(obits,bits')); + test_stats.cbits = bits'; + test_stats.obits = obits; + +endfunction + +function [dmod,cmod,omod,pass] = fmfsk_mod_test(Fs,Rs,bits,tname,M=2) + global mod_pass_fail_maxdiff; + %Run the C modulator + cmod = fmfsk_mod_c(Fs,Rs,bits); + %Set up and run the octave modulator + states.M = M; + states = fmfsk_init(Fs,Rs); + + + omod = fmfsk_mod(states,bits)'; + + dmod = cmod-omod; + pass = max(dmod)<(mod_pass_fail_maxdiff*length(dmod)); + if !pass + printf('Mod failed test %s!\n',tname); + end +endfunction + +% Random bit modulator test +% Pass random bits through the modulators and compare +function pass = test_mod_fdvbcfg_randbits + rand('state',1); + randn('state',1); + bits = rand(1,19200)>.5; + [dmod,cmod,omod,pass] = fmfsk_mod_test(48000,2400,bits,"mod fdvbcfg randbits"); + + if(!pass) + figure(1) + plot(dmod) + title("Difference between octave and C mod impl"); + end + +endfunction + +% run_sim copypasted from fsk_horus.m +% simulation of tx and rx side, add noise, channel impairments ---------------------- + +function stats = tfmfsk_run_sim(EbNodB,timing_offset=0,de=0,of=0,hpf=0,M=2) + test_frame_mode = 2; + frames = 70; + %EbNodB = 3; + %timing_offset = 0.0; % see resample() for clock offset below + %fading = 0; % modulates tx power at 2Hz with 20dB fade depth, + % to simulate balloon rotating at end of mission + df = 0; % tx tone freq drift in Hz/s + dA = 1; % amplitude imbalance of tones (note this affects Eb so not a gd idea) + + more off + rand('state',1); + randn('state',1); + + Fs = 48000; + Rbit = 2400; + + % ---------------------------------------------------------------------- + + fm_states.pre_emp = 0; + fm_states.de_emp = de; + fm_states.Ts = Fs/(Rbit*2); + fm_states.Fs = Fs; + fm_states.fc = Fs/4; + fm_states.fm_max = 3E3; + fm_states.fd = 5E3; + fm_states.output_filter = of; + fm_states = analog_fm_init(fm_states); + + % ---------------------------------------------------------------------- + + states = fmfsk_init(Fs,Rbit); + + states.verbose = 0x1; + Rs = states.Rs; + nsym = states.nsym; + Fs = states.Fs; + nbit = states.nbit; + + EbNo = 10^(EbNodB/10); + variance = states.Fs/(states.Rb*EbNo); + + % set up tx signal with payload bits based on test mode + + if test_frame_mode == 1 + % test frame of bits, which we repeat for convenience when BER testing + test_frame = round(rand(1, states.nbit)); + tx_bits = []; + for i=1:frames+1 + tx_bits = [tx_bits test_frame]; + end + end + if test_frame_mode == 2 + % random bits, just to make sure sync algs work on random data + tx_bits = round(rand(1, states.nbit*(frames+1))); + end + if test_frame_mode == 3 + % repeating sequence of all symbols + % great for initial test of demod if nothing else works, + % look for this pattern in rx_bits + + % ...10101... + tx_bits = zeros(1, states.nbit*(frames+1)); + tx_bits(1:2:length(tx_bits)) = 1; + + end + + [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)); + if(timing_offset>0) + tx = resample(tx, 1000, 1001); % simulated 1000ppm sample clock offset + end + + + noise = sqrt(variance)*randn(length(tx),1); + rx = tx + noise'; + + %Demod by analog fm + rx = analog_fm_demod(fm_states, rx); + + %High-pass filter to simulate the FM radios + if hpf>0 + printf("high-pass filtering!\n") + rx = filter(b,a,rx); + end + %rx = real(rx); + %b1 = fir2(100, [0 4000 5200 48000]/48000, [1 1 0.5 0.5]); + %rx = filter(b1,1,rx); + %[b a] = cheby2(6,40,[3000 6000]/(Fs/2)); + %rx = filter(b,a,rx); + %rx = sign(rx); + %rx(find (rx > 1)) = 1; + %rx(find (rx < -1)) = -1; + + % dump simulated rx file + + timing_offset_samples = round(timing_offset*states.Ts); + st = 1 + timing_offset_samples; + rx_bits_buf = zeros(1,2*nbit); + x_log = []; + timing_nl_log = []; + norm_rx_timing_log = []; + f_int_resample_log = []; + f_log = []; + EbNodB_log = []; + rx_bits_log = []; + rx_bits_sd_log = []; + + + test_name = sprintf("tfmfsk run sim EbNodB:%d frames:%d timing_offset:%d df:%d",EbNodB,frames,timing_offset,df); + tstats = fmfsk_demod_xt(Fs,Rbit,rx',test_name,M); + + pass = tstats.pass; + obits = tstats.obits; + cbits = tstats.cbits; + + % Figure out BER of octave and C modems + bitcnt = length(tx_bits); + rx_bits = obits; + ber = 1; + ox = 1; + for offset = (1:100) + nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))); + bern = nerr/(bitcnt-offset); + if(bern < ber) + ox = offset; + best_nerr = nerr; + end + ber = min([ber bern]); + end + offset = ox; + bero = ber; + ber = 1; + rx_bits = cbits; + ox = 1; + for offset = (1:100) + nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))); + bern = nerr/(bitcnt-offset); + if(bern < ber) + ox = offset; + best_nerr = nerr; + end + ber = min([ber bern]); + end + offset = ox; + berc = ber; + stats.berc = berc; + stats.bero = bero; + + % non-coherent BER theory calculation + % It was complicated, so I broke it up + + ms = 2; + ns = (1:ms-1); + as = (-1).^(ns+1); + bs = (as./(ns+1)); + + cs = ((ms-1)./ns); + + ds = ns.*log2(ms); + es = ns+1; + fs = exp( -(ds./es)*EbNo ); + + thrncoh = ((ms/2)/(ms-1)) * sum(bs.*((ms-1)./ns).*exp( -(ds./es)*EbNo )); + + stats.thrncoh = thrncoh; + stats.pass = pass; + endfunction + + +function pass = ebno_battery_test(timing_offset,fading,df,dA,M) + %Range of EbNodB over which to test + ebnodbrange = (3:2:13); + ebnodbs = length(ebnodbrange); + + mode = 2; + %Replication of other parameters for parcellfun + modev = repmat(mode,1,ebnodbs); + timingv = repmat(timing_offset,1,ebnodbs); + fadingv = repmat(fading,1,ebnodbs); + dfv = repmat(df,1,ebnodbs); + dav = repmat(dA,1,ebnodbs); + mv = repmat(M,1,ebnodbs); + statv = pararrayfun(floor(1.25*nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,mv); + %statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,mv); + + passv = zeros(1,length(statv)); + for ii=(1:length(statv)) + passv(ii)=statv(ii).pass; + end + + %All pass flags are '1' + pass = sum(passv)>=length(passv); + %and no tests died + pass = pass && length(passv)==ebnodbs; + passv; + assert(pass) +endfunction + +%Test with and without sample clock offset +function pass = test_timing_var(df,dA,M) + pass = ebno_battery_test(1,0,df,dA,M) + assert(pass) + pass = pass && ebno_battery_test(0,0,df,dA,M) + assert(pass) +endfunction + +%Test with and without 1 Hz/S freq drift +function pass = test_drift_var(M) + pass = test_timing_var(1,1,M) + assert(pass) + pass = pass && test_timing_var(0,1,M) + assert(pass) +endfunction + +function pass = test_fsk_battery() + pass = test_mod_horuscfg_randbits; + assert(pass) + pass = pass && test_mod_horuscfgm4_randbits; + assert(pass) + pass = pass && test_drift_var(4); + assert(pass) + pass = pass && test_drift_var(2); + assert(pass) + if pass + printf("***** All tests passed! *****\n"); + end +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); + + + statv = pararrayfun(floor(1.25*nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv); + %statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv); + + for ii = (1:length(statv)) + stat = statv(ii); + berc(ii)=stat.berc; + bero(ii)=stat.bero; + berinc(ii)=stat.thrncoh; + end + clf; + figure(M) + + semilogy(ebnodbrange, berinc,sprintf('r;%dFSK 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)) + hold off; + grid("minor"); + axis([min(ebnodbrange) max(ebnodbrange) 1E-5 1]) + legend("boxoff"); + xlabel("Eb/No (dB)"); + ylabel("Bit Error Rate (BER)") + +endfunction + + +%test_fsk_battery +%plot_fsk_bers(2) +%plot_fsk_bers(4) diff --git a/codec2-dev/src/fmfsk.c b/codec2-dev/src/fmfsk.c index e9478b50..bc3357c7 100644 --- a/codec2-dev/src/fmfsk.c +++ b/codec2-dev/src/fmfsk.c @@ -105,27 +105,6 @@ uint32_t fmfsk_nin(struct FMFSK *fmfsk){ * float mod_out[] - Buffer for N samples of modulated FMFSK * uint8_t tx_bits[] - Buffer containing Nbits unpacked bits */ -//void fmfsk_mod(struct FMFSK *fmfsk, float fmfsk_out[],uint8_t bits_in[]){ -// int i,j; -// int Ts = fmfsk->Ts; -// int N = fmfsk->N; -// -// for(i=0;iTs; - int Fs = fmfsk->Ts; + int Fs = fmfsk->Fs; int Rs = fmfsk->Rs; int nin = fmfsk->nin; int N = fmfsk->N; @@ -183,6 +162,8 @@ void fmfsk_demod(struct FMFSK *fmfsk, uint8_t rx_bits[],float fmfsk_in[]){ memcpy(&oldsamps[0] , &oldsamps[nmem-nold], sizeof(float)*nold); memcpy(&oldsamps[nold], &fmfsk_in[0] , sizeof(float)*nin ); + fprintf(stderr,"nold:%d nmem:%d nin:%d\n",nold,nmem,nin); + /* Allocate memory for filtering */ float *rx_filt = alloca(sizeof(float)*(nsym+1)*Ts); @@ -190,7 +171,7 @@ void fmfsk_demod(struct FMFSK *fmfsk, uint8_t rx_bits[],float fmfsk_in[]){ for(i=0; i<(nsym+1)*Ts; i++){ t=0; /* Integrate over some samples */ - for(j=i;j .5) next_nin += Ts/2; - if(norm_rx_timing < .5) + if(norm_rx_timing < -.5) next_nin -= Ts/2; fmfsk->nin = next_nin; @@ -250,6 +233,7 @@ void fmfsk_demod(struct FMFSK *fmfsk, uint8_t rx_bits[],float fmfsk_in[]){ for(i=0; i0 ? 1 : 0; lastv = currv; @@ -283,5 +267,6 @@ void fmfsk_demod(struct FMFSK *fmfsk, uint8_t rx_bits[],float fmfsk_in[]){ fmfsk->lodd = lastv; modem_probe_samp_f("t_norm_rx_timing",&norm_rx_timing,1); - modem_probe_samp_f("t_rx_filt",rx_filt,nsym*Ts); + modem_probe_samp_f("t_rx_filt",rx_filt,(nsym+1)*Ts); + fprintf(stderr,"rxfiltl %d\n",(nsym+1)*Ts); } diff --git a/codec2-dev/unittest/tfmfsk.c b/codec2-dev/unittest/tfmfsk.c index 877e639e..b1359624 100644 --- a/codec2-dev/unittest/tfmfsk.c +++ b/codec2-dev/unittest/tfmfsk.c @@ -168,6 +168,8 @@ int main(int argc,char *argv[]){ bitbuf = malloc(sizeof(uint8_t)*fmfsk->nbit); /* Demod-only test */ if(test_type == TEST_DEMOD){ + + fprintf(stderr,"%d\n",(fmfsk->N+fmfsk->Ts*2)); while( fread(modbuf,sizeof(float),fmfsk_nin(fmfsk),fin) == fmfsk_nin(fmfsk) ){ fmfsk_demod(fmfsk,bitbuf,modbuf); fwrite(bitbuf,sizeof(uint8_t),fmfsk->nbit,fout);