--- /dev/null
+% fm.m
+% David Rowe Dec 2014
+%
+% Analog FM Octave simulation functions.
+
+1;
+
+function fm_states = analog_fm_init(fm_states)
+
+ % FM modulator constants
+
+ 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.
+
+ Ts = fm_states.Ts;
+ desired_ncoeffs = 200;
+ ncoeffs = floor(desired_ncoeffs/Ts+1)*Ts;
+
+ % "coarse" timing offset is half filter length, we have two filters.
+ % This is the delay the two filters introduce, so we need to adjust
+ % for this when comparing tx to trx bits for BER calcs.
+
+ fm_states.nsym_delay = ncoeffs/Ts;
+
+ % input filter gets rid of excess noise before demodulator, as too much
+ % noise causes atan2() to jump around, e.g. -pi to pi. However this
+ % filter can cause harmonic distortion at very high SNRs, as it knocks out
+ % some of the FM signal spectra. This filter isn't really required for high
+ % SNRs > 20dB.
+
+ fc = (Bfm/2)/(FsOn2);
+ fm_states.bin = firls(ncoeffs,[0 fc*(1-0.05) fc*(1+0.05) 1],[1 1 0.01 0.01]);
+
+ % demoduator output filter to limit us to fm_max (e.g. 3kHz)
+
+ fc = (fm_max)/(FsOn2);
+ fm_states.bout = firls(ncoeffs,[0 0.95*fc 1.05*fc 1], [1 1 0.01 0.01]);
+
+endfunction
+
+
+function tx = analog_fm_mod(fm_states, mod)
+ Fs = fm_states.Fs;
+ fc = fm_states.fc; wc = 2*pi*fc/Fs;
+ fd = fm_states.fd; wd = 2*pi*fd/Fs;
+ nsam = length(mod);
+
+ if fm_states.pre_emp
+ mod = filter(fm_states.prede,1,mod);
+ mod = mod/max(mod); % AGC to set deviation
+ end
+
+ tx_phase = 0;
+ tx = zeros(1,nsam);
+
+ for i=0:nsam-1
+ w = wc + wd*mod(i+1);
+ tx_phase = tx_phase + w;
+ tx_phase = tx_phase - floor(tx_phase/(2*pi))*2*pi;
+ tx(i+1) = exp(j*tx_phase);
+ end
+endfunction
+
+
+function [rx_out rx_bb] = analog_fm_demod(fm_states, rx)
+ Fs = fm_states.Fs;
+ fc = fm_states.fc; wc = 2*pi*fc/Fs;
+ fd = fm_states.fd; wd = 2*pi*fd/Fs;
+ nsam = length(rx);
+ t = 0:(nsam-1);
+
+ rx_bb = rx .* exp(-j*wc*t); % down to complex baseband
+ 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.de_emp
+ rx_out = filter(1,fm_states.prede,rx_out);
+ end
+endfunction
+
+
+function sim_out = analog_fm_test(sim_in)
+ nsam = sim_in.nsam;
+ CNdB = sim_in.CNdB;
+ verbose = sim_in.verbose;
+
+ Fs = fm_states.Fs = 96000;
+ fm_max = fm_states.fm_max = 3E3;
+ fd = fm_states.fd = 5E3;
+
+ 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 = analog_fm_init(fm_states);
+ sim_out.Bfm = fm_states.Bfm;
+
+ Bfm = fm_states.Bfm;
+ m = fm_states.m; tc = fm_states.tc;
+ t = 0:(nsam-1);
+
+ fm = 1000; wm = 2*pi*fm/fm_states.Fs;
+
+ % start simulation
+
+ for ne = 1:length(CNdB)
+
+ % work out the variance we need to obtain our C/N in the bandwidth
+ % of the FM demod. The gaussian generator randn() generates noise
+ % with a bandwidth of Fs
+
+ aCNdB = CNdB(ne);
+ CN = 10^(aCNdB/10);
+ variance = Fs/(CN*Bfm);
+
+ % FM Modulator -------------------------------
+
+ mod = sin(wm*t);
+ tx = analog_fm_mod(fm_states, mod);
+
+ % Channel ---------------------------------
+
+ noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam));
+ rx = tx + noise;
+
+ % FM Demodulator
+
+ [rx_out rx_bb] = analog_fm_demod(fm_states, rx);
+
+ % notch out test tone
+
+ w = 2*pi*fm/Fs; beta = 0.99;
+ rx_notch = filter([1 -2*cos(w) 1],[1 -2*beta*cos(w) beta*beta], rx_out);
+
+ % measure power with and without test tone to determine S+N and N
+
+ settle = 1000; % filter settling time, to avoid transients
+ nsettle = nsam - settle;
+
+ sinad = (rx_out(settle:nsam) * rx_out(settle:nsam)')/nsettle;
+ nad = (rx_notch(settle:nsam) * rx_notch(settle:nsam)')/nsettle;
+
+ snr = (sinad-nad)/nad;
+ sim_out.snrdB(ne) = 10*log10(snr);
+
+ % Theory from FMTutorial.pdf, Lawrence Der, Silicon labs paper
+
+ snr_theory_dB = aCNdB + 10*log10(3*m*m*(m+1));
+ fx = 1/(2*pi*tc); W = fm_max;
+ I = (W/fx)^3/(3*((W/fx) - atan(W/fx)));
+ I_dB = 10*log10(I);
+
+ sim_out.snr_theorydB(ne) = snr_theory_dB;
+ sim_out.snr_theory_pre_dedB(ne) = snr_theory_dB + I_dB;
+
+ if verbose > 1
+ printf("modn index: %2.1f Bfm: %.0f Hz\n", m, Bfm);
+ end
+
+ if verbose > 0
+ printf("C/N: %4.1f SNR: %4.1f dB THEORY: %4.1f dB or with pre/de: %4.1f dB\n",
+ aCNdB, 10*log10(snr), snr_theory_dB, snr_theory_dB+I_dB);
+ end
+
+ if verbose > 1
+ figure(1)
+ subplot(211)
+ plot(20*log10(abs(fft(rx))))
+ title('FM Modulator Output Spectrum');
+ axis([1 length(tx) 0 100]);
+ subplot(212)
+ Rx_bb = 20*log10(abs(fft(rx_bb)));
+ plot(Rx_bb)
+ axis([1 length(tx) 0 100]);
+ title('FM Demodulator (baseband) Input Spectrum');
+
+ figure(2)
+ subplot(211)
+ plot(rx_out(settle:nsam))
+ axis([1 4000 -1 1])
+ subplot(212)
+ Rx = 20*log10(abs(fft(rx_out(settle:nsam))));
+ plot(Rx(1:10000))
+ axis([1 10000 0 100]);
+ end
+
+ end
+
+endfunction
+
+
+function run_fm_curves
+ sim_in.nsam = 96000;
+ sim_in.verbose = 1;
+ sim_in.pre_emp = 0;
+ sim_in.de_emp = 0;
+ sim_in.CNdB = -4:2:20;
+
+ sim_out = analog_fm_test(sim_in);
+
+ figure(1)
+ clf
+ plot(sim_in.CNdB, sim_out.snrdB,"r;FM Simulated;");
+ hold on;
+ plot(sim_in.CNdB, sim_out.snr_theorydB,"g;FM Theory;");
+ plot(sim_in.CNdB, sim_in.CNdB,"b; SSB Theory;");
+ hold off;
+ grid("minor");
+ xlabel("FM demod input C/N (dB)");
+ ylabel("FM demod output S/N (dB)");
+ legend("boxoff");
+
+ % C/No curves
+
+ Bfm_dB = 10*log10(sim_out.Bfm);
+ Bssb_dB = 10*log10(3000);
+
+ figure(2)
+ clf
+ plot(sim_in.CNdB + Bfm_dB, sim_out.snrdB,"r;FM Simulated;");
+ hold on;
+ plot(sim_in.CNdB + Bfm_dB, sim_out.snr_theorydB,"g;FM Theory;");
+ plot(sim_in.CNdB + Bssb_dB, sim_in.CNdB,"b; SSB Theory;");
+ hold off;
+ grid("minor");
+ xlabel("FM demod input C/No (dB)");
+ ylabel("FM demod output S/N (dB)");
+ legend("boxoff");
+
+endfunction
+
+
+function run_fm_single
+ sim_in.nsam = 96000;
+ sim_in.verbose = 2;
+ sim_in.pre_emp = 0;
+ sim_in.de_emp = 0;
+
+ sim_in.CNdB = 10;
+ sim_out = analog_fm_test(sim_in);
+end
+
+more off;
+
+%run_fm_curves
+%run_fm_single
+
% + do we need interpolation as well?
% + might leave this as pre/de not significant now
% [X] C/No curves?
-% [ ] spectrum plots
+% [X] spectrum plots or analog FM and FSK
+% [ ] figures
rand('state',1);
randn('state',1);
graphics_toolkit ("gnuplot");
+fm;
+
function sim_out = fsk_ber_test(sim_in)
Fs = 96000;
fmark = sim_in.fmark;
fm_states.pre_emp = 0;
fm_states.de_emp = 0;
fm_states.Ts = Ts;
+ fm_states.Fs = Fs;
+ fm_states.fm_max = 3E3;
+ fm_states.fd = 5E3;
fm_states = analog_fm_init(fm_states);
end
for ne = 1:length(EbNodB)
Nerrs = Terrs = Tbits = 0;
- % randn() generates noise across the entire Fs bandwidth, we want to scale
- % the noise power (i.e. the variance) so we get the Eb/No we want in the
- % bandwidth of our FSK signal, which we assume is Rs.
+ % randn() generates noise spread across the entire Fs bandwidth.
+ % The power (aka variance) of this noise is N = NoFs, or No =
+ % N/Fs. The power of each bit is C=1, so the energy of each bit
+ % is Eb=1/Rs. We want to find N as a function of Eb/No, so:
+
+ % Eb/No = (1/Rs)/(N/Fs) = Fs/(RsN)
+ % N = Fs/(Rs(Eb/No))
aEbNodB = EbNodB(ne);
EbNo = 10^(aEbNodB/10);
endfunction
-function fm_states = analog_fm_init(fm_states)
-
- % FM modulator constants
-
- fm_states.Fs = Fs = 96000; FsOn2 = Fs/2;
- fm_states.fm_max = fm_max = 3000; % max modulation freq
- fm_states.fc = 24E3; % carrier frequency
- fm_states.fd = fd = 5E3; % (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.
-
- Ts = fm_states.Ts;
- desired_ncoeffs = 200;
- ncoeffs = floor(desired_ncoeffs/Ts+1)*Ts;
-
- % "coarse" timing offset is half filter length, we have two filters.
- % This is the delay the two filters introduce, so we need to adjust
- % for this when comparing tx to trx bits for BER calcs.
-
- fm_states.nsym_delay = ncoeffs/Ts;
-
- % input filter gets rid of excess noise before demodulator, as too much
- % noise causes atan2() to jump around, e.g. -pi to pi. However this
- % filter can cause harmonic distortion at very high SNRs, as it knocks out
- % some of the FM signal spectra. This filter isn't really required for high
- % SNRs > 20dB.
-
- fc = (Bfm/2)/(FsOn2);
- fm_states.bin = firls(ncoeffs,[0 fc*(1-0.05) fc*(1+0.05) 1],[1 1 0.01 0.01]);
-
- % demoduator output filter to limit us to fm_max (e.g. 3kHz)
-
- fc = (fm_max)/(FsOn2);
- fm_states.bout = firls(ncoeffs,[0 0.95*fc 1.05*fc 1], [1 1 0.01 0.01]);
-
-endfunction
-
-
-function tx = analog_fm_mod(fm_states, mod)
- Fs = fm_states.Fs;
- fc = fm_states.fc; wc = 2*pi*fc/Fs;
- fd = fm_states.fd; wd = 2*pi*fd/Fs;
- nsam = length(mod);
-
- if fm_states.pre_emp
- mod = filter(fm_states.prede,1,mod);
- mod = mod/max(mod); % AGC to set deviation
- end
-
- tx_phase = 0;
- tx = zeros(1,nsam);
-
- for i=0:nsam-1
- w = wc + wd*mod(i+1);
- tx_phase = tx_phase + w;
- tx_phase = tx_phase - floor(tx_phase/(2*pi))*2*pi;
- tx(i+1) = exp(j*tx_phase);
- end
-endfunction
-
-
-function [rx_out rx_bb] = analog_fm_demod(fm_states, rx)
- Fs = fm_states.Fs;
- fc = fm_states.fc; wc = 2*pi*fc/Fs;
- fd = fm_states.fd; wd = 2*pi*fd/Fs;
- nsam = length(rx);
- t = 0:(nsam-1);
-
- rx_bb = rx .* exp(-j*wc*t); % down to complex baseband
- 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.de_emp
- rx_out = filter(1,fm_states.prede,rx_out);
- end
-endfunction
-
-
-function sim_out = analog_fm_test(sim_in)
- nsam = sim_in.nsam;
- CNdB = sim_in.CNdB;
- verbose = sim_in.verbose;
-
- 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 = analog_fm_init(fm_states);
- sim_out.Bfm = fm_states.Bfm;
-
- Fs = fm_states.Fs;
- Bfm = fm_states.Bfm;
- m = fm_states.m; tc = fm_states.tc; fm_max = fm_states.fm_max;
- t = 0:(nsam-1);
-
- fm = 1000; wm = 2*pi*fm/fm_states.Fs;
-
- % start simulation
-
- for ne = 1:length(CNdB)
-
- % work out the variance we need to obtain our C/N in the bandwidth
- % of the FM demod. The gaussian generator randn() generates noise
- % with a bandwidth of Fs
-
- aCNdB = CNdB(ne);
- CN = 10^(aCNdB/10);
- variance = Fs/(CN*Bfm);
-
- % FM Modulator -------------------------------
-
- mod = sin(wm*t);
- tx = analog_fm_mod(fm_states, mod);
-
- % Channel ---------------------------------
-
- noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam));
- rx = tx + noise;
-
- % FM Demodulator
-
- [rx_out rx_bb] = analog_fm_demod(fm_states, rx);
-
- % notch out test tone
-
- w = 2*pi*fm/Fs; beta = 0.99;
- rx_notch = filter([1 -2*cos(w) 1],[1 -2*beta*cos(w) beta*beta], rx_out);
-
- % measure power with and without test tone to determine S+N and N
-
- settle = 1000; % filter settling time, to avoid transients
- nsettle = nsam - settle;
-
- sinad = (rx_out(settle:nsam) * rx_out(settle:nsam)')/nsettle;
- nad = (rx_notch(settle:nsam) * rx_notch(settle:nsam)')/nsettle;
-
- snr = (sinad-nad)/nad;
- sim_out.snrdB(ne) = 10*log10(snr);
-
- % Theory from FMTutorial.pdf, Lawrence Der, Silicon labs paper
-
- snr_theory_dB = aCNdB + 10*log10(3*m*m*(m+1));
- fx = 1/(2*pi*tc); W = fm_max;
- I = (W/fx)^3/(3*((W/fx) - atan(W/fx)));
- I_dB = 10*log10(I);
-
- sim_out.snr_theorydB(ne) = snr_theory_dB;
- sim_out.snr_theory_pre_dedB(ne) = snr_theory_dB + I_dB;
-
- if verbose > 1
- printf("modn index: %2.1f Bfm: %.0f Hz\n", m, Bfm);
- end
-
- if verbose > 0
- printf("C/N: %4.1f SNR: %4.1f dB THEORY: %4.1f dB or with pre/de: %4.1f dB\n",
- aCNdB, 10*log10(snr), snr_theory_dB, snr_theory_dB+I_dB);
- end
-
- if verbose > 1
- figure(1)
- subplot(211)
- plot(20*log10(abs(fft(rx))))
- title('FM Modulator Output Spectrum');
- axis([1 length(tx) 0 100]);
- subplot(212)
- Rx_bb = 20*log10(abs(fft(rx_bb)));
- plot(Rx_bb)
- axis([1 length(tx) 0 100]);
- title('FM Demodulator (baseband) Input Spectrum');
-
- figure(2)
- subplot(211)
- plot(rx_out(settle:nsam))
- axis([1 4000 -1 1])
- subplot(212)
- Rx = 20*log10(abs(fft(rx_out(settle:nsam))));
- plot(Rx(1:10000))
- axis([1 10000 0 100]);
- end
-
- end
-
-endfunction
-
-more off;
-
function run_fsk_curves
sim_in.fmark = 1200;
sim_in.fspace = 2200;
xlabel("Eb/No (dB)");
ylabel("Bit Error Rate (BER)")
- % BER v SNR (3000 Hz noise BW and Eb=C/Rs=1/Rs)
+ % 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(sim_in.Rs/3000);
+ RsOnB_dB = 10*log10(sim_in.Rs/1);
figure(2);
clf;
semilogy(sim_in.EbNodB+RsOnB_dB, fsk_theory.BERvec,'r;FSK theory;')
semilogy(sim_in.EbNodB+RsOnB_dB, fsk_fm_sim.BERvec,'b;FSK over FM sim;')
hold off;
grid("minor");
- axis([min(sim_in.EbNodB) max(sim_in.EbNodB) 1E-4 1])
+ axis([min(sim_in.EbNodB+RsOnB_dB) max(sim_in.EbNodB+RsOnB_dB) 1E-4 1])
legend("boxoff");
- xlabel("S/N for RS=1200 bit/s and 3000 Hz noise bandwidth(dB)");
+ xlabel("C/No for Rs=1200 bit/s and 1 Hz noise bandwidth (dB)");
ylabel("Bit Error Rate (BER)")
end
-function run_fm_curves
- sim_in.nsam = 96000;
- sim_in.verbose = 1;
- sim_in.pre_emp = 0;
- sim_in.de_emp = 0;
- sim_in.CNdB = -4:2:20;
-
- sim_out = analog_fm_test(sim_in);
-
- figure(1)
- clf
- plot(sim_in.CNdB, sim_out.snrdB,"r;FM Simulated;");
- hold on;
- plot(sim_in.CNdB, sim_out.snr_theorydB,"g;FM Theory;");
- plot(sim_in.CNdB, sim_in.CNdB,"b; SSB Theory;");
- hold off;
- grid("minor");
- xlabel("FM demod input C/N (dB)");
- ylabel("FM demod output S/N (dB)");
- legend("boxoff");
-
- % C/No curves
-
- Bfm_dB = 10*log10(sim_out.Bfm);
- Bssb_dB = 10*log10(3000);
-
- figure(2)
- clf
- plot(sim_in.CNdB + Bfm_dB, sim_out.snrdB,"r;FM Simulated;");
- hold on;
- plot(sim_in.CNdB + Bfm_dB, sim_out.snr_theorydB,"g;FM Theory;");
- plot(sim_in.CNdB + Bssb_dB, sim_in.CNdB,"b; SSB Theory;");
- hold off;
- grid("minor");
- xlabel("FM demod input C/No (dB)");
- ylabel("FM demod output S/N (dB)");
- legend("boxoff");
+function run_fsk_single
+ sim_in.fmark = 1200;
+ sim_in.fspace = 2200;
+ sim_in.Rs = 1200;
+ sim_in.nsym = 1200;
+ sim_in.EbNodB = 16;
+ sim_in.fm = 1;
+ sim_in.verbose = 1;
+ fsk_sim = fsk_ber_test(sim_in);
endfunction
-function run_fm_single
- sim_in.nsam = 96000;
- sim_in.verbose = 2;
- sim_in.pre_emp = 0;
- sim_in.de_emp = 0;
-
- sim_in.CNdB = 0;
- sim_out = analog_fm_test(sim_in);
-end
%run_fsk_curves
-%run_fm_curves
-run_fm_single
+run_fsk_single
--- /dev/null
+% gmsk.m
+
+% From: https://github.com/on1arf/gmsk/blob/master/gmskmodem_codec2/API/a_dspstuff.h, which is in turn
+% from Jonathan G4KLX
+
+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,...
+2.731624380156465e-005, 1.218217140199093e-004, 4.711833994209542e-004,...
+1.580581180127418e-003, 4.598383433830095e-003, 1.160259430889949e-002,...
+2.539022692626253e-002, 4.818807833062393e-002, 7.931844341164322e-002,...
+1.132322945270602e-001, 1.401935338024111e-001, 1.505383695578516e-001,...
+1.401935338024111e-001, 1.132322945270601e-001, 7.931844341164328e-002,...
+4.818807833062393e-002, 2.539022692626253e-002, 1.160259430889949e-002,...
+4.598383433830090e-003, 1.580581180127420e-003, 4.711833994209542e-004,...
+1.218217140199093e-004, 2.731624380156465e-005, 5.312253663302753e-006,...
+8.959797186410563e-007, 1.310626978701910e-007, 1.662729407135958e-008,...
+1.829471305298363e-009, 1.745786683011426e-010, 1.444835156335356e-011,...
+1.037067381285011e-012, 6.455906007234699e-014];
+
+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.000503471198655, -0.001079755887508, -0.001671728086040, -0.002205032425392,...
+-0.002594597675000, -0.002754194565297, -0.002608210441859, -0.002104352817854,...
+-0.001225654870420, 0.000000000000000, 0.001494548041184, 0.003130012785731,...
+0.004735238379172, 0.006109242742194, 0.007040527007323, 0.007330850462455,...
+0.006821247169795, 0.005417521811131, 0.003112202160626, -0.000000000000000,...
+-0.003715739376345, -0.007727358782391, -0.011638713107503, -0.014992029537478,...
+-0.017304097563429, -0.018108937286588, -0.017003180218569, -0.013689829477969,...
+-0.008015928769710, 0.000000000000000, 0.010154104792614, 0.022059114281395,...
+0.035162729807337, 0.048781621388364, 0.062148583345584, 0.074469032280094,...
+0.084982001723750, 0.093020219991183, 0.098063819576269, 0.099782731268437,...
+0.098063819576269, 0.093020219991183, 0.084982001723750, 0.074469032280094,...
+0.062148583345584, 0.048781621388364, 0.035162729807337, 0.022059114281395,...
+0.010154104792614, 0.000000000000000, -0.008015928769710, -0.013689829477969,...
+-0.017003180218569, -0.018108937286588, -0.017304097563429, -0.014992029537478,...
+-0.011638713107503, -0.007727358782391, -0.003715739376345, -0.000000000000000,...
+0.003112202160626, 0.005417521811131, 0.006821247169795, 0.007330850462455,...
+0.007040527007323, 0.006109242742194, 0.004735238379172, 0.003130012785731,...
+0.001494548041184, 0.000000000000000, -0.001225654870420, -0.002104352817854,...
+-0.002608210441859, -0.002754194565297, -0.002594597675000, -0.002205032425392,...
+-0.001671728086040, -0.001079755887508, -0.000503471198655, -0.000000000000000,...
+0.000393913058926, 0.000661865814384, 0.000805143268199, 0.000838293462576,...
+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
+ % work out delays of filter to align bits
+ % plot eye diagrams, BERcurves, theoretical results, spectrum - will it pass thru a HT?
+endfunction