moved analog FM functions into a file that can be used elsewhere
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 15 Dec 2014 05:59:40 +0000 (05:59 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 15 Dec 2014 05:59:40 +0000 (05:59 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@1970 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/fm.m [new file with mode: 0644]
codec2-dev/octave/fsk.m
codec2-dev/octave/gmsk.m [new file with mode: 0644]

diff --git a/codec2-dev/octave/fm.m b/codec2-dev/octave/fm.m
new file mode 100644 (file)
index 0000000..0896dff
--- /dev/null
@@ -0,0 +1,257 @@
+% 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
+
index 5c9b6be5db0829b861f08ffd6c0f59a0d1715c6f..a321fc3f1fc0da71b26c42fd0ab7d469db60b423 100644 (file)
 %       + 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;
@@ -46,6 +49,9 @@ function sim_out = fsk_ber_test(sim_in)
     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
 
@@ -54,9 +60,13 @@ function sim_out = fsk_ber_test(sim_in)
   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);
@@ -174,198 +184,6 @@ function sim_out = fsk_ber_test(sim_in)
 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;
@@ -397,11 +215,11 @@ function run_fsk_curves
   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;')
@@ -410,63 +228,25 @@ function run_fsk_curves
   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
 
diff --git a/codec2-dev/octave/gmsk.m b/codec2-dev/octave/gmsk.m
new file mode 100644 (file)
index 0000000..a873313
--- /dev/null
@@ -0,0 +1,55 @@
+% 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