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.
% 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
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
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;
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
% 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,...
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,...
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
+
+
projects / freetel-svn-tracking.git / commitdiff