From bc84bea8b6c3209de7c2b1750c66b979e823a81c Mon Sep 17 00:00:00 2001 From: drowe67 Date: Tue, 16 Dec 2014 05:40:03 +0000 Subject: [PATCH] first pass a GMSK sim, but 5dB off theoretical git-svn-id: https://svn.code.sf.net/p/freetel/code@1971 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fm.m | 14 ++- codec2-dev/octave/fsk.m | 6 +- codec2-dev/octave/gmsk.m | 247 +++++++++++++++++++++++++++++++++++++-- 3 files changed, 250 insertions(+), 17 deletions(-) diff --git a/codec2-dev/octave/fm.m b/codec2-dev/octave/fm.m index 0896dff7..9a30fd01 100644 --- a/codec2-dev/octave/fm.m +++ b/codec2-dev/octave/fm.m @@ -11,13 +11,12 @@ function fm_states = analog_fm_init(fm_states) Fs = fm_states.Fs; FsOn2 = Fs/2; fm_max = fm_states.fm_max; % max modulation freq - fm_states.fc = 24E3; % carrier frequency fd = fm_states.fd; % (max) deviation fm_states.m = fd/fm_max; % modulation index fm_states.Bfm = Bfm = 2*(fd+fm_max); % Carson's rule for FM signal bandwidth fm_states.tc = tc = 50E-6; fm_states.prede = [1 -(1 - 1/(tc*Fs))]; % pre/de emp filter coeffs - + % Select length of filter to be an integer number of symbols to % assist with "fine" timing offset estimation. Set Ts to 1 for % analog modulation. @@ -43,9 +42,8 @@ function fm_states = analog_fm_init(fm_states) % demoduator output filter to limit us to fm_max (e.g. 3kHz) - fc = (fm_max)/(FsOn2); + fc = fm_max/(FsOn2); fm_states.bout = firls(ncoeffs,[0 0.95*fc 1.05*fc 1], [1 1 0.01 0.01]); - endfunction @@ -83,7 +81,9 @@ function [rx_out rx_bb] = analog_fm_demod(fm_states, rx) rx_bb = filter(fm_states.bin,1,rx_bb); rx_bb_diff = [ 1 rx_bb(2:nsam) .* conj(rx_bb(1:nsam-1))]; rx_out = (1/wd)*atan2(imag(rx_bb_diff),real(rx_bb_diff)); - rx_out = filter(fm_states.bout,1,rx_out); + if fm_states.output_filter + rx_out = filter(fm_states.bout,1,rx_out); + end if fm_states.de_emp rx_out = filter(1,fm_states.prede,rx_out); end @@ -98,10 +98,12 @@ function sim_out = analog_fm_test(sim_in) Fs = fm_states.Fs = 96000; fm_max = fm_states.fm_max = 3E3; fd = fm_states.fd = 5E3; + fm_states.fc = 24E3; fm_states.pre_emp = pre_emp = sim_in.pre_emp; fm_states.de_emp = de_emp = sim_in.de_emp; fm_states.Ts = 1; + fm_states.output_filter = 1; fm_states = analog_fm_init(fm_states); sim_out.Bfm = fm_states.Bfm; @@ -246,7 +248,7 @@ function run_fm_single sim_in.pre_emp = 0; sim_in.de_emp = 0; - sim_in.CNdB = 10; + sim_in.CNdB = 20; sim_out = analog_fm_test(sim_in); end diff --git a/codec2-dev/octave/fsk.m b/codec2-dev/octave/fsk.m index a321fc3f..74fee23f 100644 --- a/codec2-dev/octave/fsk.m +++ b/codec2-dev/octave/fsk.m @@ -50,8 +50,10 @@ function sim_out = fsk_ber_test(sim_in) fm_states.de_emp = 0; fm_states.Ts = Ts; fm_states.Fs = Fs; + fm_states.fc = Fs/4; fm_states.fm_max = 3E3; fm_states.fd = 5E3; + fm_states.output_filter = 1; fm_states = analog_fm_init(fm_states); end @@ -168,10 +170,6 @@ function sim_out = fsk_ber_test(sim_in) stem(abs(mark_int(1:100))); subplot(212) stem(abs(space_int(1:100))); - - figure(5) - clf - plot(error_positions); end if verbose diff --git a/codec2-dev/octave/gmsk.m b/codec2-dev/octave/gmsk.m index a873313b..33b86238 100644 --- a/codec2-dev/octave/gmsk.m +++ b/codec2-dev/octave/gmsk.m @@ -1,9 +1,17 @@ % gmsk.m +% David Rowe Dec 2014 +% +% GMSK modem simulation +% +% [ ] plot eye diagram +% [ ] BER curves with theoretical results +% [ ] spectrum - will it pass thru a HT? -% From: https://github.com/on1arf/gmsk/blob/master/gmskmodem_codec2/API/a_dspstuff.h, which is in turn -% from Jonathan G4KLX +% Filter coeffs From: +% https://github.com/on1arf/gmsk/blob/master/gmskmodem_codec2/API/a_dspstuff.h, +% which is in turn from Jonathan G4KLX. The demod coeffs LPF noise -gmsk_mod_coeff = [... +global gmsk_mod_coeff = [... 6.455906007234699e-014, 1.037067381285011e-012, 1.444835156335346e-011,... 1.745786683011439e-010, 1.829471305298363e-009, 1.662729407135958e-008,... 1.310626978701910e-007, 8.959797186410516e-007, 5.312253663302771e-006,... @@ -19,7 +27,7 @@ gmsk_mod_coeff = [... 1.829471305298363e-009, 1.745786683011426e-010, 1.444835156335356e-011,... 1.037067381285011e-012, 6.455906007234699e-014]; -gmsk_demod_coeff = [... +global gmsk_demod_coeff = [... -0.000153959924563, 0.000000000000000, 0.000167227768379, 0.000341615513437,... 0.000513334449696, 0.000667493753523, 0.000783901543032, 0.000838293462576,... 0.000805143268199, 0.000661865814384, 0.000393913058926, -0.000000000000000,... @@ -47,9 +55,234 @@ gmsk_demod_coeff = [... 0.000783901543032, 0.000667493753523, 0.000513334449696, 0.000341615513437,... 0.000167227768379, 0.000000000000000, -0.000153959924563]; -function tx = gmsk_mod(tx_bits) - % filter - % FM modulate, 1.2 kHz deviation, break out fm functions, specify Fs and deviation +rand('state',1); +randn('state',1); +graphics_toolkit ("gnuplot"); +fm; + +function gmsk_states = gmsk_init + gmsk_states.Fs = 48000; + gmsk_states.Rs = 4800; + gmsk_states.M = gmsk_states.Fs/gmsk_states.Rs; + global gmsk_mod_coeff; + global gmsk_demod_coeff; + + fm_states.Fs = gmsk_states.Fs; + fm_states.fc = 0; + fm_max = fm_states.fm_max = 2400; + fd = fm_states.fd = 1200; + fm_states.Ts = gmsk_states.M; + fm_states.pre_emp = fm_states.de_emp = 0; + fm_states.output_filter = 1; + gmsk_states.fm_states = analog_fm_init(fm_states); + + [x i_mod] = max(gmsk_mod_coeff); + [x i_demod] = max(gmsk_demod_coeff); + gmsk_states.filter_delay = i_mod + i_demod; + + gmsk_states.Toff = -2; + gmsk_states.dsam = dsam = gmsk_states.filter_delay + gmsk_states.fm_states.nsym_delay*gmsk_states.M; + gmsk_states.dsym = floor(dsam/gmsk_states.M); +endfunction + + +function [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits) + M = gmsk_states.M; + nsym = length(tx_bits); + nsam = nsym*M; + global gmsk_mod_coeff; + + tx_symbols = zeros(1,nsam); + tx_symbols(1:M:nsam) = -1 + 2*tx_bits; + + tx_filt = filter(gmsk_mod_coeff, 1, tx_symbols); + tx_filt = tx_filt / max(tx_filt); + tx = analog_fm_mod(gmsk_states.fm_states, tx_filt); + % work out delays of filter to align bits % plot eye diagrams, BERcurves, theoretical results, spectrum - will it pass thru a HT? endfunction + + +function [rx_bits rx_filt] = gmsk_demod(gmsk_states, rx) + M = gmsk_states.M; + nsam = length(rx); + nsym = floor(nsam/M); + global gmsk_demod_coeff; + + [rx_bb fm_bb] = analog_fm_demod(gmsk_states.fm_states, rx); + rx_filt = filter(gmsk_demod_coeff, 1, rx_bb); + + rx_bits = zeros(1,nsym); + Toff = gmsk_states.Toff; + dsam = gmsk_states.dsam; + rx_bits = rx_filt(dsam+1+Toff:M:length(rx_filt)) > 0; + + figure(1) + plot(fm_bb) +endfunction + + +function sim_out = gmsk_test(sim_in) + nsym = sim_in.nsym; + EbNodB = sim_in.EbNodB; + verbose = sim_in.verbose; + + gmsk_states = gmsk_init; + M = gmsk_states.M; + Fs = gmsk_states.Fs; + Rs = gmsk_states.Rs; + Bfm = gmsk_states.fm_states.Bfm; + dsam = gmsk_states.dsam; + + for ne = 1:length(EbNodB) + aEbNodB = EbNodB(ne); + EbNo = 10^(aEbNodB/10); + variance = Fs/(Rs*EbNo); + + tx_bits = round(rand(1, nsym)); + [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits); + nsam = length(tx); + + noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); + rx = tx + noise; + + [rx_bits rx_filt] = gmsk_demod(gmsk_states, rx); + + l = length(rx_bits); + error_positions = xor(rx_bits(1:l), tx_bits(1:l)); + Nerrs = sum(error_positions); + TERvec(ne) = Nerrs; + BERvec(ne) = Nerrs/l; + + if verbose > 0 + printf("EbNo dB: %3.1f Nerrs: %d BER: %f BER Theory: %f\n", aEbNodB, Nerrs, BERvec(ne), 0.5*erfc(sqrt(0.75*EbNo))); + end + + if verbose > 1 + + figure(2); + nplot = 16; + clf; + subplot(211) + stem(tx_symbols(1:M:nplot*M)) + axis([1 nplot -1 1]) + subplot(212) + Toff = gmsk_states.Toff; + stem(rx_filt(dsam+1+Toff:M:dsam+nplot*M+Toff)) + axis([1 nplot -1 1]) + + figure(3); + clf; + subplot(211) + plot(tx_filt(21+1:21+nplot*M)) + axis([1 nplot*M -1 1]); + subplot(212) + d = gmsk_states.filter_delay + gmsk_states.fm_states.nsym_delay*M; + plot(rx_filt(d+1:d+nplot*M)) + axis([1 nplot*M -1 1]); + + figure(4); + clf + subplot(211); + f = fft(rx); + Tx = 20*log10(abs(f)); + plot(Tx) + grid; + title('GMSK Demodulator Input Spectrum'); + axis([1 Fs/2 0 80]) + + subplot(212) + f = fft(tx); + f = f(1:length(f)/2); + cs = cumsum(abs(f).^2); + plot(cs) + hold on; + x = 0.99; + tots = x*sum(abs(f).^2); + xpercent_pwr = find(cs > tots); + bw = 2*xpercent_pwr(1); + plot([1 Fs/2],[tots tots],'r') + plot([bw/2 bw/2],[0 tots],'r') + hold off; + title("Cumulative Power"); + grid; + + printf("Bfm: %4.0fHz %3.0f%% power bandwidth %4.0fHz = %3.2f*Rb\n", Bfm, x*100, bw, bw/Rs); + + figure(5) + clf; + subplot(211) + Tx_filt = 20*log10(abs(fft(tx_filt))); + plot(Tx_filt) + axis([1 10000 0 80]); + grid; + subplot(212) + Rx_filt = 20*log10(abs(fft(rx_filt))); + plot(Rx_filt) + axis([1 10000 0 80]); + grid; + end + end + + sim_out.TERvec = TERvec; + sim_out.BERvec = BERvec; + sim_out.Rs = gmsk_states.Rs; +endfunction + +function run_gmsk_single + sim_in.nsym = 4800; + sim_in.EbNodB = 120; + sim_in.verbose = 2; + + sim_out = gmsk_test(sim_in); +endfunction + +function run_gmsk_curves + sim_in.nsym = 4800; + sim_in.EbNodB = 5:10; + sim_in.verbose = 1; + + gmsk_sim = gmsk_test(sim_in); + Rs = gmsk_sim.Rs; + EbNo = 10 .^ (sim_in.EbNodB/10); + alpha = 0.75; + gmsk_theory.BERvec = 0.5*erfc(sqrt(alpha*EbNo)); + + % BER v Eb/No curves + + figure(1); + clf; + semilogy(sim_in.EbNodB, gmsk_theory.BERvec,'r;GMSK theory;') + hold on; + semilogy(sim_in.EbNodB, gmsk_sim.BERvec,'g;GMSK sim;') + hold off; + grid("minor"); + axis([min(sim_in.EbNodB) max(sim_in.EbNodB) 1E-4 1]) + legend("boxoff"); + xlabel("Eb/No (dB)"); + ylabel("Bit Error Rate (BER)") + + % BER v C/No (1 Hz noise BW and Eb=C/Rs=1/Rs) + % Eb/No = (C/Rs)/(1/(N/B)) + % C/N = (Eb/No)*(Rs/B) + + RsOnB_dB = 10*log10(Rs/1); + figure(2); + clf; + semilogy(sim_in.EbNodB+RsOnB_dB, gmsk_theory.BERvec,'r;GMSK theory;') + hold on; + semilogy(sim_in.EbNodB+RsOnB_dB, gmsk_sim.BERvec,'g;GMSK sim;') + hold off; + grid("minor"); + axis([min(sim_in.EbNodB+RsOnB_dB) max(sim_in.EbNodB+RsOnB_dB) 1E-4 1]) + legend("boxoff"); + xlabel("C/No for Rs=4800 bit/s and 1 Hz noise bandwidth (dB)"); + ylabel("Bit Error Rate (BER)") + +endfunction + +run_gmsk_single +%run_gmsk_curves + + -- 2.25.1