diff options
| author | drowe67 <[email protected]> | 2023-07-07 13:55:18 +0930 |
|---|---|---|
| committer | David Rowe <[email protected]> | 2023-07-07 13:55:18 +0930 |
| commit | 87025ac7d2a23192bf124ce8684626f821b099ea (patch) | |
| tree | 6946606f4eb57c85f9b692d6cd4901e189cdd619 /octave | |
| parent | 70c5595be1e5e0250a4bc5fc89f0adce861b4765 (diff) | |
rm-ed a bunch of octave files
Diffstat (limited to 'octave')
| -rw-r--r-- | octave/2400ab_frame_design.ods | bin | 19194 -> 0 bytes | |||
| -rw-r--r-- | octave/apsk_ser.m | 56 | ||||
| -rw-r--r-- | octave/c2wideband_map | 8 | ||||
| -rw-r--r-- | octave/closed_quant_slope.m | 3 | ||||
| -rw-r--r-- | octave/cma.m | 114 | ||||
| -rw-r--r-- | octave/codec2_demo.m | 108 | ||||
| -rw-r--r-- | octave/fm.m | 484 | ||||
| -rw-r--r-- | octave/fm_radio_filt_model.txt | 8 | ||||
| -rw-r--r-- | octave/fmfsk.m | 346 | ||||
| -rw-r--r-- | octave/fsk4_dmr.m | 540 | ||||
| -rw-r--r-- | octave/fsk_basic.m | 60 | ||||
| -rw-r--r-- | octave/fsk_cml.m | 132 | ||||
| -rw-r--r-- | octave/fsk_cml_sam.m | 376 | ||||
| -rwxr-xr-x | octave/fsk_demod_BER_test.py | 610 | ||||
| -rw-r--r-- | octave/fsk_llr_plot.m | 51 | ||||
| -rw-r--r-- | octave/fsk_llr_test.m | 294 | ||||
| -rw-r--r-- | octave/fsk_lock_down.m | 393 | ||||
| -rw-r--r-- | octave/fsk_v_afsk.m | 232 | ||||
| -rw-r--r-- | octave/fskdemodgui.py | 220 | ||||
| -rw-r--r-- | octave/gmsk.m | 1021 | ||||
| -rw-r--r-- | octave/hp_filt.m | 12 |
21 files changed, 0 insertions, 5068 deletions
diff --git a/octave/2400ab_frame_design.ods b/octave/2400ab_frame_design.ods Binary files differdeleted file mode 100644 index 015fc6b..0000000 --- a/octave/2400ab_frame_design.ods +++ /dev/null diff --git a/octave/apsk_ser.m b/octave/apsk_ser.m deleted file mode 100644 index c357a64..0000000 --- a/octave/apsk_ser.m +++ /dev/null @@ -1,56 +0,0 @@ -function [svec] = apsk_ser(Esvec, M) -% Return approx (uncoded) SER for 16 or 32 APSK, or plot SERs if no args -% USAGE apsk_ser(); OR SERvec = apsk_ser(Esvec_dBs, M) -% VK5DSP, July 2020 - -% Fit 3rd order polys to data in Figure 3 of -% "Perf anal of APSK mod for DVB-S2 trans over nonlinear channels" -% by Wonjin Sung et al; -% Published in Int. J. Satellite Communications Networking, 2009 -% (see Fig 1 for constellation diagrams and bit mappings) - -% from the Figure data, we get 2 polynomials ... -%{ -Fig3_Es = [10 12.5 15 17.5 20 21]; -Fig3_16_SER = [2e-1 8e-2 2.5e-2 3e-3 1e-4 2e-5]; -Fig3_32_SER = [4e-1 2.8e-1 1.5e-1 4e-2 1e-2 4e-3 1.05e-5]; - -p16 = polyfit(Fig3_Es, log10(Fig3_16_SER), 3) -p32 = polyfit([Fig3_Es 25], log10(Fig3_32_SER), 3) -%} - -p16 = [-0.001816729857582 0.053322736530042 -0.653699203208837 2.315143765468362]; -p32 = [-0.001418894437779 0.049903575028137 -0.669130206501662 2.743251206655785]; - -if nargin==0 - EStest = 10:21; - EStest2 = 10:25; - SER16 = polyval(p16, EStest); - SER32 = polyval(p32, EStest2); - figure(55); - semilogy(EStest, 10.^SER16, 'b', EStest2, 10.^SER32, 'g') - title('Approx Symb Error Rates for APSK') - xlabel('Es/N0 (dB)') - ylabel('SER'); grid on; - legend('16APSK','32APSK'); legend('boxoff') -elseif nargin~=2, - error('usage is apsk_ser(Esvec, M)'); - -else - if (M~=16)&&(M~=32) - error('M must be 16 or 32') - end - if min(Esvec)<8, - error('Es/No values should be > 8 dB'); - end - if (M==16) && (max(Esvec)>23), - error('Es/No values should be < 23 dB'); - end - if (M==32) && (max(Esvec)>27), - error('Es/No values should be < 27 dB'); - end - - if M==16, svec = 10.^polyval(p16, Esvec); end - if M==32, svec = 10.^polyval(p32, Esvec); end - -end diff --git a/octave/c2wideband_map b/octave/c2wideband_map deleted file mode 100644 index 946aa19..0000000 --- a/octave/c2wideband_map +++ /dev/null @@ -1,8 +0,0 @@ - 1.00000000e+00 3.00000000e+00 7.00000000e+00 8.00000000e+00 1.10000000e+01 1.50000000e+01 1.90000000e+01 2.60000000e+01 2.10000000e+01 2.40000000e+01 2.00000000e+01 3.00000000e+01 3.80000000e+01 4.80000000e+01 2.90000000e+01 3.20000000e+01 4.20000000e+01 6.40000000e+01 6.30000000e+01 5.40000000e+01 5.60000000e+01 5.00000000e+01 7.20000000e+01 9.10000000e+01 7.80000000e+01 6.70000000e+01 5.70000000e+01 7.00000000e+01 7.60000000e+01 1.04000000e+02 - 2.00000000e+00 5.00000000e+00 1.20000000e+01 2.20000000e+01 2.30000000e+01 3.30000000e+01 3.10000000e+01 4.30000000e+01 3.40000000e+01 3.90000000e+01 4.40000000e+01 4.10000000e+01 4.90000000e+01 6.00000000e+01 5.10000000e+01 7.70000000e+01 9.70000000e+01 9.00000000e+01 1.14000000e+02 8.10000000e+01 1.21000000e+02 1.13000000e+02 8.00000000e+01 9.90000000e+01 1.18000000e+02 1.09000000e+02 8.90000000e+01 1.02000000e+02 1.07000000e+02 1.19000000e+02 - 4.00000000e+00 9.00000000e+00 1.60000000e+01 2.70000000e+01 4.00000000e+01 4.50000000e+01 4.70000000e+01 5.50000000e+01 4.60000000e+01 6.10000000e+01 6.50000000e+01 6.60000000e+01 7.50000000e+01 8.30000000e+01 6.80000000e+01 1.06000000e+02 9.40000000e+01 1.27000000e+02 1.36000000e+02 1.38000000e+02 1.37000000e+02 1.41000000e+02 8.80000000e+01 1.50000000e+02 1.59000000e+02 1.47000000e+02 1.33000000e+02 1.28000000e+02 1.90000000e+02 1.31000000e+02 - 6.00000000e+00 1.80000000e+01 2.80000000e+01 5.20000000e+01 8.40000000e+01 6.20000000e+01 7.30000000e+01 7.90000000e+01 7.10000000e+01 9.80000000e+01 8.50000000e+01 9.30000000e+01 8.70000000e+01 1.12000000e+02 1.15000000e+02 1.10000000e+02 1.29000000e+02 1.99000000e+02 1.42000000e+02 1.70000000e+02 2.18000000e+02 1.62000000e+02 1.74000000e+02 2.05000000e+02 1.43000000e+02 1.64000000e+02 1.72000000e+02 1.75000000e+02 2.02000000e+02 2.16000000e+02 - 1.00000000e+01 2.50000000e+01 3.50000000e+01 6.90000000e+01 7.40000000e+01 9.20000000e+01 9.60000000e+01 1.35000000e+02 1.11000000e+02 1.56000000e+02 9.50000000e+01 1.05000000e+02 1.26000000e+02 1.34000000e+02 1.63000000e+02 2.15000000e+02 1.49000000e+02 1.48000000e+02 1.53000000e+02 2.28000000e+02 1.80000000e+02 2.30000000e+02 2.11000000e+02 2.13000000e+02 2.40000000e+02 2.24000000e+02 2.09000000e+02 2.22000000e+02 2.26000000e+02 1.73000000e+02 - 1.30000000e+01 3.70000000e+01 5.30000000e+01 1.08000000e+02 8.60000000e+01 1.17000000e+02 1.03000000e+02 1.00000000e+02 1.24000000e+02 1.22000000e+02 2.01000000e+02 1.92000000e+02 2.03000000e+02 2.19000000e+02 1.91000000e+02 1.69000000e+02 2.20000000e+02 2.14000000e+02 1.39000000e+02 2.04000000e+02 2.36000000e+02 1.79000000e+02 1.81000000e+02 2.00000000e+02 1.93000000e+02 2.21000000e+02 2.35000000e+02 1.87000000e+02 2.08000000e+02 2.17000000e+02 - 1.40000000e+01 3.60000000e+01 5.90000000e+01 1.16000000e+02 1.30000000e+02 1.78000000e+02 1.40000000e+02 1.32000000e+02 1.57000000e+02 1.52000000e+02 1.68000000e+02 1.60000000e+02 1.23000000e+02 1.97000000e+02 1.83000000e+02 2.34000000e+02 1.44000000e+02 2.23000000e+02 1.65000000e+02 1.67000000e+02 2.31000000e+02 1.96000000e+02 1.71000000e+02 2.33000000e+02 1.82000000e+02 2.25000000e+02 1.88000000e+02 2.37000000e+02 2.27000000e+02 2.10000000e+02 - 1.70000000e+01 5.80000000e+01 8.20000000e+01 1.01000000e+02 1.20000000e+02 1.86000000e+02 1.46000000e+02 1.25000000e+02 1.51000000e+02 1.94000000e+02 1.61000000e+02 1.84000000e+02 1.58000000e+02 1.89000000e+02 1.77000000e+02 1.95000000e+02 1.45000000e+02 1.66000000e+02 2.06000000e+02 1.54000000e+02 1.98000000e+02 2.12000000e+02 2.38000000e+02 2.29000000e+02 1.85000000e+02 1.76000000e+02 2.39000000e+02 1.55000000e+02 2.07000000e+02 2.32000000e+02 diff --git a/octave/closed_quant_slope.m b/octave/closed_quant_slope.m deleted file mode 100644 index 7d5a465..0000000 --- a/octave/closed_quant_slope.m +++ /dev/null @@ -1,3 +0,0 @@ -function b = closed_quant_slope(b) - b(1) = max(0.5, b(1)); -end diff --git a/octave/cma.m b/octave/cma.m deleted file mode 100644 index a5a2195..0000000 --- a/octave/cma.m +++ /dev/null @@ -1,114 +0,0 @@ -% cma.m -% -% Constant modulus equaliser example from: -% -% http://dsp.stackexchange.com/questions/23540/matlab-proper-estimation-of-weights-and-how-to-calculate-mse-for-qpsk-signal-f -% -% Adapted to run bpsk and fsk signals - - rand('seed',1); - randn('seed',1); - - N = 5000; % # symbols - h = [1 0 0 0 0 0 0.0 0.5]; % simulation of HF multipath channel impulse response - h = h/norm(h); - Le = 20; % equalizer length - mu = 1E-3; % step size - snr = 30; % snr in dB - M = 10; % oversample rate, e.g. Rs=400Hz at Fs=8000Hz - - tx_type = "fsk"; % select modulation type here "bpsk" or "fsk" - - if strcmp(tx_type, "bpsk") - s0 = round( rand(N,1) )*2 - 1; % BPSK signal - s0M = zeros(N*M,1); % oversampled BPSK signal - k = 1; - for i=1:M:N*M - s0M(i:i+M-1) = s0(k); - k ++; - end - end - - if strcmp(tx_type, "fsk") - tx_bits = round(rand(1,N)); - - % continuous phase FSK modulator - - w1 = pi/4; - w2 = pi/2; - tx_phase = 0; - tx = zeros(M*N,1); - - for i=1:N - for k=1:M - if tx_bits(i) - tx_phase += w2; - else - tx_phase += w1; - end - tx((i-1)*M+k) = exp(j*tx_phase); - end - end - - s0M = tx; - end - - s = filter(h,1,s0M); % filtered signal - - % add Gaussian noise at desired snr - - n = randn(N*M,1); - vs = var(s); - vn = vs*10^(-snr/10); - n = sqrt(vn)*n; - r = s + n; % received signal - - e = zeros(N*M,1); % error - w = zeros(Le,1); % equalizer coefficients - w(Le)=1; % actual filter taps are flipud(w)! - - yd = zeros(N*M,1); - - for i = 1:N*M-Le, - x = r(i:Le+i-1); - y = w'*x; - yd(i)=y; - e(i) = abs(y).^2 - 1; - w = w - mu * e(i) * real(conj(y) * x); - end - - np = 100; % # sybmols to plot (last np will be plotted); np < N! - - figure(1); clf; - %subplot(211), plot( 1:np, e(N-np+1-Le+1:N-Le+1).*e(N-np+1-Le+1:N-Le+1)), title('error') - subplot(211), plot(e.*e), title('error'); - subplot(212), stem(conv(flipud(w),h)), title('equalized channel impulse response') - - figure(2); clf; - subplot(311) - plot(1:np, s0M(N-np+1:N)) - title('transmitted, received, and equalized signal') - subplot(312) - plot(1:np, r(N-np+1:N)) - subplot(313) - plot(1:np, yd(N-np+1-Le+1:N-Le+1)) - - figure(3); clf; - h1 = freqz(h); - h2 = freqz(flipud(w)); - h3 = freqz(conv(flipud(w),h)); - subplot(311); plot(20*log10(abs(h1))); - title('channel, equaliser, combined freq resp') - subplot(312); plot(20*log10(abs(h2))); - subplot(313); plot(20*log10(abs(h3))); - - figure(4); - subplot(211) - plot(20*log10(abs(fft(s0M)))) - axis([1 length(s0M) 0 80]); - grid; - subplot(212) - plot(20*log10(abs(fft(s)))) - axis([1 length(s0M) 0 80]); - grid; - diff --git a/octave/codec2_demo.m b/octave/codec2_demo.m deleted file mode 100644 index 0f3950b..0000000 --- a/octave/codec2_demo.m +++ /dev/null @@ -1,108 +0,0 @@ -% Copyright David Rowe 2012 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% codec2_demo.m - -% Designed as an educational tool to explain the operation of Codec 2 -% for conference and user group presentations on a projector. An -% alternative to static overhead slides. -% -% Derived from codec2-dev/octave/plamp.m -% -% usage: -% octave:1> plamp("../src/hts2a",40) -% -% Then press: -% c - to cycle through the wavform being displayed on the figure -% n - next frame -% b - back one frame -% -% tip: hold down n or b to animate the display -% -% The text files used as input are generated using c2sim: -% -% /codec2-dev/src$ c2sim ../raw/hts2a.raw --dump hts2a -% -% The Codec 2 README explains how to build c2sim with dump files -% enabled. - -function codec2_demo(samname, f) - - sn_name = strcat(samname,"_sn.txt"); - Sn = load(sn_name); - - sw_name = strcat(samname,"_sw.txt"); - Sw = load(sw_name); - - model_name = strcat(samname,"_model.txt"); - model = load(model_name); - - figure(1); - - k = ' '; - wf = "Sn"; - do - - if strcmp(wf,"Sn") - clf; - s = [ Sn(2*f-1,:) Sn(2*f,:) ]; - plot(s); - axis([1 length(s) -20000 20000]); - end - - if (strcmp(wf,"Sw")) - clf; - plot((0:255)*4000/256, Sw(f,:),";Sw;"); - end - - if strcmp(wf,"SwAm") - Wo = model(f,1); - L = model(f,2); - Am = model(f,3:(L+2)); - plot((0:255)*4000/256, Sw(f,:),";Sw;"); - hold on; - plot((1:L)*Wo*4000/pi, 20*log10(Am),"+;Am;r"); - axis([1 4000 -10 80]); - hold off; - end - - if strcmp(wf,"Am") - Wo = model(f,1); - L = model(f,2); - Am = model(f,3:(L+2)); - plot((1:L)*Wo*4000/pi, 20*log10(Am),"+;Am;r"); - axis([1 4000 -10 80]); - end - - % interactive menu - - printf("\rframe: %d menu: n-next b-back w-cycle window q-quit", f); - fflush(stdout); - k = kbhit(); - if (k == 'n') - f = f + 1; - end - if (k == 'b') - f = f - 1; - end - if (k == 'w') - if strcmp(wf,"Sn") - next_wf = "Sw"; - end - if strcmp(wf,"Sw") - next_wf = "SwAm"; - end - if strcmp(wf,"SwAm") - next_wf = "Am"; - end - if strcmp(wf,"Am") - next_wf = "Sn"; - end - wf = next_wf; - end - - until (k == 'q') - printf("\n"); - -endfunction diff --git a/octave/fm.m b/octave/fm.m deleted file mode 100644 index e25432f..0000000 --- a/octave/fm.m +++ /dev/null @@ -1,484 +0,0 @@ -% fm.m -% David Rowe Dec 2014 -% -% Analog FM Octave simulation functions. - -1; - -graphics_toolkit ("gnuplot"); - -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 - 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 - fm_states.ph_dont_limit = 0; % Limit rx delta-phase - - % 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 fm_fir_coeff_file(fm_states, filename) - global gt_alpha5_root; - global Nfilter; - - f=fopen(filename,"wt"); - - fprintf(f,"/* Generated by fm_fir_coeff_file() Octave function in fm.m */\n\n"); - fprintf(f,"const float bin[]={\n"); - for m=1:length(fm_states.bin)-1 - fprintf(f," %g,\n", fm_states.bin(m)); - endfor - fprintf(f," %g\n};\n\n", fm_states.bin(length(fm_states.bin))); - - fprintf(f,"const float bout[]={\n"); - for m=1:length(fm_states.bout)-1 - fprintf(f," %g,\n", fm_states.bout(m)); - endfor - fprintf(f," %g\n};\n", fm_states.bout(length(fm_states.bout))); - - fclose(f); -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); - - % differentiate first, in rect domain, then find angle, this puts - % signal on the positive side of the real axis - - rx_bb_diff = [ 1 rx_bb(2:nsam) .* conj(rx_bb(1:nsam-1))]; - rx_out = atan2(imag(rx_bb_diff),real(rx_bb_diff)); - - % limit maximum phase jumps, to remove static type noise at low SNRs - if !fm_states.ph_dont_limit - rx_out(find(rx_out > wd)) = wd; - rx_out(find(rx_out < -wd)) = -wd; - end - rx_out *= (1/wd); - - 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 -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.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; - - 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 = 20; - sim_out = analog_fm_test(sim_in); -end - - -function fm_mod_file(file_name_out, file_name_in, CNdB) - fm_states.Fs = 48000; - fm_states.fm_max = 3E3; - fm_states.fd = 5E3; - fm_states.fc = fm_states.Fs/4; - fm_states.pre_emp = 0; - fm_states.de_emp = 0; - fm_states.Ts = 1; - fm_states.output_filter = 1; - fm_states = analog_fm_init(fm_states); - - if nargin == 1 - nsam = fm_states.Fs * 10; - t = 0:(nsam-1); - fm = 1000; wm = 2*pi*fm/fm_states.Fs; - mod = sin(wm*t); - else - fin = fopen(file_name_in,"rb"); - mod = fread(fin,"short")'; - mod /= 32767; - fclose(fin); - end - tx = analog_fm_mod(fm_states, mod); - - if (nargin == 3) - % Optionally add some noise - - CN = 10^(CNdB/10); - variance = fm_states.Fs/(CN*fm_states.Bfm); - tx += sqrt(variance)*randn(1,length(tx)); - end - - tx_out = tx*16384; - fout = fopen(file_name_out,"wb"); - fwrite(fout, tx_out, "short"); - fclose(fout); -endfunction - - -function fm_demod_file(file_name_out, file_name_in) - fin = fopen(file_name_in,"rb"); - rx = fread(fin,"short")'; - rx = rx(100000:length(rx)); % strip of wave header - fclose(fin); - - Fs = fm_states.Fs = 48000; - fm_max = fm_states.fm_max = 3E3; - fd = fm_states.fd = 5E3; - fm_states.fc = 12E3; - - fm_states.pre_emp = 0; - fm_states.de_emp = 1; - fm_states.Ts = 1; - fm_states.output_filter = 1; - fm_states = analog_fm_init(fm_states); - - [rx_out rx_bb] = analog_fm_demod(fm_states, rx); - - rx_out *= 20000; - fout = fopen(file_name_out,"wb"); - fwrite(fout, rx_out, "short"); - fclose(fout); - - figure(1) - subplot(211) - plot(rx) - subplot(212) - plot(20*log10(abs(fft(rx)))) - title('FM Dmodulator Input Spectrum'); - - figure(2) - subplot(211) - Rx_bb = 20*log10(abs(fft(rx_bb))); - plot(Rx_bb) - title('FM Demodulator (baseband) Input Spectrum'); - - subplot(212) - plot(20*log10(abs(fft(rx_out)))) - title('FM Dmodulator Output Spectrum'); - - figure(3) - plot(rx_out) - title('FM Dmodulator Output'); - - % estimate SNR, C/No etc - - npower_window = 1024; - rx_power = conv(rx.^2,ones(1,npower_window))/(npower_window); - rx_power_dB = 10*log10(rx_power); - figure; - subplot(211) - plot(rx); - subplot(212) - plot(rx_power_dB); - axis([1 length(rx_power) max(rx_power_dB)-9 max(rx_power_dB)+1]) - grid("minor") - - % estimate FM demod output SNR if a 1000 Hz tone is present - - w = 2*pi*1000/Fs; beta = 0.99; - rx_notch = filter([1 -2*cos(w) 1],[1 -2*beta*cos(w) beta*beta], rx_out); - - rx_out_power = conv(rx_out.^2,ones(1,npower_window))/(npower_window); - rx_out_power_dB = 10*log10(rx_out_power); - rx_notch_power = conv(rx_notch.^2,ones(1,npower_window))/(npower_window); - rx_notch_power_dB = 10*log10(rx_notch_power); - figure; - plot(rx_out_power_dB,'r;FM demod output power;'); - hold on; - plot(rx_notch_power_dB,'b;1000 Hz notch filter output power;'); - plot(rx_out_power_dB-rx_notch_power_dB,'g;1000 Hz tone SNR;'); - hold off; - legend("boxoff"); - ylabel('dB'); - xlabel('Time (samples)'); - grid("minor") - -endfunction - - -% generate filter coeffs for C implementation of FM demod - -function make_coeff_file - fm_states.Fs = 44400; - fm_states.fm_max = 3E3; - fm_states.fd = 5E3; - fm_states.fc = fm_states.Fs/4; - - fm_states.pre_emp = 0; - fm_states.de_emp = 0; - fm_states.Ts = 1; - fm_states.output_filter = 1; - fm_states = analog_fm_init(fm_states); - - fm_fir_coeff_file(fm_states, "fm_fir_coeff.h") -endfunction - -function test_fm_modulator - fm_states.Fs = 48000; - fm_states.fm_max = 3E3; - fm_states.fd = 5E3; - %fm_states.fc = fm_states.Fs/4; - fm_states.fc = 0; - - fm_states.pre_emp = 0; - fm_states.de_emp = 0; - fm_states.Ts = 1; - fm_states.output_filter = 1; - fm_states = analog_fm_init(fm_states); - - test_t = [1:(fm_states.Fs*10)]; - test_freq1 = 2*pi*3000/fm_states.Fs; - test_freq2 = 2*pi*1000/fm_states.Fs; - - test_sig = .5*sin(test_t*test_freq1) + .5*sin(test_t*test_freq2); - %test_sig = zeros(1,length(test_t)); - %test_sig = ones(1,length(test_t)); - - ftsig = fopen("fm_test_sig.raw","wb"); - fwrite(ftsig,test_sig*16384,"short"); - fclose(ftsig); - - system("../fm_test fm_test_sig.raw fm_test_out.raw"); - ftmod = fopen("fm_test_out.raw","r"); - test_mod_p = rot90(fread(ftmod,"short"))/16384; - test_mod_r = test_mod_p(1:2:length(test_mod_p)); - test_mod_i = test_mod_p(2:2:length(test_mod_p)); - test_mod = test_mod_r .+ i*test_mod_i; - fclose(ftmod); - - comp_mod = analog_fm_mod(fm_states,test_sig); - - figure(1) - comp_mod_real = real(comp_mod); - size(comp_mod_real) - size(test_mod) - mod_diff = zeros(1,length(test_mod)); - mod_diff = test_mod .- comp_mod; - plot(test_t,real(test_mod .- comp_mod),test_t,imag(test_mod .- comp_mod)); - -endfunction - -more off; - -%run_fm_curves -%fm_demod_file("ssb_fm_out.raw","~/Desktop/ssb_fm.wav") -%fm_demod_file("ssb25_fm_de.raw", "~/Desktop/ssb25db.wav") -%run_fm_single -%make_coeff_file -%fm_mod_file("fm_1000.raw"); -%test_fm_modulator diff --git a/octave/fm_radio_filt_model.txt b/octave/fm_radio_filt_model.txt deleted file mode 100644 index 368f7e2..0000000 --- a/octave/fm_radio_filt_model.txt +++ /dev/null @@ -1,8 +0,0 @@ -# Created by Octave 4.0.0, Wed Feb 10 20:14:16 2016 CST <baobrien@baobrien-d-desktop> -# name: filt -# type: matrix -# rows: 1 -# columns: 1001 - 4.111934884608988e-05 4.529266752036107e-05 5.012294778018435e-05 5.387673663540811e-05 5.656595617636202e-05 5.754247327156636e-05 5.745029328387444e-05 5.810382275274602e-05 6.038440324493526e-05 6.151519405194188e-05 5.995001465356246e-05 5.965124091807071e-05 5.894036887788694e-05 5.282714911187688e-05 4.16286613756742e-05 2.733281113424976e-05 1.151301392357209e-05 -5.011952829637155e-06 -2.038060526743049e-05 -2.88950898384027e-05 -3.140205230802841e-05 -2.680377848432233e-05 -1.425766597021347e-05 2.754137907668088e-06 2.533067638032289e-05 5.098405707491412e-05 7.524935743597888e-05 9.490108720790089e-05 0.0001114408649087926 0.0001232907290105822 0.0001298558067129026 0.0001314146843973354 0.000126874153326206 0.0001211256678148737 0.0001128586850331393 0.0001021842594143064 9.340521285980458e-05 8.465588232787989e-05 7.602754602719386e-05 7.476083494605245e-05 7.795597606716329e-05 8.085052742218461e-05 8.855370570715267e-05 9.993747197784049e-05 0.0001100627878797455 0.0001157727179778551 0.0001211159686312227 0.0001276193402990368 0.0001274951568358215 0.0001237937282219144 0.0001199864757482567 0.0001145773485240607 0.0001078643622955885 9.949247921871633e-05 9.282321723585223e-05 8.670764514445847e-05 8.236867203976522e-05 8.178131856442866e-05 8.653490906422325e-05 9.474089964535933e-05 0.0001046664881685255 0.0001171407871684955 0.0001284375906324497 0.0001435371988464706 0.0001538019728144609 0.000157821986861445 0.0001617032194688335 0.000161139381508245 0.0001573044997290387 0.0001505912434581831 0.0001444545226667355 0.0001371895926370354 0.0001282916901414378 0.0001164777490003531 0.0001017868860228118 8.38857960955751e-05 6.397168698635467e-05 4.505424718694719e-05 2.609404178225459e-05 9.390120446350423e-06 -6.139147851498066e-06 -1.74535589572342e-05 -2.432568386408825e-05 -2.719609151945925e-05 -2.559556014458576e-05 -1.629614429158625e-05 -1.3631941893713e-06 1.383091083546116e-05 3.467670361993851e-05 6.180780852724694e-05 8.838658249276322e-05 0.0001119097915371812 0.0001372537482217498 0.0001573246833069575 0.0001685794934648446 0.0001736521706833196 0.0001743894362264053 0.000170440134727034 0.0001570034141396783 0.0001415291461137907 0.0001293030383250257 0.0001170129681552746 0.0001056747747162765 9.587608125739278e-05 8.209630166614193e-05 6.930393501059805e-05 6.085367944553192e-05 5.275816447026444e-05 4.194241650880381e-05 3.420239452264183e-05 3.304591574725996e-05 3.024288241806293e-05 3.012154305643325e-05 3.738917134601837e-05 4.770388200482866e-05 5.781977436721907e-05 7.420695272211161e-05 8.723833080263831e-05 9.505250826889752e-05 0.0001065158066456628 0.0001104312822684124 0.0001098275912662366 0.000107206024810334 0.0001025516069602837 9.894340278646091e-05 9.407503918864603e-05 8.7319365554045e-05 8.141197975598374e-05 8.422482885728253e-05 9.626406072217444e-05 0.0001047837385634293 0.0001197702730190546 0.0001441800101209614 0.000160456934552367 0.0001774138786603596 0.0001926926553236415 0.0001984206468214209 0.0001945597427125206 0.0001836754090288533 0.0001644147124457892 0.0001450923278130216 0.000131637876080013 0.0001194268883117051 0.0001142849115928036 0.0001090067194513346 0.0001167532618894625 0.0001293962113881341 0.0001304604164937581 0.0001310955784948294 0.0001353721674033472 0.0001341681487126995 0.0001198264422032539 0.0001012747411089002 7.153140827090248e-05 4.191084339680094e-05 2.303431628630918e-05 1.138159045954137e-05 1.214910078795732e-05 2.029103100753587e-05 3.914274819454478e-05 7.436668389413293e-05 0.0001131985409445768 0.0001541458098186209 0.0001981874567952984 0.0002293068332902554 0.0002574810386944711 0.0002797524248574673 0.000280098997160137 0.0002691199306180641 0.0002483396512410882 0.0002130987320308693 0.0001696438030090705 0.0001233684249499195 8.43926061428125e-05 5.198681974621152e-05 3.567320177887548e-05 4.949705291163776e-05 7.787515033943651e-05 0.0001300149878038 0.0002020713027294515 0.0002783857988869775 0.0003565185832637648 0.0004241895316775589 0.000479626912439881 0.0005159892932585141 0.0005326203262494115 0.0005366426137323655 0.0005240387399286967 0.0005041051313799223 0.0004864550530904512 0.0004727638664349475 0.0004645950604270155 0.0004548102499574072 0.000455527247961606 0.0004612992161328969 0.0004674556394527062 0.0004716258022734013 0.0004668087315020519 0.0004622641649603435 0.0004465585690774155 0.000425258701401412 0.00039751658715728 0.0003635800988525798 0.0003328960310700836 0.0003048849901663823 0.0002810007495450609 0.0002673721567286335 0.0002702604512501591 0.0002862572870604029 0.0003084406671027542 0.0003355248671444351 0.0003761847353904122 0.0004148377173810094 0.0004485558688712593 0.0004800892227209435 0.0005013700140284393 0.0005138034260364562 0.000523675449512398 0.0005317635475002124 0.0005322207128124348 0.0005292033760978139 0.0005269773541343302 0.0005175977188382754 0.0004993895797658196 0.0004858979454010794 0.0004684001252653537 0.0004587034220581103 0.0004638037190751786 0.0004761823498486873 0.0004885372867673723 0.0005043751480265481 0.0005446800229664075 0.0005890158012063288 0.0006231865628459465 0.0006616542116343469 0.0006972231165628774 0.0007040149509799145 0.000712178888077443 0.0007101214361108593 0.0006931453291674526 0.0006924927262763824 0.0006853703272953059 0.0007004172297637321 0.0007466476213180197 0.00080183074513318 0.0008653125079193471 0.000932413333045206 0.0009780458869915588 0.001000137248327646 0.0009991649016365171 0.0009565650123854935 0.0008785166874477056 0.0007712153995398795 0.0006440612144425213 0.0005014886072902373 0.0003647207877911045 0.0002314452035966864 0.0001215619643926978 5.102070582077005e-05 3.53881850929961e-06 5.024961716295417e-06 4.952203678213442e-05 0.0001260092656341149 0.0002274300947131321 0.0003513571295408856 0.0004918678731923972 0.0006212972022330388 0.0007279470267812717 0.0008134213275591345 0.000856084224659261 0.0008367582390177951 0.0007963231371937705 0.0007308920792511363 0.0006278320617335954 0.0005155411192978977 0.0004049366284488046 0.0003146483534597321 0.000232809201512697 0.0001611481391302676 0.0001276356087955276 0.0001258477694353633 0.000148410784426943 0.0001791264427517458 0.0002152831436118296 0.0002722059220446592 0.0003148588847685937 0.0003235563379009122 0.0003234961169067191 0.0003141523386480817 0.0002970986813756301 0.0002753731712067278 0.0002391073092942144 0.0002217724458095815 0.0002500563490184254 0.0002860821491880679 0.0003319979542964821 0.0004127132318400748 0.0005101781526047186 0.0006116301200159012 0.0007154307983328388 0.0008127224793692451 0.0008838704404048489 0.0009291650300966136 0.0009480193252716987 0.0009397344198433433 0.0009095689814133252 0.0008547924134192681 0.0007841007737282723 0.0007018037456432335 0.0006005948006725787 0.0005078784471340359 0.0004406418213749735 0.0003714124757418396 0.0003264057688551507 0.0003047322000871897 0.0002906911596356952 0.0002861738822980829 0.0002743316277907311 0.0002636405310368626 0.0002526665586205466 0.000226716276058895 0.0001875320324724495 0.0001636750239487869 0.0001372388125105794 0.0001222573562860838 0.0001269254248804013 0.0001139495649892648 0.0001106414381361491 0.0001196292706110303 0.0001189746129582091 0.0001112327955701268 0.0001296685790607515 0.0001643819608291414 0.0001878863422254414 0.0002344588678557017 0.0003203860709214426 0.0004064572654023451 0.0004940737860301076 0.000596014360602431 0.0006749848938371238 0.0007389096824511 0.0007792584703839535 0.0007965798134060041 0.0008056540909816971 0.0007849711152884964 0.000736062831596933 0.0006709962921592219 0.000603317655045247 0.0005383093164802969 0.0004900142937176716 0.0004551531490340078 0.000449534106387221 0.0004708650906301616 0.0004934795110626268 0.0005345974954791178 0.0006078660984573032 0.000702762715896311 0.0007775394691744047 0.0008818221619807017 0.001000679407785467 0.001092707546060794 0.001210613677564552 0.001312089019158298 0.001405270807169912 0.001490114661507714 0.001531394363834887 0.001555410794560944 0.001543704930528568 0.001470479573974331 0.001389350245308104 0.00128125012918577 0.001130987351527693 0.001007863016474597 0.0008564828937716682 0.0007074501288004987 0.0006213523466697201 0.0005713540643662852 0.0005583807693821872 0.0006196820100433617 0.0007361817663542591 0.0008554464717298295 0.0009758354537658702 0.001030406744297961 0.001035680319931336 0.001012121155964727 0.0009548313236383478 0.0008790095856581987 0.0008222269458738788 0.0007806888647762658 0.00070288663268893 0.0006454619358970265 0.0006054134743706602 0.000585536509799107 0.0005761856239940797 0.0005868746914752406 0.0006209602772533917 0.0006504193524742758 0.0006871616100603014 0.0007025287393104163 0.0007171995117985032 0.0007172937860848917 0.0006961439382863224 0.0006718024855287201 0.0006254758975065462 0.0005416137326971845 0.0004555547815183704 0.0003867779579310328 0.0003290178707355808 0.0002953031361596086 0.000255119760230264 0.0002488539502124204 0.0002677640501578254 0.0002652005008809371 0.0002906090024497945 0.0003166398471483209 0.0003089920401697223 0.0003114866028324989 0.0003133184110010134 0.0002840418132516435 0.0002386702280686958 0.0001892280180500572 0.0001115506030916713 1.983929251944428e-07 -0.0001282921891725907 -0.0002827272707108253 -0.0004368740799478221 -0.0005848618113809967 -0.0007194441046966012 -0.0008196524633909468 -0.0009048202133206396 -0.0009695404302984057 -0.00101680303425177 -0.00104331484535875 -0.001044670960196941 -0.001096308220052877 -0.001198928001528127 -0.001307622162194671 -0.001475710104821832 -0.001686034525628763 -0.001932017480014 -0.002224364121912755 -0.002509340882950428 -0.002795659773184621 -0.00309134921771924 -0.003339621248074061 -0.003527096039668997 -0.003713374797631359 -0.003865457109591278 -0.003984544055673171 -0.004135912398823348 -0.004315826807472575 -0.004496951439455919 -0.004688515788482826 -0.004898735004401356 -0.005110796780538788 -0.005339127441411028 -0.00553035210769078 -0.005705695069882694 -0.005870819376498557 -0.006004875496650137 -0.006110777442311181 -0.006164210930893388 -0.00618756905448379 -0.006096936349455484 -0.005939592020122337 -0.005693193052915928 -0.005354094079859258 -0.005027974684133828 -0.004692376525298168 -0.004435709093992968 -0.004274376328604747 -0.004208157300760209 -0.004223103959974963 -0.004327093561393999 -0.004497183833286816 -0.004613763151773346 -0.004718655531703288 -0.00476162995946275 -0.004733215833039038 -0.004717818557505252 -0.004685430263054647 -0.004754957289104762 -0.004977412914427918 -0.005345713165123564 -0.005956911360135976 -0.006682313612074329 -0.007507498625754038 -0.008207774465767762 -0.008390747806307929 -0.007774338884105022 -0.005813672313667084 -0.002020675163340676 0.003934699513891164 0.01218867179800687 0.0227189593674642 0.03474656425430857 0.0474559430852135 0.05969919927081149 0.06940044807963783 0.07690251527634465 0.08229165974889782 0.07690251527634465 0.06940044807963783 0.05969919927081149 0.04745594308521349 0.03474656425430857 0.02271895936746419 0.01218867179800687 0.003934699513891165 -0.002020675163340676 -0.005813672313667083 -0.007774338884105022 -0.008390747806307929 -0.008207774465767762 -0.007507498625754038 -0.006682313612074331 -0.005956911360135976 -0.005345713165123564 -0.004977412914427918 -0.004754957289104762 -0.004685430263054647 -0.004717818557505251 -0.004733215833039039 -0.00476162995946275 -0.004718655531703288 -0.004613763151773345 -0.004497183833286817 -0.004327093561393998 -0.004223103959974966 -0.004208157300760209 -0.004274376328604745 -0.00443570909399297 -0.004692376525298168 -0.005027974684133829 -0.00535409407985926 -0.005693193052915932 -0.005939592020122337 -0.006096936349455483 -0.006187569054483791 -0.006164210930893387 -0.006110777442311182 -0.006004875496650137 -0.005870819376498557 -0.005705695069882694 -0.00553035210769078 -0.005339127441411028 -0.005110796780538787 -0.004898735004401356 -0.004688515788482826 -0.004496951439455919 -0.004315826807472577 -0.004135912398823348 -0.003984544055673172 -0.003865457109591278 -0.00371337479763136 -0.003527096039668997 -0.003339621248074061 -0.00309134921771924 -0.00279565977318462 -0.002509340882950424 -0.002224364121912757 -0.001932017480014 -0.001686034525628767 -0.001475710104821831 -0.001307622162194672 -0.001198928001528127 -0.001096308220052873 -0.001044670960196942 -0.00104331484535875 -0.001016803034251772 -0.0009695404302984054 -0.00090482021332064 -0.0008196524633909476 -0.0007194441046966014 -0.0005848618113809963 -0.000436874079947822 -0.0002827272707108247 -0.0001282921891725906 1.9839292519393e-07 0.0001115506030916713 0.0001892280180500569 0.0002386702280686963 0.0002840418132516444 0.0003133184110010132 0.0003114866028324986 0.0003089920401697223 0.0003166398471483207 0.0002906090024497948 0.0002652005008809363 0.000267764050157825 0.0002488539502124205 0.0002551197602302644 0.000295303136159608 0.0003290178707355802 0.0003867779579310346 0.0004555547815183697 0.0005416137326971845 0.000625475897506546 0.0006718024855287199 0.0006961439382863226 0.0007172937860848909 0.0007171995117985036 0.0007025287393104162 0.0006871616100603015 0.000650419352474276 0.0006209602772533919 0.0005868746914752408 0.0005761856239940797 0.0005855365097991073 0.0006054134743706607 0.0006454619358970264 0.0007028866326889299 0.0007806888647762659 0.000822226945873879 0.0008790095856581987 0.0009548313236383482 0.001012121155964728 0.001035680319931336 0.001030406744297961 0.0009758354537658703 0.0008554464717298296 0.0007361817663542591 0.0006196820100433625 0.0005583807693821864 0.000571354064366285 0.0006213523466697193 0.0007074501288005006 0.0008564828937716656 0.001007863016474597 0.001130987351527696 0.001281250129185769 0.001389350245308103 0.001470479573974332 0.001543704930528569 0.001555410794560944 0.001531394363834887 0.001490114661507714 0.001405270807169912 0.001312089019158298 0.001210613677564552 0.001092707546060794 0.001000679407785467 0.0008818221619807018 0.0007775394691744045 0.0007027627158963109 0.0006078660984573031 0.000534597495479118 0.0004934795110626265 0.0004708650906301616 0.0004495341063872209 0.0004551531490340077 0.0004900142937176715 0.0005383093164802969 0.000603317655045247 0.0006709962921592213 0.0007360628315969327 0.0007849711152884964 0.0008056540909816968 0.0007965798134060039 0.000779258470383954 0.0007389096824511003 0.0006749848938371232 0.0005960143606024314 0.00049407378603011 0.0004064572654023451 0.000320386070921443 0.0002344588678557023 0.0001878863422254417 0.0001643819608291412 0.0001296685790607518 0.0001112327955701267 0.0001189746129582092 0.0001196292706110302 0.0001106414381361493 0.000113949564989265 0.0001269254248804011 0.000122257356286084 0.0001372388125105794 0.0001636750239487874 0.0001875320324724492 0.0002267162760588951 0.0002526665586205469 0.0002636405310368628 0.000274331627790731 0.0002861738822980834 0.0002906911596356955 0.0003047322000871902 0.0003264057688551512 0.0003714124757418389 0.0004406418213749726 0.0005078784471340354 0.000600594800672577 0.0007018037456432339 0.0007841007737282728 0.0008547924134192688 0.0009095689814133239 0.0009397344198433432 0.0009480193252716988 0.0009291650300966141 0.0008838704404048497 0.0008127224793692449 0.0007154307983328389 0.0006116301200159011 0.0005101781526047185 0.0004127132318400751 0.0003319979542964823 0.0002860821491880682 0.0002500563490184256 0.0002217724458095814 0.0002391073092942146 0.0002753731712067281 0.0002970986813756304 0.0003141523386480819 0.0003234961169067191 0.0003235563379009124 0.0003148588847685941 0.0002722059220446591 0.0002152831436118294 0.0001791264427517463 0.0001484107844269434 0.0001258477694353638 0.0001276356087955284 0.0001611481391302664 0.0002328092015126938 0.0003146483534597321 0.0004049366284488059 0.000515541119297898 0.0006278320617335942 0.0007308920792511362 0.0007963231371937704 0.0008367582390177949 0.000856084224659261 0.0008134213275591343 0.0007279470267812722 0.0006212972022330388 0.0004918678731923973 0.0003513571295408855 0.0002274300947131321 0.000126009265634115 4.952203678213429e-05 5.024961716295127e-06 3.538818509299518e-06 5.102070582076996e-05 0.0001215619643926977 0.0002314452035966865 0.0003647207877911044 0.0005014886072902379 0.0006440612144425216 0.0007712153995398789 0.0008785166874477058 0.0009565650123854939 0.0009991649016365179 0.001000137248327646 0.0009780458869915582 0.0009324133330452072 0.0008653125079193493 0.0008018307451331798 0.000746647621318018 0.0007004172297637317 0.0006853703272953076 0.0006924927262763831 0.0006931453291674521 0.0007101214361108594 0.0007121788880774427 0.000704014950979915 0.0006972231165628779 0.000661654211634347 0.0006231865628459464 0.0005890158012063288 0.0005446800229664074 0.0005043751480265481 0.0004885372867673724 0.0004761823498486873 0.0004638037190751788 0.0004587034220581102 0.0004684001252653539 0.0004858979454010797 0.0004993895797658201 0.0005175977188382752 0.0005269773541343305 0.0005292033760978143 0.0005322207128124349 0.0005317635475002128 0.000523675449512398 0.000513803426036456 0.0005013700140284385 0.0004800892227209432 0.0004485558688712597 0.0004148377173810098 0.0003761847353904126 0.0003355248671444355 0.0003084406671027549 0.0002862572870604022 0.000270260451250159 0.0002673721567286339 0.0002810007495450609 0.000304884990166382 0.0003328960310700833 0.00036358009885258 0.0003975165871572802 0.000425258701401412 0.0004465585690774157 0.0004622641649603433 0.0004668087315020517 0.0004716258022734012 0.0004674556394527063 0.0004612992161328973 0.0004555272479616062 0.000454810249957407 0.0004645950604270157 0.0004727638664349478 0.000486455053090451 0.0005041051313799224 0.0005240387399286972 0.0005366426137323651 0.0005326203262494113 0.0005159892932585147 0.0004796269124398812 0.0004241895316775592 0.0003565185832637645 0.0002783857988869775 0.0002020713027294521 0.0001300149878038001 7.787515033943645e-05 4.94970529116387e-05 3.567320177887565e-05 5.198681974621122e-05 8.439260614281255e-05 0.0001233684249499196 0.0001696438030090706 0.0002130987320308694 0.0002483396512410884 0.000269119930618064 0.0002800989971601369 0.0002797524248574673 0.0002574810386944714 0.0002293068332902553 0.0001981874567952985 0.0001541458098186209 0.0001131985409445768 7.436668389413293e-05 3.914274819454482e-05 2.02910310075359e-05 1.214910078795737e-05 1.138159045954142e-05 2.303431628630914e-05 4.191084339680085e-05 7.153140827090229e-05 0.0001012747411089004 0.0001198264422032543 0.0001341681487126989 0.0001353721674033473 0.0001310955784948295 0.0001304604164937583 0.0001293962113881343 0.0001167532618894624 0.0001090067194513343 0.0001142849115928038 0.0001194268883117053 0.0001316378760800131 0.0001450923278130218 0.0001644147124457893 0.0001836754090288532 0.0001945597427125205 0.000198420646821421 0.0001926926553236415 0.0001774138786603597 0.0001604569345523672 0.0001441800101209615 0.0001197702730190544 0.0001047837385634294 9.626406072217448e-05 8.422482885728265e-05 8.141197975598381e-05 8.7319365554045e-05 9.4075039188646e-05 9.894340278646086e-05 0.0001025516069602838 0.0001072060248103339 0.0001098275912662369 0.0001104312822684126 0.000106515806645663 9.505250826889828e-05 8.723833080263867e-05 7.420695272211161e-05 5.781977436721884e-05 4.770388200482809e-05 3.7389171346018e-05 3.012154305643296e-05 3.024288241806269e-05 3.30459157472601e-05 3.42023945226418e-05 4.194241650880383e-05 5.275816447026444e-05 6.085367944553199e-05 6.930393501059805e-05 8.209630166614197e-05 9.587608125739296e-05 0.0001056747747162764 0.0001170129681552745 0.0001293030383250257 0.0001415291461137908 0.0001570034141396783 0.0001704401347270341 0.0001743894362264053 0.0001736521706833197 0.0001685794934648447 0.0001573246833069577 0.00013725374822175 0.0001119097915371812 8.838658249276298e-05 6.180780852724682e-05 3.467670361993845e-05 1.383091083546099e-05 -1.363194189371131e-06 -1.629614429158657e-05 -2.559556014458577e-05 -2.719609151945906e-05 -2.432568386408844e-05 -1.745355895723421e-05 -6.139147851498011e-06 9.390120446350515e-06 2.609404178225465e-05 4.505424718694725e-05 6.397168698635465e-05 8.388579609557507e-05 0.0001017868860228118 0.0001164777490003533 0.0001282916901414377 0.0001371895926370353 0.0001444545226667357 0.0001505912434581831 0.0001573044997290388 0.0001611393815082451 0.0001617032194688336 0.0001578219868614449 0.000153801972814461 0.0001435371988464706 0.0001284375906324499 0.0001171407871684956 0.0001046664881685255 9.474089964535937e-05 8.653490906422336e-05 8.178131856442878e-05 8.236867203976517e-05 8.670764514445867e-05 9.282321723585233e-05 9.949247921871623e-05 0.0001078643622955885 0.0001145773485240607 0.0001199864757482572 0.0001237937282219144 0.0001274951568358209 0.0001276193402990368 0.0001211159686312227 0.0001157727179778552 0.0001100627878797457 9.993747197784049e-05 8.855370570715271e-05 8.085052742218457e-05 7.795597606716331e-05 7.47608349460523e-05 7.602754602719387e-05 8.465588232787989e-05 9.340521285980456e-05 0.0001021842594143064 0.0001128586850331394 0.0001211256678148737 0.0001268741533262061 0.0001314146843973354 0.0001298558067129026 0.0001232907290105823 0.0001114408649087925 9.490108720790098e-05 7.524935743597891e-05 5.098405707491407e-05 2.533067638032283e-05 2.754137907668097e-06 -1.425766597021359e-05 -2.680377848432208e-05 -3.140205230802843e-05 -2.889508983840295e-05 -2.038060526743043e-05 -5.011952829637147e-06 1.151301392357218e-05 2.733281113424968e-05 4.16286613756742e-05 5.28271491118769e-05 5.894036887788697e-05 5.965124091807064e-05 5.995001465356243e-05 6.151519405194192e-05 6.038440324493523e-05 5.810382275274606e-05 5.745029328387439e-05 5.754247327156637e-05 5.656595617636194e-05 5.387673663540817e-05 5.012294778018432e-05 4.529266752036116e-05 4.11193488460899e-05 - - diff --git a/octave/fmfsk.m b/octave/fmfsk.m deleted file mode 100644 index 4b3cc91..0000000 --- a/octave/fmfsk.m +++ /dev/null @@ -1,346 +0,0 @@ -% -% fmfsk.m -% Author: Brady O'Brien 3 Feb 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 veRbion 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 <http://www.gnu.org/licenses/>. - -% mancyfsk.m modem, extracted and made suitable for C implementation - -fm; -pkg load signal; -pkg load parallel; -1; - -% Init fmfsk modem -%Fs is sample frequency -%Rb is pre-manchester bit rate -function states = fmfsk_init(Fs,Rb) - assert(mod(Fs,Rb*2)==0); - - %Current fixed processing buffer size, in non-ME bits - nbit = 96; - - states.Rb = Rb; - states.Rs = Rb*2; % Manchester-encoded bitrate - states.Fs = Fs; - states.Ts = Fs/states.Rs; - states.N = nbit*2*states.Ts; - states.nin = states.N; % Samples in the next demod cycle - states.nstash = states.Ts*2; % How many samples to stash away between proc cycles for timing adjust - states.nmem = states.N+(4*states.Ts); - states.nsym = nbit*2; - states.nbit = nbit; - - %tates.nsym = floor(states.Rs*.080); - %states.nbit = floor(states.Rb*.080) - %Old sample memory - - states.oldsamps = zeros(1,states.nmem); - - %Last sampled-stream output, for odd bitstream generation - states.lastint = 0; - - %Some stats - states.norm_rx_timing = 0; - -endfunction - -%Generate a stream of manchester-coded bits to be sent -% to any ordinary FM modulator or VCO or something -function tx = fmfsk_mod(states,inbits) - Ts = states.Ts; - tx = zeros(1,length(inbits)*2); - for ii = 1:length(inbits) - st = 1 + (ii-1)*Ts*2; - md = st+Ts-1; - en = md+Ts; - if inbits(ii)==0 - tx(st:md) = -ones(1,Ts); - tx(md+1:en) = ones(1,Ts); - else - tx(st:md) = ones(1,Ts); - tx(md+1:en) = -ones(1,Ts); - end - end -endfunction - -%Demodulate a bag of bits from the output of an FM demodulator -% This function produces nbits output bits and takes states.nin samples -function [rx_bits states] = fmfsk_demod(states,rx) - Ts = states.Ts; - Fs = states.Fs; - Rs = states.Rs; - nin = states.nin; - N = states.N; - nsym = states.nsym; - nbits = states.nsym/2; - nmem = states.nmem; - nstash = states.nstash; - - nold = nmem-nin; - ssamps = states.oldsamps; - - - %Shift in nin samples - ssamps(1:nold) = ssamps(nmem-nold+1:nmem); - ssamps(nold+1:nmem) = rx; - states.oldsamps = ssamps; - - rx_filt = zeros(1,(nsym+1)*Ts); - %Integrate Ts input samples at every offset - %This is the same thing as filtering with a filter of all ones - % out to Ts. - % It's implemented like this for ease of C-porting - for ii=(1:(nsym+1)*Ts) - st = ii; - en = st+Ts-1; - rx_filt(ii) = sum(ssamps(st:en)); - end - states.rx_filt = rx_filt; - % Fine timing estimation ------------------------------------------------------ - - % Estimate fine timing using line at Rs/2 that Manchester encoding provides - % We need this to sync up to Manchester codewords. - Np = length(rx_filt); - w = 2*pi*(Rs)/Fs; - x = (rx_filt .^ 2) * exp(-j*w*(0:Np-1))'; - norm_rx_timing = angle(x)/(2*pi)-.42; - - rx_timing = round(norm_rx_timing*Ts); - - %If rx timing is too far out, ask for more or less sample the next time - % around to even it all out - next_nin = N; - if norm_rx_timing > -.2; - next_nin += Ts/2; - end - if norm_rx_timing < -.65; - next_nin -= Ts/2; - end - - states.nin = next_nin; - states.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 - lastv = states.lastint; - for ii = (0:nsym-1) - currv = rx_filt(sample_offset+(ii*Ts)+1); - 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; - else - apodd += abs(mdiff); - rx_odd( floor(ii/2)+1 ) = mbit; - end - end - states.symsamp = symsamp; - % Decide on the correct ME alignment - if(apeven>apodd) - rx_bits = rx_even; - else - rx_bits = rx_odd; - end - - states.lastint = lastv; -endfunction - -% run_sim copypasted from fsk_horus.m -% simulation of tx and rx side, add noise, channel impairments ---------------------- - -function fmfsk_run_sim(EbNodB,timing_offset=0,de=0,of=0,hpf=0) - 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 - - load fm_radio_filt_model.txt - - [b, a] = cheby1(4, 1, 300/Fs, 'high'); % 300Hz HPF to simulate FM radios - - tx_pmod = fmfsk_mod(states, tx_bits); - 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 = filter(filt,1,rx); - - figure(4) - title("Spectrum of rx-ed signal after FM demod and FM radio channel"); - plot(20*log10(abs(fft(rx)))) - figure(5) - title("Time domain of rx-ed signal after FM demod and FM radio channel"); - plot(rx) - %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 = []; - - for f=1:frames - - % extract nin samples from input stream - - nin = states.nin; - en = st + states.nin - 1; - sf = rx(st:en); - st += nin; - - % demodulate to stream of bits - - [rx_bits states] = fmfsk_demod(states, sf); - - 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]; - - end - - ber = 1; - ox = 1; - rx_bits = rx_bits_log; - bitcnt = length(tx_bits); - 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; - figure(3); - plot(xor(rx_bits(ox:length(rx_bits)),tx_bits(1:length(rx_bits)+1-ox))) - - printf("BER: %f Errors: %d Bits:%d\n",ber,best_nerr,bitcnt-offset); - - endfunction - - -% demodulate a file of 8kHz 16bit short samples -------------------------------- - - - - diff --git a/octave/fsk4_dmr.m b/octave/fsk4_dmr.m deleted file mode 100644 index e6ccb7a..0000000 --- a/octave/fsk4_dmr.m +++ /dev/null @@ -1,540 +0,0 @@ -% fsk4.m -% -% Brady O'Brien October 2015 -% -% 4FSK modem attempt from the DMR spec - -graphics_toolkit("gnuplot"); - -fm; % analog FM modulator functions - -pkg load signal; - -% Init function for modem ------------------------------------------------------------ - -function fsk4_states = fsk4_init(fsk4_states,fsk4_info) - Fs = fsk4_states.Fs = 48000; %Sample rate - Rs = fsk4_states.Rs = fsk4_info.rs; %Symbol rate - M = fsk4_states.M = fsk4_states.Fs/fsk4_states.Rs; %Samples per symbol - - % Set up 4FSK raised cosine filter. This probably screws up perf if we were using - % optimal mod and dmeods but helps performance when using nasty old analog FM mods - % and demods - - empty_filter = [zeros(1,99) 1]; - - rf = (0:(Fs/2)); - %If there's no filter with this modem configuration, don't bother generating one - if fsk4_info.no_filter - fsk4_states.tx_filter = empty_filter; - fsk4_states.rx_filter = empty_filter; - else - fsk4_states.tx_filter = fir2(400 ,rf/(Fs/2),fsk4_info.tx_filt_resp(rf)); - fsk4_states.rx_filter = fir2(400 ,rf/(Fs/2),fsk4_info.rx_filt_resp(rf)); - endif - - %fsk4_states.tx_filter = fsk4_states.rx_filter = [zeros(1,99) 1]; - %Set up the 4FSK symbols - fsk4_states.symmap = fsk4_info.syms / fsk4_info.max_dev; - - fm_states.Ts = M; - fm_states.Fs = Fs; - fm_states.fc = 0; - fm_states.fm_max = fsk4_info.max_dev*2; - fm_states.fd = fsk4_info.max_dev; - fm_states.pre_emp = fm_states.de_emp = 0; - fm_states.output_filter = 0; - fm_states.ph_dont_limit = 1; - fsk4_states.fm_states = analog_fm_init(fm_states); - fsk4_states.modinfo = fsk4_info; - fsk4_states.verbose = 0; -endfunction - -%Integrate over data and dump every M samples -function d = idmp(data, M) - d = zeros(1,length(data)/M); - for i = 1:length(d) - d(i) = sum(data(1+(i-1)*M:i*M)); - end -endfunction - - -% DMR modulator ---------------------------------------------------------- - -function [tx, tx_filt, tx_stream] = fsk4_mod(fsk4_states, tx_bits) - verbose = fsk4_states.verbose - - hbits = tx_bits(1:2:length(tx_bits)); - lbits = tx_bits(2:2:length(tx_bits)); - %Pad odd bit lengths - if(length(hbits)!=length(lbits)) - lbits = [lbits 0]; - end - tx_symbols = lbits + hbits*2 + 1; - M = fsk4_states.M; - nsym = length(tx_symbols); - nsam = nsym*M; - - tx_stream = zeros(1,nsam); - for i=1:nsym - tx_stream(1+(i-1)*M:i*M) = fsk4_states.symmap(tx_symbols(i)); - end - tx_filt = filter(fsk4_states.tx_filter, 1, tx_stream); - tx = analog_fm_mod(fsk4_states.fm_states, tx_filt); - - if verbose - figure(10); - plot(20*log10(abs(fft(tx)))) - title("Spectrum of modulated 4FSK") - endif - -endfunction - - -% Integrate and Dump 4FSK demod ---------------------------------------------------- - -function bits = fsk4_demod_thing(fsk4_states, rx) - - M = fsk4_states.M; - Fs = fsk4_states.Fs; - verbose = fsk4_states.verbose; - t = (0:length(rx)-1); - symup = fsk4_states.modinfo.syms; - - % Integrator is like an FIR filter with rectangular window coeffs. - % This has some nasty side lobes so lets limit the overall amount - % of noise getting in. tx filter just happens to work, but I imagine - % other LPF would as well. - - Fs = fsk4_states.Fs; - rf = (0:(Fs/2)); - rx_filter_a = fir1(100 ,.2); - rx_filter_b = fsk4_states.rx_filter; - rx_filter_n = [zeros(1,99) 1]; - - rx = filter(rx_filter_b, 1, rx); - - sym1m = exp(-j*2*pi*(symup(1)/Fs)*t).*rx; - sym2m = exp(-j*2*pi*(symup(2)/Fs)*t).*rx; - sym3m = exp(-j*2*pi*(symup(3)/Fs)*t).*rx; - sym4m = exp(-j*2*pi*(symup(4)/Fs)*t).*rx; - - % this puppy found by experiment between 1 and M. Will vary with different - % filter impulse responses, as delay will vary. f you add M to it coarse - % timing will adjust by 1. - - fine_timing = 54; - - sym1m = idmp(sym1m(fine_timing:length(sym1m)),M); sym1m = (real(sym1m).^2+imag(sym1m).^2); - sym2m = idmp(sym2m(fine_timing:length(sym2m)),M); sym2m = (real(sym2m).^2+imag(sym2m).^2); - sym3m = idmp(sym3m(fine_timing:length(sym3m)),M); sym3m = (real(sym3m).^2+imag(sym3m).^2); - sym4m = idmp(sym4m(fine_timing:length(sym4m)),M); sym4m = (real(sym4m).^2+imag(sym4m).^2); - - - figure(2); - nsym = 500; - %subplot(411); plot(sym1m(1:nsym)) - %subplot(412); plot(sym2m(1:nsym)) - %subplot(413); plot(sym3m(1:nsym)) - %subplot(414); plot(sym4m(1:nsym)) - plot((1:nsym),sym1m(1:nsym),(1:nsym),sym2m(1:nsym),(1:nsym),sym3m(1:nsym),(1:nsym),sym4m(1:nsym)) - - [x iv] = max([sym1m; sym2m; sym3m; sym4m;]); - bits = zeros(1,length(iv*2)); - figure(3); - hist(iv); - for i=1:length(iv) - bits(1+(i-1)*2:i*2) = [[0 0];[0 1];[1 0];[1 1]](iv(i),(1:2)); - end -endfunction - -function dat = bitreps(in,M) - dat = zeros(1,length(in)*M); - for i=1:length(in) - dat(1+(i-1)*M:i*M) = in(i); - end -endfunction - -% Minimal Running Disparity, 4 symbol encoder -% This is a simple 1 bit to 1 symbol encoding for 4fsk modems built -% on old fashoned FM radios. -function syms = mrd4(bits) - syms = zeros(1,length(bits)); - rd=0; - lastsym=0; - for n = (1:length(bits)) - bit = bits(n); - sp = [1 3](bit+1); %Map a bit to a +1 or +3 - [x,v] = min(abs([rd+sp rd-sp])); %Select +n or -n, whichever minimizes disparity - ssel = [sp -sp](v); - if(ssel == lastsym)ssel = -ssel;endif %never run 2 of the same syms in a row - syms(n) = ssel; %emit the symbol - rd = rd + ssel; %update running disparity - lastsym = ssel; %remember this symbol for next time - end -endfunction - -% Minimal Running Disparity, 8 symbol encoder -% This is a simple 2 bit to 1 symbol encoding for 8fsk modems built -% on old fashoned FM radios. -function syms = mrd8(bits) - bitlen = length(bits); - if mod(bitlen,2) == 1 - bits = [bits 0] - endif - - syms = zeros(1,length(bits)*.5); - rd=0; - lastsym=0; - for n = (1:2:length(bits)) - bit = (bits(n)*2)+bits(n+1); - sp = [1 3 7 5](bit+1); %Map a bit to a +1 or +3 - [x,v] = min(abs([rd+sp rd-sp])); %Select +n or -n, whichever minimizes disparity - ssel = [sp -sp](v); - if(ssel == lastsym)ssel = -ssel;endif %never run 2 of the same syms in a row - syms((n+1)/2) = ssel; %emit the symbol - rd = rd + ssel; %update running disparity - lastsym = ssel; %remember this symbol for next time - end -endfunction - -% "Manchester 4" encoding -function syms = mane4(bits) - syms = zeros(1,floor(bits/2)*2); - for n = (1:2:length(bits)) - bit0 = bits(n); - bit1 = bits(n+1); - sel = 2*bit0+bit1+1; - syms(n:n+1) = [[3 -3];[-3 3];[1 -1];[-1 1]]( sel,(1:2) ); - end -endfunction - -function out = fold_sum(in,l) - sublen = floor(length(in)/l); - out = zeros(1,l); - for i=(1:sublen) - v = in(1+(i-1)*l:i*l); - out = out + v; - end -endfunction - -function [bits err rxphi] = fsk4_demod_fmrid(fsk4_states, rx, enable_fine_timing = 0) - %Demodulate fsk signal with an analog fm demod - rxd = analog_fm_demod(fsk4_states.fm_states,rx); - - M = fsk4_states.M; - verbose = fsk4_states.verbose; - %This is the ideal fine timing, assuming the same offset in nfbert - fine_timing = 61; - - %This is meant to be adjusted by the fine timing estimator. comment out for - %ideal timing - %fine_timing = 59; - - %RRC filter to get rid of some of the noise - rxd = filter(fsk4_states.rx_filter, 1, rxd); - - %Try and figure out where sampling should happen over 30 symbol periods - diffsel = fold_sum(abs(diff( rxd(3001:3001+(M*30)) )),10); - - if verbose - figure(11); - plot(diffsel); - title("Fine timing estimation"); - endif - - %adjust fine timing - [v iv] = min(diffsel); - if enable_fine_timing - fine_timing = 59 + iv; - endif - rxphi = iv; - - %sample symbols - sym = rxd(fine_timing:M:length(rxd)); - - if verbose - figure(4) - plot(sym(1:1000)); - title("Sampled symbols") - endif - %eyediagram(afsym,2); - % Demod symbol map. I should probably figure a better way to do this. - % After sampling, the furthest symbols tend to be distributed about .80 - - % A little cheating to demap the symbols - % Take a histogram of the sampled symbols, find the center of the largest distribution, - % and correct the symbol map to match it - [a b] = hist(abs(sym),50); - [a ii] = max(a); - %grmax = abs(b(ii)); - %grmax = (grmax<.65)*.65 + (grmax>=.65)*grmax; - grmax = .84; - dmsyms = rot90(fsk4_states.symmap*grmax) - (dmsyms(2)+dmsyms(1))/2 - - if verbose - figure(2) - hist(abs(sym),200); - title("Sampled symbol histogram") - endif - - %demap the symbols - [err, symout] = min(abs(sym-dmsyms)); - - if verbose - figure(3) - hist(symout); - title("De-mapped symbols") - endif - - bits = zeros(1,length(symout)*2); - %Translate symbols back into bits - - for i=1:length(symout) - bits(1+(i-1)*2:i*2) = [[1 1];[1 0];[0 1];[0 0]](symout(i),(1:2)); - end -endfunction - -% Frequency response of the DMR raised cosine filter -% from ETSI TS 102 361-1 V2.2.1 page 111 -dmr.tx_filt_resp = @(f) sqrt(1.0*(f<=1920) - cos((pi*f)/1920).*1.0.*(f>1920 & f<=2880)); -dmr.rx_filt_resp = dmr.tx_filt_resp; -dmr.max_dev = 1944; -dmr.syms = [-1944 -648 1944 648]; -dmr.rs = 4800; -dmr.no_filter = 0; -dmr.demod_fx = @fsk4_demod_fmrid; -global dmr_info = dmr; - - -% No-filter 4FSK 'ideal' parameters -nfl.tx_filt_resp = @(f) 1; -nfl.rx_filt_resp = nfl.tx_filt_resp; -nfl.max_dev = 7200; -%nfl.syms = [-3600 -1200 1200 3600]; -nfl.syms = [-7200,-2400,2400,7200]; -nfl.rs = 4800; -nfl.no_filter = 1; -nfl.demod_fx = @fsk4_demod_thing; -global nflt_info = nfl; - -%Some parameters for the NXDN filters -nxdn_al = .2; -nxdn_T = 416.7e-6; -nxdn_fl = ((1-nxdn_al)/(2*nxdn_T)); -nxdn_fh = ((1+nxdn_al)/(2*nxdn_T)); - -%Frequency response of the NXDN filters -% from NXDN TS 1-A V1.3 page 13 -% Please note : NXDN not fully implemented or tested -nxdn_H = @(f) 1.0*(f<nxdn_fl) + cos( (nxdn_T/(4*nxdn_al))*(2*pi*f-(pi*(1-nxdn_al)/nxdn_T)) ).*(f<=nxdn_fh & f>nxdn_fl); -nxdn_P = @(f) (f<=nxdn_fh & f>0).*((sin(pi*f*nxdn_T))./(.00001+(pi*f*nxdn_T))) + 1.0*(f==0); -nxdn_D = @(f) (f<=nxdn_fh & f>0).*((pi*f*nxdn_T)./(.00001+sin(pi*f*nxdn_T))) + 1.0*(f==0); - -nxdn.tx_filt_resp = @(f) nxdn_H(f).*nxdn_P(f); -nxdn.rx_filt_resp = @(f) nxdn_H(f).*nxdn_D(f); -nxdn.rs = 4800; -nxdn.max_dev = 1050; -nxdn.no_filter = 0; -nxdn.syms = [-1050,-350,350,1050]; -nxdn.demod_fx = @fsk4_demod_fmrid; -global nxdn_info = nxdn; - -% Bit error rate test ---------------------------------------------------------- -% Params - aEsNodB - EbNo in decibels -% - timing_offset - how far the fine timing is offset -% - bitcnt - how many bits to check -% - demod_fx - demodulator function -% Returns - ber - the measured BER -% - thrcoh - theory BER of a coherent demod -% - thrncoh - theory BER of non-coherent demod -function [ber thrcoh thrncoh] = nfbert(aEsNodB,modem_config, bitcnt=100000, timing_offset = 10) - - rand('state',1); - randn('state',1); - - %How many bits should this test run? - bitcnt = 120000; - - test_bits = [zeros(1,100) rand(1,bitcnt)>.5]; %Random bits. Pad with zeros to prime the filters - fsk4_states.M = 1; - fsk4_states = fsk4_init(fsk4_states,modem_config); - - %Set this to 0 to cut down on the plotting - fsk4_states.verbose = 1; - Fs = fsk4_states.Fs; - Rb = fsk4_states.Rs * 2; % Multiply symbol rate by 2, since we have 2 bits per symbol - - tx = fsk4_mod(fsk4_states,test_bits); - - %add noise here - %shamelessly copied from gmsk.m - EsNo = 10^(aEsNodB/10); - EbNo = EsNo - variance = Fs/(Rb*EbNo); - nsam = length(tx); - noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); - rx = tx*exp(j*pi/2) + noise; - - rx = rx(timing_offset:length(rx)); - - rx_bits = modem_config.demod_fx(fsk4_states,rx); - ber = 1; - - %thing to account for offset from input data to output data - %No preamble detection yet - ox = 1; - for offset = (1:100) - nerr = sum(xor(rx_bits(offset:length(rx_bits)),test_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; - printf("\ncoarse timing: %d nerr: %d\n", offset, best_nerr); - - % Coherent BER theory - thrcoh = erfc(sqrt(EbNo)); - - % non-coherent BER theory calculation - % It was complicated, so I broke it up - - ms = 4; - 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 )); - -endfunction - -% RX fine timing estimation playground -function rxphi = fine_ex(timing_offset = 1) - global dmr_info; - global nxdn_info; - global nflt_info; - - rand('state',1); - randn('state',1); - - bitcnt = 12051; - test_bits = [zeros(1,100) rand(1,bitcnt)>.5]; %Random bits. Pad with zeros to prime the filters - t_vec = [0 0 1 1]; - %test_bits = repmat(t_vec,1,ceil(24000/length(t_vec))); - - - fsk4_states.M = 1; - fsk4_states = fsk4_init(fsk4_states,dmr_info); - Fs = fsk4_states.Fs; - Rb = fsk4_states.Rs * 2; %Multiply symbol rate by 2, since we have 2 bits per symbol - - tx = fsk4_mod(fsk4_states,test_bits); - - %add noise here - %shamelessly copied from gmsk.m - %EsNo = 10^(aEsNodB/10); - %EbNo = EsNo - %variance = Fs/(Rb*EbNo); - %nsam = length(tx); - %noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); - %rx = tx*exp(j*pi/2) + noise; - rx = tx; - rx = rx(timing_offset:length(rx)); - - [rx_bits biterr rxphi] = fsk4_demod_fmrid(fsk4_states,rx); - ber = 1; - - %thing to account for offset from input data to output data - %No preamble detection yet - ox = 1; - for offset = (1:100) - nerr = sum(xor(rx_bits(offset:length(rx_bits)),test_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; - printf("\ncoarse timing: %d nerr: %d\n", offset, best_nerr); - -endfunction - -%Run over a wide range of offsets and make sure fine timing makes sense -function fsk4_rx_phi(socket) - %pkg load parallel - offrange = [100:200]; - [a b c phi] = pararrayfun(1.25*nproc(),@nfbert,10*length(offrange),offrange); - - close all; - figure(1); - clf; - plot(offrange,phi); -endfunction - - -% Run this function to compare the theoretical 4FSK modem performance -% with our DMR modem simulation - -function fsk4_ber_curves - global dmr_info; - global nxdn_info; - global nflt_info; - - EbNodB = 1:20; - bers_tco = bers_real = bers_tnco = bers_idealsim = ones(1,length(EbNodB)); - - %vectors of the same param to pass into pararrayfun - dmr_infos = repmat(dmr_info,1,length(EbNodB)); - nflt_infos = repmat(nflt_info,1,length(EbNodB)); - thing = @fsk4_demod_thing; - - % Lovely innovation by Brady to use all cores and really speed up the simulation - - %try - pkg load parallel - bers_idealsim = pararrayfun(floor(1.25*nproc()),@nfbert,EbNodB,nflt_infos); - [bers_real,bers_tco,bers_tnco] = pararrayfun(floor(1.25*nproc()),@nfbert,EbNodB,dmr_infos); - %catch - % printf("You should install package parallel. It'll make this run way faster\n"); - % for ii=(1:length(EbNodB)); - %[bers_real(ii),bers,tco(ii),bers_tnco(ii)] = nfbert(EbNodB(ii)); - % end - %end_try_catch - - close all - figure(1); - clf; - semilogy(EbNodB, bers_tnco,'r;4FSK non-coherent theory;') - hold on; - - semilogy(EbNodB, bers_tco,'b;4FSK coherent theory;') - semilogy(EbNodB, bers_real ,'g;4FSK DMR simulation;') - semilogy(EbNodB, bers_idealsim, 'v;FSK4 Ideal Non-coherent simulation;') - hold off; - grid("minor"); - axis([min(EbNodB) max(EbNodB) 1E-5 1]) - legend("boxoff"); - xlabel("Eb/No (dB)"); - ylabel("Bit Error Rate (BER)") - -endfunction - - - - - - - - diff --git a/octave/fsk_basic.m b/octave/fsk_basic.m deleted file mode 100644 index 753f4ae..0000000 --- a/octave/fsk_basic.m +++ /dev/null @@ -1,60 +0,0 @@ -% fsk_basic.m -% David Rowe 30 sep 2016 -% -% Basic non-coherent FSK modem simulation to illustrate principles -% and compare to ideal - -rand('seed',1); -randn('seed',1); - -Fs = 9600; % sample rate -f1 = 1200; -f2 = 2400; -Rs = 1200; % symbol rate -Ts = Fs/Rs; % length of each symbol in samples -Nbits = 10000; -EbNodB = 9; - -tx_bits = round(rand(1,Nbits)); - -% continuous phase FSK modulator - -w1 = 2*pi*f1/Fs; -w2 = 2*pi*f2/Fs; -tx_phase = 0; -tx = zeros(1,Ts*Nbits); - -for i=1:Nbits - for k=1:Ts - if tx_bits(i) - tx_phase += w2; - else - tx_phase += w1; - end - tx((i-1)*Ts+k) = exp(j*tx_phase); - end -end - -% AWGN channel noise - -EbNo = 10^(EbNodB/10); -variance = Fs/(Rs*EbNo); -noise = sqrt(variance/2)*(randn(1,Nbits*Ts) + j*randn(1,Nbits*Ts)); -rx = tx + noise; - -% integrate and dump demodulator - -rx_bits = zeros(1,Nbits); -for i=1:Nbits - arx_symb = rx((i-1)*Ts + (1:Ts)); - filt1 = sum(exp(-j*w1*(1:Ts)) .* arx_symb); - filt2 = sum(exp(-j*w2*(1:Ts)) .* arx_symb); - rx_bits(i) = filt2 > filt1; -end - -Nerrors = sum(xor(tx_bits, rx_bits)); -ber = Nerrors/Nbits; -printf("EbNodB: %4.1f Nerrors: %d BER: %1.3f\n", EbNodB, Nerrors, ber); - - - diff --git a/octave/fsk_cml.m b/octave/fsk_cml.m deleted file mode 100644 index 25ca790..0000000 --- a/octave/fsk_cml.m +++ /dev/null @@ -1,132 +0,0 @@ -% Test MFSK at symbol level in CML using RA LDPC codes, Bill, June 2020
-% Simulate in AWGM and plot BERs. Assumes that CML is installed!
-%
-% If required setup numb of codewords (Ncw) and channel bits/sym (bps, 1 to 4)
-% may also select FEC code with Ctype (1 to 3); use plt to for debug plots of LLR values
-
-% July 1 version allows M=8 and pads required vectors if required to makeup 3rd bit
-
-ldpc;
-
-%rand('seed',1);
-%randn('seed',1);
-format short g
-more off
-init_cml();
-
-if exist('Ncw')==0, Ncw=100, end %setup defaults
-if exist('plt')==0, plt=0; end
-if exist('bps')==0, bps=4; end
-if exist('Ctype')==0, Ctype=1, end
-
-if Ctype==1
- load H_256_768_22.txt; HRA = H_256_768_22; % rate 1/3
-elseif Ctype==2
- load H_256_512_4.mat; HRA=H; % K=256, rate 1/2 code
- % above code might be improved -- but still works better than rate 1/3
-elseif Ctype==3
- load HRAa_1536_512.mat; % rate 3/4, N=2k code
-else
- error('bad Ctype');
-end
-
-M=2^bps; nos =0; clear res
-modulation = 'FSK'; mod_order=M; mapping = 'gray';
-code_param = ldpc_init_user(HRA, modulation, mod_order, mapping);
-
-[Nr Nc] = size(HRA);
-Nbits = Nc - Nr;
-Krate = (Nc-Nr)/Nc
-framesize = Nc;
-%{
-[H_rows, H_cols] = Mat2Hrows(HRA);
-code_param.H_rows = H_rows;
-code_param.H_cols = H_cols;
-code_param.P_matrix = [];
-%}
-
-
-S =CreateConstellation('FSK', M);
-if M==2, Ebvec=[7:0.2: 8.7], end
-if M==4, Ebvec=[4:0.25: 7], end
-if M==8, Ebvec =[3.5: 0.25: 5.5]; end
-if M==16, Ebvec=[2.5:0.25:4.8]; end
-
-disp(['Symbol-based ' num2str(M) 'FSK sim with K=' ...
- num2str(Nbits) ', code rate=' num2str(Krate) ', #CWs=' num2str(Ncw)])
-
-
-
-% if M=8, for 3 bits/symbol, may need to pad codeword with bits ...
-if floor(Nc/bps)*bps ~= Nc
- Npad = ceil(Nc/bps) *bps-Nc
- disp('padding codeword')
-else
- Npad=0;
-end
-
-for Eb = Ebvec
-
- Ec = Eb + 10*log10(Krate);
- Es = Ec + 10*log10(bps);
- Eslin = 10^(Es/10); %Es/N0 = 1/2k_n^2
-
- Terrs =0;
- for nn = 1:Ncw
-
- txbits = randi(2,1,Nbits) -1;
-
- codeword = LdpcEncode( txbits, code_param.H_rows, code_param.P_matrix );
- code_param.code_bits_per_frame = length( codeword );
- code_param.data_bits_per_frame = length(txbits);
- Nsymb = (code_param.code_bits_per_frame+Npad)/bps;
-
- if Npad; codeword = [codeword zeros(1,Npad)]; end
- Tx = Modulate(codeword, S);
-
- kn = sqrt(1/(2*Eslin));
- Rx = Tx + kn * (randn(size(Tx)) + j*randn(size(Tx)));
-
- SNRlin = Eslin; % Valenti calls this SNR, but seems to be Es/N0
- symL = DemodFSK(Rx, SNRlin, 2); %demod type is nonCOH, without estimate amplitudes
- bitL = Somap(symL);
-
- if Npad, bitL(end-Npad+1:end)=[]; end
- if plt>0, figure(110); hist(bitL); title('bit LLRs')
- figure(111); hist(bitL); title('Sym Ls'), pause,
- end
- max_it =100; decoder_type =0;
-
- [x_hat, PCcnt] = MpDecode( -bitL, code_param.H_rows, code_param.H_cols, ...
- max_it, decoder_type, 1, 1);
- Niters = sum(PCcnt~=0);
- detected_data = x_hat(Niters,:);
- error_positions = xor( detected_data(1:code_param.data_bits_per_frame), txbits );
-
- Nerrs = sum( error_positions);
- if plt>1, figure(121); plot(error_positions); Nerrs, end
- Terrs = Terrs + Nerrs;
- end
-
- BER = Terrs/ (Ncw*Nbits);
-
- %HDs = (sign(bitL)+1)/2;
- %NerrsHD = sum(codeword~=HDs);
- %BER_HD = Nerrs/Nbits;
-
- nos = nos+1;
- disp('Eb Nerrs BER')
- res(nos, :) = [Eb, Terrs, BER]
-
-end
-figure(91)
-semilogy(res(:,1), res(:,3), '-x'); grid on; hold on;
-%semilogy(res(:,1), berawgn(res(:,1), 'fsk', M, 'noncoherent'), 'g');
-title([num2str(M) 'FSK BER with LDPC FEC'])
-xlabel('Eb/N0'); ylabel('BER')
-
-
-
-
-
-
diff --git a/octave/fsk_cml_sam.m b/octave/fsk_cml_sam.m deleted file mode 100644 index f7d02f9..0000000 --- a/octave/fsk_cml_sam.m +++ /dev/null @@ -1,376 +0,0 @@ -%fsk_llr_sam Test MFSK using David's sample-based mod/demod with CML LLR routines -% using 2k rate 3/4 RA LDPC code. Bill, July 2020 -% -%LLR conversion options: (Ltype) -% 1 David's original M=2 SD to LLR routine -% 2 Bill's pdf-based M=2 or 4 SD to LLR -% 3 Some simple HD conversions -% 4 CML approach using symbol likelihoods, then converted to bit LLRs -% -% Results from this sim are stored in "res" -- use fsk_llr_plot to see BER figs -% Use "plt" to see some useful plots (for selected Ltypes) : -% eg 1 is SD pdf histograms; 2 is Rx PSD; 3 is bit LLR histograms -% -% Adjust Evec and Nbits as required before running. - -#{ - Example 1 2FSK: - octave:4> fsk_cml_sam - octave:5> fsk_llr_plot - - Example 2 4FSK: - octave:4> M=4; fsk_cml_sam - octave:4> fsk_llr_plot -#} - -ldpc; - -% define Rician pdf -% note that Valenti uses an approximation that avoids Bessel evaluation -function y = rice(x,v,s) - s2 = s*s; - y = (x / s2) .* exp(-0.5 * (x.^2 + v.^2)/s2) .* besseli(0, x*v/s2); -endfunction - -function plot_pdf(v,s) - x=(0:0.1:2*v); - y= rice(x, v, s); - figure(201); hold on - plot(x,y,'g'); - %title('Rician pdf: signal carrier') - y= rice(x, 0, s); - plot(x,y,'b'); - title('Rician pdf: signal and noise-only carriers') - pause(0.01); -endfunction - -% single Eb/No point simulation --------------- -function [raw_ber rx_filt rx_bits tx_symbols demapper sig_est SNRest v_est] = ... - run_single(tx_bits, M, EcNodB, plt) - % Ec/N0 is per channel bit - bps = log2(M); % bits per symbol - Ts = 16; % length of each symbol in samples - - if length(tx_bits)==1 - Nbits = tx_bits; - tx_bits = randi(2,1,Nbits)-1; % make random bits - end - - Nbits = length(tx_bits); - Nsymbols = Nbits/log2(M); - tx_symbols = zeros(1,Nsymbols); - - mapper = bps:-1:1; - % look up table demapper from symbols to bits (hard decision) - demapper=zeros(M,bps); - for m=1:M - for b=1:bps - if bitand(m-1,b) demapper(m,bps-b+1) = 1; end - end - end - - % continuous phase mFSK modulator - - w(1:M) = 2*pi*(1:M)/Ts; - tx_phase = 0; - tx = zeros(1,Ts*Nsymbols); - - for s=1:Nsymbols - bits_for_this_symbol = tx_bits(bps*(s-1)+1:bps*s); - symbol_index = bits_for_this_symbol * mapper' + 1; - tx_symbols(s) = symbol_index; - assert(demapper(symbol_index,:) == bits_for_this_symbol); - for k=1:Ts - tx_phase = tx_phase + w(symbol_index); - tx((s-1)*Ts+k) = exp(j*tx_phase); - end - end - - % AWGN channel noise - - EsNodB = EcNodB + 10*log10(bps); - EsNo = 10^(EsNodB/10); - variance = Ts/EsNo; - noise = sqrt(variance/2)*(randn(1,Nsymbols*Ts) + j*randn(1,Nsymbols*Ts)); - rx = tx + noise; - - if plt==2, % check the Spectrum - [psd,Fpsd] =pwelch(rx,128,0.5,128,Ts); - figure(110); plot(Fpsd,10*log10(psd)); - title('Rx Signal: PSD '); - xlabel('Freq/Rs'); - %figure(111);plot(unwrap(arg(tx))); - pause(0.01); - end - - - % integrate and dump demodulator - - rx_bits = zeros(1,Nbits); - rx_filt = zeros(Nsymbols,M); - rx_pows = zeros(1,M); - rx_nse_pow = 0.0; rx_sig_pow =0.0; - for s=1:Nsymbols - arx_symb = rx((s-1)*Ts + (1:Ts)); - for m=1:M - r= sum(exp(-j*w(m)*(1:Ts)) .* arx_symb); - rx_pows(m)= r * conj(r); - rx_filt(s,m) = abs(r); - end - [tmp symbol_index] = max(rx_filt(s,:)); - rx_sig_pow = rx_sig_pow + rx_pows(symbol_index); - rx_pows(symbol_index)=[]; - rx_nse_pow = rx_nse_pow + sum(rx_pows)/(M-1); - rx_bits(bps*(s-1)+1:bps*s) = demapper(symbol_index,:); - end - - % using Rxpower = v^2 + sigmal^2 - rx_sig_pow = rx_sig_pow/Nsymbols; - rx_nse_pow = rx_nse_pow/Nsymbols; - v_est = sqrt(rx_sig_pow-rx_nse_pow); - SNRest = rx_sig_pow/rx_nse_pow; - sig_est = sqrt(rx_nse_pow/2); % for Rayleigh: 2nd raw moment = 2 .sigma^2 - Kest = rx_sig_pow/(2.0*sig_est^2) -1.0; - - Nerrors = sum(xor(tx_bits, rx_bits)); - raw_ber = Nerrors/Nbits; - printf('Ec (dB): %4.1f M: %2d Total bits: %5d HD Errs: %6d (Raw) BER: %1.3f\n', ... - EcNodB, M, Nbits, Nerrors, raw_ber); - if plt==1, plot_hist(rx_filt,tx_symbols, M); end - -endfunction - -% simulate one codeword of 2 or 4 FSK -------------------------- -function [Nerrors raw_ber EcNodB] = run_single_ldpc(M, Ltype, Nbits,EbNodB, plt, HRA) - - disp([num2str(M) 'FSK coded test ... ']) - if M==2 - bps = 1; modulation = 'FSK'; mod_order=2; mapping = 'gray'; - elseif M==4 - bps = 2; modulation = 'FSK'; mod_order=4; mapping = 'gray'; - else - error('sorry - bad value of M!'); - end - decoder_type = 0; max_iterations = 100; - - Hsize=size(HRA); - Krate = (Hsize(2)-Hsize(1))/Hsize(2); % - EcNodB = EbNodB + 10*log10(Krate); - code_param = ldpc_init_user(HRA, modulation, mod_order, mapping); - Nframes = floor(Nbits/code_param.data_bits_per_frame); - Nbits = Nframes*code_param.data_bits_per_frame; - - % Encoder - data_bits = round(rand(1,code_param.data_bits_per_frame)); - tx_bits = []; - for f=1:Nframes; - codeword_bits = LdpcEncode(data_bits, code_param.H_rows, code_param.P_matrix); - tx_bits = [tx_bits codeword_bits]; - end - - % modem/channel simulation - [raw_ber rx_filt rx_bits tx_symbols demapper sig_est SNRlin v_est] = ... - run_single(tx_bits,M,EcNodB, plt ); - - % Decoder - Nerrors = 0; - for f=1:Nframes - st = (f-1)*code_param.coded_bits_per_frame/bps + 1; - en = st + code_param.coded_bits_per_frame/bps - 1; - if or(Ltype==1, Ltype==3) - if bps==1, - sd = rx_filt(st:en,1) - rx_filt(st:en,2); - % OR ind = rx_filt(st:en,1) > rx_filt(st:en,2); - % llr = ind'*2 -1; % this works but use SNR scaling - if Ltype==3, HD=1; else, HD = 0; end - llr = sd_to_llr(sd, HD)'; % David's orig 2FSK routine - end - if bps==2, - if Ltype==3, - llr = mfsk_hd_to_llrs(rx_filt(st:en,:), demapper); - else - error('Ltype =1 not provided for coded 4FSK'); - end - end - end - if Ltype==2, % SDs are converted to LLRs - % use the SD amp estimates, try Bill's new LLR routine - if plt==1, plot_pdf(v_est, sig_est); end - llr = mfsk_sd_to_llrs(rx_filt(st:en,:), demapper, v_est, sig_est); - end - if Ltype==4, - % use CML demod: non-coherent; still seems to need amplitude estimate - symL = DemodFSK(1/v_est*rx_filt(st:en,:)', SNRlin, 1); - llr = -Somap(symL); % now convert to bit LLRs - end - if plt==3, figure(204); hist(llr); - title('Histogram LLR for decoder inputs'); pause(0.01); end - - [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, ... - max_iterations, decoder_type, 1, 1); - Niters = sum(PCcnt~=0); - detected_data = x_hat(Niters,:); - Nerrors = Nerrors + sum(xor(data_bits, detected_data(1:code_param.data_bits_per_frame))); - end - - ber = Nerrors/Nbits; - printf('Eb (dB): %4.1f Data Bits: %4d Num CWs: %5d Tot Errs: %5d (Cod) BER: %1.3f\n', ... - EbNodB,Nbits , Nframes, Nerrors, ber); -endfunction - -function plot_hist(rx_filt,tx_symbols, M) - % more general version of previous fn; - plots histograms for any Tx patterns - Smax = 36; - X = 0:Smax-1; - H = zeros(1,Smax); H2 = zeros(1,Smax); s2=0.0; - for m = 1:M - ind = tx_symbols==m; - ind2 = tx_symbols~=m; - H= H+ hist(rx_filt(ind,m),X); - H2= H2+ hist(rx_filt(ind2,m),X); - x=rx_filt(ind2,m); - s2 =s2 + sum(x(:).^2)/length(x); - end - disp('noise RMS is '); sqrt(s2/4) - figure(207); clf, plot(X,H); - title([num2str(M) 'FSK pdf for rx=tx symbol']) - hold on, plot(X,H2,'g'); - title([num2str(M) 'FSK pdf for rx!=tx symbol']) - pause(0.1); -endfunction - - -% 2FSK SD -> LLR mapping that we used for Wenet SSTV system -function llr = sd_to_llr(sd, HD=0) % original 2FSK + HD option - sd = sd / mean(abs(sd)); - x = sd - sign(sd); - sumsq = sum(x.^2); - summ = sum(x); - mn = summ/length(sd); - estvar = sumsq/length(sd) - mn*mn; - estEsN0 = 1/(2* estvar + 1E-3); - if HD==0, - llr = 4 * estEsN0 * sd; - else - llr = 8 * estEsN0 * sign(sd); %%%% *4 >> *8 - endif -endfunction - - - -function llrs = mfsk_sd_to_llrs(SD, map, v, sig) - % Array of MFSK SD is converted to bit LLR vector according to mapping 'map' - % SD is M elements wide; map is M by log2(M); v and sig are params of Rician pdf - % Bill (VK5DSP) for Codec2, version 24/6/20 - - Llim = 1e-7; - XX =1.0; % check LLR scaling? - Smap = size(map); - Ssd = size(SD); - if Smap(1) == 2 - llrs = zeros(1, Ssd(1)); - bps = 1; - assert(2^bps == Ssd(2), 'wrong SD length?') - % note that when v=0, the pdf reverts to Rayleigh - for kk = 1: Ssd(1) - Prx_1 = rice(SD(kk,2), v, sig) * rice(SD(kk,1),0,sig); - Prx_0 = rice(SD(kk,2), 0, sig) * rice(SD(kk,1),v,sig); - - if or(isnan(Prx_1), Prx_1<Llim), Prx_1 = Llim; end - if or(isnan(Prx_0), Prx_0<Llim), Prx_0 = Llim; end - llrs(kk) = -XX * log(Prx_1/Prx_0); - %% note inversion of LLRs - endfor - else - if map==[0 0; 0 1; 1 0; 1 1] % a specific case for 4FSK - llrs = zeros(2, Ssd(1)); - for kk = 1: Ssd(1) - % pn_m is prob of sub car n given bit m is transmitted - p1_0 = rice(SD(kk,1), 0, sig); p1_1 = rice(SD(kk,1), v, sig); - p2_0 = rice(SD(kk,2), 0, sig); p2_1 = rice(SD(kk,2), v, sig); - p3_0 = rice(SD(kk,3), 0, sig); p3_1 = rice(SD(kk,3), v, sig); - p4_0 = rice(SD(kk,4), 0, sig); p4_1 = rice(SD(kk,4), v, sig); - - % second bit - Prx_1 = p1_0*p3_0*(p2_1*p4_0 + p2_0*p4_1); - Prx_0 = p2_0*p4_0*(p1_1*p3_0 + p1_0*p3_1); - - if or(isnan(Prx_1), Prx_1<Llim), Prx_1 = Llim; end - if or(isnan(Prx_0), Prx_0<Llim), Prx_0 = Llim; end - llrs(2,kk) = -XX*log(Prx_1/Prx_0); - - % 1st bit (LHS) - Prx_1 = p1_0*p2_0*(p3_1*p4_0 + p3_0*p4_1); - Prx_0 = p3_0*p4_0*(p1_1*p2_0 + p1_0*p2_1); - - if or(isnan(Prx_1), Prx_1<Llim), Prx_1 = Llim; end - if or(isnan(Prx_0), Prx_0<Llim), Prx_0 = Llim; end - llrs(1,kk) = -XX* log(Prx_1/Prx_0); - endfor - llrs= llrs(:); % convert to single vector of LLRs - llrs = llrs'; - else - disp('the mapping is '); map - error('not sure how that mapping works!!') - endif - endif -endfunction - -function llrs = mfsk_hd_to_llrs(sd, map, SNRlin=4) -% for M=4, compare these results to SD method of mfsk_sd_to_llrs - [x, ind] = max(sd'); % find biggest SD - b2 = map(ind', :)'; % get bit pairs - b1 = b2(:)'; % convert to row vector - b1 = b1*2 -1; % convert to -1/ +1 - llrs = -4*SNRlin*b1; % approx scaling; default is ~6dB SNR -endfunction - - -% main ------------------------------------------------------------------------------ - -format short -more off -init_cml(); - -if exist('Ctype')==0, Ctype=1, end - -if Ctype==1 - load H_256_768_22.txt; HRA = H_256_768_22; % rate 1/3 -elseif Ctype==2 - load H_256_512_4.mat; HRA=H; % rate 1/2 code -elseif Ctype==3 - load HRAa_1536_512.mat; % rate 3/4 2k code -else - error('bad Ctype'); -end - -Hsize=size(HRA); -printf('using LDPC code length: %5d: Code rate: %5.2f\n',length(HRA), (Hsize(2)-Hsize(1))/Hsize(2)) - -% store results in array "res" and plot afterwards -% retain prev sims? -if exist('keep_res')==0, nrun = 0; clear res; end - -if exist('Nbits')==0, Nbits = 20000, end -if exist('Evec')==0, Evec = [5: 0.5: 9], end -if exist('plt')==0. plt=0; end - -if exist('M')==0, M=2, end - -for Ltype = [ 2 4] % << add other Ltype options here, if required - for Eb = Evec - if and(M==4, Ltype==1) - disp('skip original SD >> LLRs in 4FSK') - else - [Nerr raw_ber Ec] = run_single_ldpc(M, Ltype, Nbits, Eb, plt, HRA); - nrun = nrun+1; res(nrun,:) = [Eb Ec M Ltype Nbits Nerr raw_ber]; - endif - end -end -disp('results stored in array "res" - plot with fsk_cml_plot') -disp(' Eb Ec M Ltype #Bits #Errs Raw-BER') -for n = 1: nrun - printf(' %5.1f %5.1f %3d %3d %6d %6d %5.4f \n', res(n,:)) -end - - diff --git a/octave/fsk_demod_BER_test.py b/octave/fsk_demod_BER_test.py deleted file mode 100755 index f783cf0..0000000 --- a/octave/fsk_demod_BER_test.py +++ /dev/null @@ -1,610 +0,0 @@ -#!/usr/bin/env python3 -# -# Perform automated Eb/N0 testing of the C-implementation of fsk_mod / fsk_demod -# -# Based on the analysis performed here: https://github.com/projecthorus/radiosonde_auto_rx/blob/master/auto_rx/test/notes/2019-03-03_generate_lowsnr_validation.md -# -# Copyright (C) 2020 Mark Jessop <[email protected]> -# Released under GNU GPL v3 or later -# -# Requirements: -# - csdr must be installed and available on the path. https://github.com/ha7ilm/csdr -# - The following utilities from codec2 need to be built: -# - fsk_get_test_bits, fsk_put_test_bits -# - fsk_mod, fsk_demod -# - Create the directories: 'samples' and 'generated' in this directory (octave) -# -import json -import logging -import os -import time -import traceback -import subprocess -import sys -import argparse - -import numpy as np -import matplotlib.pyplot as plt -import scipy.signal -import scipy.interpolate - - -# Variables you will want to adjust: - -# Eb/N0 Range to test: -# Default: 0 through 5 dB in 0.5 db steps, then up to 20 db in 1db steps. -EBNO_RANGE = np.append(np.arange(0, 5, 0.5), np.arange(5, 20.5, 1)) - -# Baud rates to test: -BAUD_RATES = [100, 50, 25] - -# Order of the FSK signal (2 or 4) -FSK_ORDER = 4 - -# Test Length (bits) -TEST_LENGTH = 2e4 - -# Pseudorandom sequence length to generate test frames. -# NOTE: BER results are quite dependent on the frame length and threshold parameters. -FRAME_LENGTH = 2000 - -# Frame threshold detection. This has the effect of setting an upper bound on the BER. -FRAME_THRESHOLD = 0.4 - -# Allow a reduction in 'expected' bits of this value, as we expect the modem to need -# some time to 'spin up' the estimators. -FRAME_IGNORE = FRAME_LENGTH - -# IF sample rate -SAMPLE_RATE = 48000 - -# Frequency estimator limits -ESTIMATOR_LOWER_LIMIT = 100 -ESTIMATOR_UPPER_LIMIT = int(SAMPLE_RATE/2 - 1000) - -# Frequency of the low tone (Hz) -LOW_TONE = 2000 - -# Tone spacing (Hz) -TONE_SPACING = 270 - -# Mask Estimator -MASK_ESTIMATOR = True - -# Switch to 'Low Bit-Rate' mode below this baud rate. -#LBR_BREAK_POINT = 600 # No more LBR mode - -# Halt simulation for a particular baud rate when the BER drops below this level. -BER_BREAK_POINT = 1e-4 - -# If enabled, calculate Frequency Estimator error -FEST_ERROR = True - -# Frequency estimator error calculation threshold (*Rs) -FEST_THRESHOLD = 0.2 - -# Enable doppler shift. -# NOTE: This will apply up to +/- 6kHz of doppler shift to the test signal, -# emulating transmission through a LEO linear transponder. -# You will need to set the modem centre frequency and parameters such that -# the modem signal will be contained within a 1-22 kHz modem RX passband. -# For 100 baud Horus Binary testing, I use: -# LOW_TONE = 10000 -# TONE_SPACING = 250 -# The TEST_LENGTH setting must also be long enough so that the test modem file -# is at least 780 seconds long. For 100 baud 4FSK, a TEST_LENGTH of 2e6 is enough. -DOPPLER_ENABLED = False -DOPPLER_FILE = "doppler.npz" # generate using sat_doppler.py - - -STATS_OUTPUT = True - -# Where to place the initial test samples. -SAMPLE_DIR = "./samples" - -# Where to place the generated low-SNR samples. -GENERATED_DIR = "./generated" - -# Location of the codec2 utils -CODEC2_UTILS = "../build/src" - - -THEORY_EBNO = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] -THEORY_BER_4 = [ - 0.22934, - 0.18475, - 0.13987, - 0.09772, - 0.06156, - 0.03395, - 0.01579, - 0.00591, - 0.00168, - 3.39e-4, - 4.44e-5, - # 3.38e-6, - # 1.30e-7, - # 2.16e-9, - # 1.23e-11, - # 1.85e-14, - # 5.13e-18, - # 1.71e-22 -] -THEORY_BER_2 = [ - 0.30327, - 0.26644, - 0.22637, - 0.18438, - 0.14240, - 0.10287, - 0.06831, - 0.04080, - 0.02132, - 0.00942, - 0.00337, -] - -# -# Functions to read files and add noise. -# - - -def load_sample(filename, loadreal=True): - # If loading real samples (which is what fsk_mod outputs), apply a hilbert transform to get an analytic signal. - if loadreal: - return scipy.signal.hilbert(np.fromfile(filename, dtype="f4")) - else: - return np.fromfile(filename, dtype="c8") - - -def save_sample(data, filename): - # We have to make sure to convert to complex64.. - data.astype(dtype="c8").tofile(filename) - - # TODO: Allow saving as complex s16 - see view solution here: https://stackoverflow.com/questions/47086134/how-to-convert-a-numpy-complex-array-to-a-two-element-float-array - - -def apply_doppler(data, dopplerfile, fs=48000): - """ Apply a doppler curve to an input data stream """ - - npzfile = np.load(dopplerfile) - - _time = npzfile["arr_0"] - _doppler = npzfile["arr_1"] - - if len(data) < _time[-1] * fs: - print("Input data length too short - use more bits!") - - # Clip data length if its too long for the input doppler data - if len(data) > _time[-1] * fs: - data = data[: int(_time[-1] * fs)] - - # Interpolate the doppler data - _interp = scipy.interpolate.interp1d(_time, _doppler, kind="cubic") - - _timesteps = np.arange(0, len(data)) / fs - _interp_doppler = _interp(_timesteps) - - # This isn't working properly. - phase = np.cumsum(_interp_doppler / fs) - mixed = data * np.exp(1.0j * 2.0 * np.pi * phase) - - return mixed - - -def calculate_variance(data, threshold=-100.0): - # Calculate the variance of a set of radiosonde samples. - # Optionally use a threshold to limit the sample the variance - # is calculated over to ones that actually have sonde packets in them. - - _data_log = 20 * np.log10(np.abs(data)) - - # MSE is better than variance as a power estimate, as it counts DC - data_thresh = data[_data_log > threshold] - return np.mean(data_thresh * np.conj(data_thresh)) - - -def add_noise( - data, - variance, - baud_rate, - ebno, - fs=96000, - bitspersymbol=1.0, - normalise=True, - real=False, -): - # Add calibrated noise to a sample. - - # Calculate Eb/No in linear units. - _ebno = 10.0 ** ((ebno) / 10.0) - - # Calculate the noise variance we need to add - _noise_variance = variance * fs / (baud_rate * _ebno * bitspersymbol) - - # If we are working with real samples, we need to halve the noise contribution. - if real: - _noise_variance = _noise_variance * 0.5 - - # Generate complex random samples - _rand_i = np.sqrt(_noise_variance / 2.0) * np.random.randn(len(data)) - _rand_q = np.sqrt(_noise_variance / 2.0) * np.random.randn(len(data)) - - _noisy = data + (_rand_i + 1j * _rand_q) - - if normalise: - # print("Normalised to 1.0") - return _noisy / np.max(np.abs(_noisy)) - else: - return _noisy - - -def generate_lowsnr(sample, outfile, fs, baud, ebno, order): - """ Generate a low SNR test file """ - - if order == 2: - _bits_per_symbol = 1 - else: - _bits_per_symbol = 2 - - _var = calculate_variance(sample) - - _noisy = add_noise(sample, _var, baud, ebno, fs, _bits_per_symbol) - - save_sample(_noisy, outfile) - - return outfile - - -# -# Functions to deal with codec2 utils -# - - -def generate_fsk(baud): - """ Generate a set of FSK data """ - - # Calculate the number of bits we need to generate to get out desired test length. - if FSK_ORDER == 2: - _num_bits = TEST_LENGTH * baud - else: - _num_bits = TEST_LENGTH * baud * 2 - - _num_bits = TEST_LENGTH - - _filename = "%s/fsk_%d_%d_%d_f.bin" % (SAMPLE_DIR, FSK_ORDER, SAMPLE_RATE, baud) - - # Generate the command we need to make: - _cmd = ( - "%s/fsk_get_test_bits - %d %d | %s/fsk_mod %d %d %d %d %d - - | csdr convert_s16_f > %s" - % ( - CODEC2_UTILS, - _num_bits, - FRAME_LENGTH, - CODEC2_UTILS, - FSK_ORDER, - SAMPLE_RATE, - baud, - LOW_TONE, - TONE_SPACING, - _filename, - ) - ) - - print(_cmd) - - print("Generating test signal: %d-FSK, %d baud" % (FSK_ORDER, baud)) - - # Run the command. - try: - _start = time.time() - _output = subprocess.check_output(_cmd, shell=True, stderr=None) - _output = _output.decode() - except: - # traceback.print_exc() - _output = "error" - - _runtime = time.time() - _start - - print("Finished generating test signal.") - - return _filename - - -def process_fsk( - filename, baud, complex_samples=True, override_bits=None, stats=False, statsfile="" -): - """ Run a fsk file through fsk_demod """ - - _estim_limits = "-b %d -u %d " % (ESTIMATOR_LOWER_LIMIT, ESTIMATOR_UPPER_LIMIT) - - if MASK_ESTIMATOR: - _mask = "--mask %d " % TONE_SPACING - else: - _mask = "" - - if complex_samples: - _cpx = "--cs16 " - else: - _cpx = "" - - if stats: - _stats_file = GENERATED_DIR + "/" + statsfile + ".stats" - _stats = "--stats=50 " - else: - _stats = "" - _stats_file = None - - _cmd = "cat %s | csdr convert_f_s16 | %s/fsk_demod %s%s%s%s%d %d %d - - " % ( - filename, - CODEC2_UTILS, - _mask, - _estim_limits, - _cpx, - _stats, - FSK_ORDER, - SAMPLE_RATE, - baud, - ) - - if stats: - _cmd += "2> %s" % _stats_file - - _cmd += "| %s/fsk_put_test_bits -f %d -t %.2f - 2>&1" % ( - CODEC2_UTILS, - FRAME_LENGTH, - FRAME_THRESHOLD, - ) - - # print("Processing %s" % filename) - - print(_cmd) - - # Run the command. - try: - _start = time.time() - _output = subprocess.check_output(_cmd, shell=True) - _output = _output.decode() - except subprocess.CalledProcessError as e: - _output = e.output.decode() - except: - traceback.print_exc() - _output = "error" - print("Run failed!") - return (-1, _stats_file) - - _runtime = time.time() - _start - - # Try to grab last line of the stderr output - try: - _last_line = _output.split("\n")[-3] - except: - # Lack of a line indicates that we have decoded no data. return a BER of 1. - print("No bits decoded.") - return (1.0, _stats_file) - - # Detect no decoded bits when feeding in custom put_bits parameters. - if "Using" in _last_line: - print("No bits decoded.") - return (1.0, _stats_file) - - # Example line: - # [0009] BER 0.000, bits tested 18000, bit errors 0 errs: 0 - # [0009] BER 0.000, bits tested 18000, bit errors 0 - # PASS - # - - # split into fields - _fields = _last_line.split() - - # Extract number of bits and errors - _bits = float(_fields[5][:-1]) # remove the trailing comma - _errors = float(_fields[8]) - - print("Bits: %d, Errors: %d, Raw BER: %.8f" % (_bits, _errors, _errors / _bits)) - - if override_bits != None: - if _bits < override_bits: - print("Demod got %d bits, but we sent %d bits." % (_bits, override_bits)) - _errors += override_bits - _bits - - # Calculate and return BER - _ber = _errors / _bits - - if _ber > 1.0: - _ber = 1.0 - - return (_ber, _stats_file) - - -def read_stats(filename, sps = 50): - """ Read in a statistics file, and re-organise it for easier calculations """ - - _output = { - 'ebno': [], - 'f1_est': [], - 'f2_est': [], - 'f3_est': [], - 'f4_est': [], - 'ppm': [], - 'time': [] - } - - with open(filename, 'r') as _f: - for _line in _f: - if _line[0] != '{': - continue - - try: - _data = json.loads(_line) - except Exception as e: - #print("Line parsing error: %s" % str(e)) - continue - - _output['ebno'].append(_data['EbNodB']) - _output['f1_est'].append(_data['f1_est']) - _output['f2_est'].append(_data['f2_est']) - - if 'f3_est' in _data: - _output['f3_est'].append(_data['f3_est']) - _output['f4_est'].append(_data['f4_est']) - - _output['ppm'].append(_data['ppm']) - - if _output['time'] == []: - _output['time'] = [0] - else: - _output['time'].append(_output['time'][-1]+1.0/sps) - - return _output - - -def freq_est_error(data, Rs): - """ Calculate the frequency estimator error """ - - _threshold = FEST_THRESHOLD*Rs - - _total_points = len(data['f1_est'])*FSK_ORDER - - _errors = 0 - - _errors += np.sum(np.abs(np.array(data['f1_est'])-LOW_TONE) > _threshold) - _errors += np.sum(np.abs(np.array(data['f2_est'])-LOW_TONE-TONE_SPACING) > _threshold) - - if FSK_ORDER == 4: - _errors += np.sum(np.abs(np.array(data['f3_est'])-LOW_TONE-TONE_SPACING*2) > _threshold) - _errors += np.sum(np.abs(np.array(data['f4_est'])-LOW_TONE-TONE_SPACING*3) > _threshold) - - - return _errors/_total_points - - -if __name__ == "__main__": - - parser = argparse.ArgumentParser(description="FSK modem BER simulations") - parser.add_argument("--test", action="store_true", help="run automated test") - args = parser.parse_args() - - if args.test: - # test the AWGN channel simulation. We use BPSK that's phase shifted to exercise the - # complex maths a bit - - nb_bits = 100000 - EbNo = 4 - tx_bits = np.random.randint(2, size=nb_bits) - tx_bpsk_symbols = (2 * tx_bits - 1) * np.exp(1j * np.pi / 3) - tx_power = calculate_variance(tx_bpsk_symbols) - - # check calculate_variance() - assert tx_power < 1.1 and tx_power > 0.9 - - # BPSK modem simulation - rx_bpsk_symbols = add_noise(tx_bpsk_symbols, tx_power, 1, EbNo, 1, 1) - rx_bpsk_symbols = rx_bpsk_symbols * np.exp(-1j * np.pi / 3) - rx_bits = np.array([1 if np.real(s) > 0 else 0 for s in rx_bpsk_symbols]) - nb_errors = np.sum(np.bitwise_xor(tx_bits, rx_bits)) - ber = nb_errors / nb_bits - - # set limit of +/- 0.25dB on BER - EbNo_lin_upper = 10 ** ((EbNo - 0.25) / 10) - bpsk_ber_theory_upper = 0.5 * scipy.special.erfc(np.sqrt(EbNo_lin_upper)) - EbNo_lin_lower = 10 ** ((EbNo + 0.25) / 10) - bpsk_ber_theory_lower = 0.5 * scipy.special.erfc(np.sqrt(EbNo_lin_lower)) - print( - "nb_errors: %d ber: %4.3f ber_lower_limit: %4.3f ber_upper_limit: %4.3f" - % (nb_errors, ber, bpsk_ber_theory_lower, bpsk_ber_theory_upper) - ) - assert ber < bpsk_ber_theory_upper and ber > bpsk_ber_theory_lower - print("AWGN channel simulation test PASSED!") - exit() - - plot_data = {} - - for _baud in BAUD_RATES: - - _file = generate_fsk(_baud) - - print("Loading file and converting to complex.") - _sample = load_sample(_file) - - if DOPPLER_ENABLED: - print("Applying Doppler.") - _sample = apply_doppler(_sample, DOPPLER_FILE) - print("Done.") - _override_bits = _baud * (len(_sample) / SAMPLE_RATE) - FRAME_IGNORE - else: - _override_bits = TEST_LENGTH - FRAME_IGNORE - - _temp_file = "%s/temp.bin" % GENERATED_DIR - - _ebnos = [] - _bers = [] - _fest_err = [] - - for _ebno in EBNO_RANGE: - generate_lowsnr(_sample, _temp_file, SAMPLE_RATE, _baud, _ebno, FSK_ORDER) - - _ber, _stats_file = process_fsk( - _temp_file, - _baud, - override_bits=_override_bits, - stats=STATS_OUTPUT, - statsfile="fsk_%d_%.1f" % (_baud, _ebno), - ) - - print("%.1f, %.8f" % (_ebno, _ber)) - - _ebnos.append(_ebno) - _bers.append(_ber) - - if FEST_ERROR: - _stats = read_stats(_stats_file) - _fest_err.append(freq_est_error(_stats, _baud)) - - # Halt the simulation if the BER drops below our break point. - if _ber < BER_BREAK_POINT: - break - - plot_data[_baud] = {"baud": _baud, "ebno": _ebnos, "ber": _bers, "fest_err":_fest_err} - print(plot_data[_baud]) - - # plt.semilogy(plot_data[_baud]['ebno'], plot_data[_baud]['ber'], label="Simulated - %d bd" % _baud) - - plt.figure() - - print(plot_data) - - for _b in plot_data: - plt.semilogy( - plot_data[_b]["ebno"], plot_data[_b]["ber"], label="Simulated - %d bd" % _b - ) - - if FSK_ORDER == 2: - plt.semilogy(THEORY_EBNO, THEORY_BER_2, label="Theory") - else: - plt.semilogy(THEORY_EBNO, THEORY_BER_4, label="Theory") - - plt.xlabel("Eb/N0 (dB)") - plt.ylabel("BER") - - # Crop plot to reasonable limits - plt.ylim(1e-3, 1) - plt.xlim(0, 10) - - plt.title("fsk_demod %d-FSK BER performance" % FSK_ORDER) - plt.grid() - plt.legend() - - if FEST_ERROR: - plt.figure() - - for _b in plot_data: - plt.plot( - plot_data[_b]["ebno"], plot_data[_b]["fest_err"], label="Simulated - %d bd" % _b - ) - - plt.title("Frequency Estimator Error") - plt.grid() - plt.legend() - - plt.show() diff --git a/octave/fsk_llr_plot.m b/octave/fsk_llr_plot.m deleted file mode 100644 index f1b1a9f..0000000 --- a/octave/fsk_llr_plot.m +++ /dev/null @@ -1,51 +0,0 @@ -% Plot some results from FSK LLR tests -% Assume array "res" contains rows of simulation results::: -% Eb Ec M Ltype Nbits Nerr BERraw -% (some uncoded rows might contain -1 to indicate val is not applicable) - -figure(102); clf; hold on; - -%uncoded results -sub = res(res(:,4)==-1 & res(:,3)==2, :) -semilogy(sub(:,1), sub(:,7), 'k+--') -sub = res(res(:,4)==-1 & res(:,3)==4, :) -semilogy(sub(:,1), sub(:,7), 'k--') - -leg=[]; -% coded results -for M = [2 4 ] - -if M==2, lt = '-+'; else lt='-x'; end - - sub = res(res(:,4)==1 & res(:,3)==M, :) - if length(sub)>0, - semilogy(sub(:,1), sub(:,6)./sub(:,5), ['k' lt]) - leg= [leg; 'Orig LLRs']; - end - - -sub = res(res(:,4)==2 & res(:,3)==M, :) -if length(sub)>0, - semilogy(sub(:,1), sub(:,6)./sub(:,5), ['g' lt]) - leg= [leg; ' PDF LLRs']; -end - -sub = res(res(:,4)==3 & res(:,3)==M, :) -if length(sub)>0 - semilogy(sub(:,1), sub(:,6)./sub(:,5), ['b' lt]) - leg= [leg; ' HD LLRs']; -end - -sub = res(res(:,4)==4 & res(:,3)==M, :) -semilogy(sub(:,1), sub(:,6)./sub(:,5), ['m' lt]) -leg= [leg; ' CML LLRs']; -endfor - -ylabel('BER') -xlabel('Eb/N0 (Info Bits; dB)') - title('MFSK LLR test (+is 2FSK, xis 4FSK') - legend(leg) - legend('boxoff'); - -if exist('plotname'), print -dpng plotname; disp('saved'); end - diff --git a/octave/fsk_llr_test.m b/octave/fsk_llr_test.m deleted file mode 100644 index 3c21c64..0000000 --- a/octave/fsk_llr_test.m +++ /dev/null @@ -1,294 +0,0 @@ -% fsk_llr_test.m -% -% 2/4FSK simulation to develop LLR estimation algorithms for 4FSK/LDPC modems -% Modified version of David's fsk_llr.m; Bill - -#{ - TODO - The 'v' param of the Ricean pdf is the signal-only amplitude: genie value=16 - In practice, given varying input levels, this value needs to be estimated. - - A small scaling factor seems to improve 2FSK performance -- probably the 'sig' - estimate can be improved. - - Only tested with short code -- try a longer one! - - Simulation should be updated to exit Eb after given Nerr reached - -#} - -ldpc; - -% define Rician pdf -function y = rice(x,v,s) - s2 = s*s; - y = (x / s2) .* exp(-0.5 * (x.^2 + v.^2)/s2) .* besseli(0, x*v/s2); -endfunction - -function plot_pdf(v,s) - x=(0:0.1:2*v); - y= rice(x, v, s); - figure(201); hold on - plot(x,y,'g'); - %title('Rician pdf: signal carrier') - y= rice(x, 0, s); - plot(x,y,'b'); - title('Rician pdf: signal and noise-only carriers') - pause(0.01); -endfunction - -% single Eb/No point simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -function [raw_ber rx_filt rx_bits tx_symbols demapper sig_est ] = run_single(tx_bits, M, EcNodB, plt=0) - % Ec/N0 is per channel bit - bps = log2(M); % bits per symbol - Ts = 16; % length of each symbol in samples - - if length(tx_bits)==1 - Nbits = tx_bits; - tx_bits = randi(2,1,Nbits)-1; % make random bits - endif - - Nbits = length(tx_bits); - Nsymbols = Nbits/log2(M); - tx_symbols = zeros(1,Nsymbols); - - mapper = bps:-1:1; - % look up table demapper from symbols to bits (hard decision) - demapper=zeros(M,bps); - for m=1:M - for b=1:bps - if bitand(m-1,b) demapper(m,bps-b+1) = 1; end - end - end - - % continuous phase mFSK modulator - - w(1:M) = 2*pi*(1:M)/Ts; - tx_phase = 0; - tx = zeros(1,Ts*Nsymbols); - - for s=1:Nsymbols - bits_for_this_symbol = tx_bits(bps*(s-1)+1:bps*s); - symbol_index = bits_for_this_symbol * mapper' + 1; - tx_symbols(s) = symbol_index; - assert(demapper(symbol_index,:) == bits_for_this_symbol); - for k=1:Ts - tx_phase += w(symbol_index); - tx((s-1)*Ts+k) = exp(j*tx_phase); - end - end - - % AWGN channel noise - - - EsNodB = EcNodB + 10*log10(bps) - EsNo = 10^(EsNodB/10); - variance = Ts/EsNo; - noise = sqrt(variance/2)*(randn(1,Nsymbols*Ts) + j*randn(1,Nsymbols*Ts)); - rx = tx + noise; - - if plt==2, % check the Spectrum - [psd,Fpsd] =pwelch(rx,128,0.5,128,Ts); - figure(110); plot(Fpsd,10*log10(psd)); - title('Rx Signal: PSD '); - xlabel('Freq/Rs'); - %figure(111);plot(unwrap(arg(tx))); - pause(0.01); - endif - - - % integrate and dump demodulator - - rx_bits = zeros(1,Nbits); - rx_filt = zeros(Nsymbols,M); - rx_pows = zeros(1,M); - rx_nse_pow = 0.0; rx_sig_pow =0.0; - for s=1:Nsymbols - arx_symb = rx((s-1)*Ts + (1:Ts)); - for m=1:M - r= sum(exp(-j*w(m)*(1:Ts)) .* arx_symb); - rx_pows(m)= r * conj(r); - rx_filt(s,m) = abs(r); - end - [tmp symbol_index] = max(rx_filt(s,:)); - rx_sig_pow = rx_sig_pow + rx_pows(symbol_index); - rx_pows(symbol_index)=[]; - rx_nse_pow = rx_nse_pow + sum(rx_pows)/(M-1); - rx_bits(bps*(s-1)+1:bps*s) = demapper(symbol_index,:); - end - % using Rxpower = v^2 + sigmal^2 - - rx_sig_pow = rx_sig_pow/Nsymbols; - rx_nse_pow = rx_nse_pow/Nsymbols; - sig_est = sqrt(rx_nse_pow/2) % for Rayleigh: 2nd raw moment = 2 .sigma^2 - Kest = rx_sig_pow/(2.0*sig_est^2) -1.0 - - Nerrors = sum(xor(tx_bits, rx_bits)); - raw_ber = Nerrors/Nbits; - printf("EcNodB: %4.1f M: %2d Uncoded Nbits: %5d Nerrors: %4d (Raw) BER: %1.3f\n", ... - EcNodB, M, Nbits, Nerrors, raw_ber); - if plt==1, plot_hist(rx_filt,tx_symbols, M); end - -endfunction - - -% Plot histograms of Rx filter outputs -function plot_hist(rx_filt,tx_symbols, M) - % more general version of previous fn; - plots histograms for any Tx patterns - Smax = 36; - X = 0:Smax-1; - H = zeros(1,Smax); H2 = zeros(1,Smax); s2=0.0; - for m = 1:M - ind = tx_symbols==m; - ind2 = tx_symbols~=m; - H= H+ hist(rx_filt(ind,m),X); - H2= H2+ hist(rx_filt(ind2,m),X); - x=rx_filt(ind2,m); - s2 =s2 + sum(x(:).^2)/length(x); - end - disp('noise RMS is '); sqrt(s2/4) - figure(1); clf; plot(X,H); - title([num2str(M) 'FSK pdf for rx=tx symbol']) - figure(2); clf; plot(X,H2); - title([num2str(M) 'FSK pdf for rx!=tx symbol']) - pause(0.1); - -endfunction - -% 2FSK SD -> LLR mapping that we used for Wenet SSTV system -function llr = sd_to_llr(sd, HD=0) % original 2FSK + HD option - sd = sd / mean(abs(sd)); - x = sd - sign(sd); - sumsq = sum(x.^2); - summ = sum(x); - mn = summ/length(sd); - estvar = sumsq/length(sd) - mn*mn; - estEsN0 = 1/(2* estvar + 1E-3); - if HD==0, - llr = 4 * estEsN0 * sd; - else - llr = 4 * estEsN0 * sign(sd); - endif -endfunction - - -% single point LDPC encoded frame simulation, using 2FSK as a tractable starting point -% Note: ~/cml/matCreateConstellation.m has some support for FSK - can it do 4FSK? - -%%%%%%%%%%%%%%%%%%%%%%%%% - function [Nerrors raw_ber EcNodB] = run_single_ldpc(M, Ltype, Nbits,EbNodB, plt=0) - - disp([num2str(M) 'FSK coded test ... ']) - if M==2 - bps = 1; modulation = 'FSK'; mod_order=2; mapping = 'gray'; - elseif M==4 - bps = 2; modulation = 'FSK'; mod_order=4; mapping = 'gray'; - else - error('sorry - bad value of M!'); - endif - decoder_type = 0; max_iterations = 100; - - load H_256_768_22.txt - Krate = 1/3; - EcNodB = EbNodB + 10*log10(Krate); - code_param = ldpc_init_user(H_256_768_22, modulation, mod_order, mapping); - Nframes = floor(Nbits/code_param.data_bits_per_frame) - Nbits = Nframes*code_param.data_bits_per_frame - - % Encoder - data_bits = round(rand(1,code_param.data_bits_per_frame)); - tx_bits = []; - for f=1:Nframes; - codeword_bits = LdpcEncode(data_bits, code_param.H_rows, code_param.P_matrix); - tx_bits = [tx_bits codeword_bits]; - end - %tx_bits = zeros(1,length(tx_bits)); - - % modem/channel simulation - [raw_ber rx_filt rx_bits tx_symbols demapper sig_est ] = run_single(tx_bits,M,EcNodB, 0 ); - - % Decoder - Nerrors = 0; - for f=1:Nframes - st = (f-1)*code_param.coded_bits_per_frame/bps + 1; - en = st + code_param.coded_bits_per_frame/bps - 1; - - if or(Ltype==1, Ltype==3) - if bps==1, - sd = rx_filt(st:en,1) - rx_filt(st:en,2); - % OR ind = rx_filt(st:en,1) > rx_filt(st:en,2); - % llr = ind'*2 -1; % this works but use SNR scaling - if Ltype==3, HD=1; else, HD = 0; endif - llr = sd_to_llr(sd, HD)'; - endif - if bps==2, - if Ltype==3, - llr = mfsk_hd_to_llrs(rx_filt(st:en,:), demapper); - else - error('Ltype =1 not provided for coded 4FSK'); - endif - endif - endif - if Ltype==2, % SDs are converted to LLRs - v=16; - if plt==1, plot_pdf(v, sig_est); endif - llr = mfsk_sd_to_llrs(rx_filt(st:en,:), demapper, v, sig_est); - endif - - [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, ... - max_iterations, decoder_type, 1, 1); - Niters = sum(PCcnt!=0); - detected_data = x_hat(Niters,:); - Nerrors += sum(xor(data_bits, detected_data(1:code_param.data_bits_per_frame))); - endfor - ber = Nerrors/Nbits; - printf("EbNodB: %4.1f Coded Nbits: %5d Nerrors: %4d BER: %1.3f\n", EbNodB, Nbits, Nerrors, ber); -endfunction - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%rand('seed',1); -%randn('seed',1); -format short -more off -init_cml(); - -% store results in array "res" and plot afterwards -% comment the following line if you want to retain prev sims -nrun = 0; clear res; - -Nbits = 20000; plt=0; - -#{ -disp(' uncoded runs') -for M= [2 4] -for Eb = [6:10] - raw_ber = run_single(Nbits, M, Eb, plt) % 2fsk coded - nrun = nrun+1; res(nrun,:) = [Eb Eb M -1 Nbits -1 raw_ber] -endfor -endfor -#} - -disp(' coded runs '); - -M=2, -for Ltype = [1 2 3] -for Eb = [7: 0.5: 9] - [Nerr raw_ber Ec] = run_single_ldpc(M, Ltype, Nbits, Eb, plt) - nrun = nrun+1; res(nrun,:) = [Eb Ec M Ltype Nbits Nerr raw_ber] -endfor -endfor - -M=4, %v=16; -for Ltype = [2 3] -for Eb = [8.0 8.3 8.6 ] - [Nerr raw_ber Ec] = run_single_ldpc(M, Ltype, Nbits, Eb, plt) - nrun = nrun+1; res(nrun,:) = [Eb Ec M Ltype Nbits Nerr raw_ber] -endfor -endfor - - -date = datestr(now) -save 'mfsk_test_res.mat' res date - diff --git a/octave/fsk_lock_down.m b/octave/fsk_lock_down.m deleted file mode 100644 index 85ac034..0000000 --- a/octave/fsk_lock_down.m +++ /dev/null @@ -1,393 +0,0 @@ -% fsk_lock_down.m -% David Rowe April 2020 -% -% tests for "Lock down" waveform, low Eb/No 4FSK - -fsk_lib; - -% "Genie" (ie perfect) timing estimate, we just want to see impact on BER from frequency est errors -function rx_bits = simple_fsk_demod(states, rx, f) - M = states.M; Ts = states.Ts; Fs = states.Fs; N = states.nin; - - nsymb = states.nin/states.Ts; - rx_filter = zeros(states.nin,M); - rx_symbols = zeros(nsymb,M); - - % Down convert each tone. We can use any time index for down - % conversion as it's non-coherent - for m=1:M - phase = exp(-j*2*pi*(1:N)'*f(m)/Fs); - rx_filter(:,m) = rx .* phase; - end - - % sum energy for each symbol - for s=1:nsymb - st_sym = (s-1)*Ts+1; en_sym = s*Ts; - for m=1:M - rx_symbols(s,m) = sum(rx_filter(st_sym:en_sym,m)); - end - end - - % map symbols back to bits - rx_bits = []; - for s=1:nsymb - [tone_max tone_index] = max(rx_symbols(s,:)); - arx_bits = dec2bin(tone_index - 1, states.bitspersymbol) - '0'; - rx_bits = [rx_bits arx_bits]; - end - -end - -% set up "lock down" waveform -function [states M bits_per_frame] = lock_down_init(Rs,Fs,df,tx_tone_separation=250) - M = 4; - states = fsk_init(Fs,Rs,M,P=8,nsym=100); - bits_per_frame = 512; - states.tx_real = 0; % complex signal - states.tx_tone_separation = tx_tone_separation; - states.ftx = -2.5*states.tx_tone_separation + states.tx_tone_separation*(1:M); - states.fest_fmin = -Fs/2; - states.fest_fmax = +Fs/2; - states.df = df; - - % cumulative PDF, cdf(x) probability of 0....x errors in frame - cdf = binocdf(1:bits_per_frame, bits_per_frame, 0.3); - % our valid frame threshold is 50% probability, so if we get this many errors - % we have a 50% chance it's a valid frame - nerrs_valid = find(cdf>=0.5)(1); - % our invalid frame threshold is 99% probability, so very unlikely to - % get this many errors - nerrs_invalid = find(cdf>=0.99)(1); - states.ber_valid_thresh = nerrs_valid/bits_per_frame; - states.ber_invalid_thresh = nerrs_invalid/bits_per_frame; -end - -% run a test at an Eb/No point, measure how many dud freq estimates using both algorithms -function [states f_log f_log2 num_dud1 num_dud2 ber ber2] = freq_run_test(EbNodB = 10, num_frames=10, Fs=8000, Rs=100, df=0) - [states M bits_per_frame] = lock_down_init(Rs,Fs,df); - N = states.N; - - EbNo = 10^(EbNodB/10); - variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol); - - nbits = bits_per_frame*num_frames; - tx_bits = round(rand(1,nbits)); - tx = fsk_mod(states, tx_bits); - noise = sqrt(variance/2)*randn(length(tx),1) + j*sqrt(variance/2)*randn(length(tx),1); - rx = tx + noise; - run_frames = floor(length(rx)/N)-1; - st = 1; f_log = []; f_log2 = []; rx_bits = []; rx_bits2 = []; - for f=1:run_frames - - % extract nin samples from input stream - nin = states.nin; - en = st + states.nin - 1; - - % due to nin variations it's possible to overrun buffer - if en < length(rx) - sf = rx(st:en); - states = est_freq(states, sf, states.M); - arx_bits = simple_fsk_demod(states, sf, states.f); - rx_bits = [rx_bits arx_bits]; - arx_bits = simple_fsk_demod(states, sf, states.f2); - rx_bits2 = [rx_bits2 arx_bits]; - f_log = [f_log; states.f]; f_log2 = [f_log2; states.f2]; - st += nin; - end - end - - % ignore start up transient - startup = 1; % TODO make this sensible/proportional so its scales across Rs - if num_frames > startup - tx_bits = tx_bits(startup*bits_per_frame:end); - rx_bits = rx_bits(startup*bits_per_frame:end); - rx_bits2 = rx_bits2(startup*bits_per_frame:end); - end - - % measure BER - nerrors = sum(xor(tx_bits(1:length(rx_bits)),rx_bits)); ber = nerrors/nbits; - nerrors2 = sum(xor(tx_bits(1:length(rx_bits2)),rx_bits2)); ber2 = nerrors2/nbits; - - % Lets say that for a valid freq estimate, all four tones must be within 0.1*Rs of their tx frequency - num_dud1 = 0; num_dud2 = 0; - for i=1:length(f_log) - if sum(abs(f_log(i,:)-states.ftx) > 0.1*states.Rs) - num_dud1++; - end - if sum(abs(f_log2(i,:)-states.ftx) > 0.1*states.Rs) - num_dud2++; - end - end -end - -function freq_run_single(EbNodB = 3, num_frames = 10) - [states f_log f_log2 num_dud1 num_dud2 ber ber2] = freq_run_test(EbNodB, num_frames); - - percent_dud1 = 100*num_dud1/length(f_log); - percent_dud2 = 100*num_dud2/length(f_log); - printf("EbNodB: %4.2f dB tests: %3d duds1: %3d %5.2f %% duds2: %3d %5.2f %% ber1: %4.3f ber2: %4.3f\n", - EbNodB, length(f_log), num_dud1, percent_dud1, num_dud2, percent_dud2, ber, ber2) - - figure(1); clf; - ideal=ones(length(f_log),1)*states.ftx; - plot((1:length(f_log)),ideal(:,1),'bk;ideal;') - hold on; plot((1:length(f_log)),ideal(:,2:states.M),'bk'); hold off; - hold on; - plot(f_log(:,1), 'linewidth', 2, 'b;peak;'); - plot(f_log(:,2:states.M), 'linewidth', 2, 'b'); - plot(f_log2(:,1),'linewidth', 2, 'r;mask;'); - plot(f_log2(:,2:states.M),'linewidth', 2, 'r'); - hold off; - xlabel('Time (frames)'); ylabel('Frequency (Hz)'); - title(sprintf("EbNo = %4.2f dB", EbNodB)); - print("fsk_freq_est_single.png", "-dpng") - - figure(2); clf; - errors = (f_log - states.ftx)(:); - ind = find(abs(errors) < 100); - errors2 = (f_log2 - states.ftx)(:); - ind2 = find(abs(errors2) < 100); - if length(ind) - subplot(211); hist(errors(ind),50) - end - if length(ind2) - subplot(212); hist(errors2(ind2),50) - end -end - - -% test peak and mask algorithms side by side -function freq_run_curve_peak_mask - - EbNodB = 0:9; - m4fsk_ber_theory = [0.23 0.18 0.14 0.09772 0.06156 0.03395 0.01579 0.00591 0.00168 3.39E-4]; - percent_log = []; ber_log = []; - for ne = 1:length(EbNodB) - [states f_log f_log2 num_dud1 num_dud2 ber ber2] = freq_run_test(EbNodB(ne), 10); - percent_dud1 = 100*num_dud1/length(f_log); - percent_dud2 = 100*num_dud2/length(f_log); - percent_log = [percent_log; [percent_dud1 percent_dud2]]; - ber_log = [ber_log; [ber ber2]]; - printf("EbNodB: %4.2f dB tests: %3d duds1: %3d %5.2f %% duds2: %3d %5.2f %% ber1: %4.3f ber2: %4.3f\n", - EbNodB(ne), length(f_log), num_dud1, percent_dud1, num_dud2, percent_dud2, ber, ber2) - end - - figure(1); clf; plot(EbNodB, percent_log(:,1), 'linewidth', 2, '+-;peak;'); grid; - hold on; plot(EbNodB, percent_log(:,2), 'linewidth', 2, 'r+-;mask;'); hold off; - xlabel('Eb/No (dB)'); ylabel('% Errors'); - title(sprintf("Fs = %d Rs = %d df = %3.2f", states.Fs, states.Rs, states.df)); - print("fsk_freq_est_errors.png", "-dpng") - - figure(2); clf; semilogy(EbNodB, m4fsk_ber_theory, 'linewidth', 2, 'bk+-;theory;'); grid; - hold on; semilogy(EbNodB, ber_log(:,1), 'linewidth', 2, '+-;peak;'); - semilogy(EbNodB, ber_log(:,2), 'linewidth', 2, 'r+-;mask;'); hold off; - xlabel('Eb/No (dB)'); ylabel('BER'); - title(sprintf("Fs = %d Rs = %d df = %3.2f", states.Fs, states.Rs, states.df)); - print("fsk_freq_est_ber.png", "-dpng") -end - - -function freq_run_curve_mask(Fs,Rs) - EbNodB = 0:9; - m4fsk_ber_theory = [0.23 0.18 0.14 0.09772 0.06156 0.03395 0.01579 0.00591 0.00168 3.39E-4]; - figure(1); clf; semilogy(EbNodB, m4fsk_ber_theory, 'linewidth', 2, 'bk+-;theory;'); grid; - xlabel('Eb/No (dB)'); ylabel('BER'); - title(sprintf("Mask: Fs = %d Hz Rs = %d Hz", Fs, Rs)); - hold on; - - for df=-0.01:0.01:0.01 - ber_log = []; - for ne = 1:length(EbNodB) - [states f_log f_log2 num_dud1 num_dud2 ber ber2] = freq_run_test(EbNodB(ne), 100, Fs, Rs, df*Rs); - ber_log = [ber_log; [ber ber2]]; - printf("Fs: %d Rs: %d df %3.2f EbNodB: %4.2f dB tests: %3d ber: %4.3f\n", - Fs, Rs, df, EbNodB(ne), length(f_log), ber2) - end - semilogy(EbNodB, ber_log(:,2), 'linewidth', 2, sprintf("+-;df=% 3.2f Hz/s;",df*Rs)); - end - hold off; - print(sprintf("fsk_freq_est_ber_%d_%d.png",Fs,Rs), "-dpng") -end - - -% Run a complete modem (freq and timing estimators running) at a -% single Eb/No point. At low Eb/No the estimators occasionally fall -% over so we get complete junk, we consider that case a packet error -% and exclude it from the BER estimation. - -function [states ber per] = modem_run_test(EbNodB = 10, num_frames=10, Fs=8000, Rs=100, df=0, plots=0, spreadHz=0,tx_tone_separation=250) - [states M bits_per_frame] = lock_down_init(Rs, Fs, df, tx_tone_separation); - N = states.N; - if plots; states.verbose = 0x4; end - EbNo = 10^(EbNodB/10); - variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol); - - nbits = bits_per_frame*num_frames; - test_frame = round(rand(1,bits_per_frame)); tx_bits = []; - for f=1:num_frames - tx_bits = [tx_bits test_frame]; - end - - tx = fsk_mod(states, tx_bits); - noise = sqrt(variance/2)*randn(length(tx),1) + j*sqrt(variance/2)*randn(length(tx),1); - if spreadHz - % just use phase part of doppler spread, not interested in amplitude fading - spread = doppler_spread(spreadHz, Fs, round(1.1*length(tx))); - spread = exp(j*arg(spread(1:length(tx)))); - rx = tx.*rot90(spread) + noise; - else - rx = tx + noise; - end - run_frames = floor(length(rx)/N)-1; - st = 1; f_log = []; f_log2 = []; rx_bits = []; rx_bits2 = []; - for f=1:run_frames - - % extract nin samples from input stream - nin = states.nin; - en = st + states.nin - 1; - - % due to nin variations it's possible to overrun buffer - if en < length(rx) - sf = rx(st:en); - states = est_freq(states, sf, states.M); states.f = states.f2; - [arx_bits states] = fsk_demod(states, sf); - rx_bits = [rx_bits arx_bits]; - f_log = [f_log; states.f]; - st += nin; - end - end - - num_frames=floor(length(rx_bits)/bits_per_frame); - log_nerrs = []; num_frames_rx = 0; - for f=1:num_frames-1 - st = (f-1)*bits_per_frame + 1; en = (f+1)*bits_per_frame; - states = ber_counter(states, test_frame, rx_bits(st:en)); - log_nerrs = [log_nerrs states.nerr]; - if states.ber_state; num_frames_rx++; end - end - if states.Terrs - printf("Fs: %d Rs: %d df % 3.2f sp: %2.1f EbNo: %4.2f ftx: %3d frx: %3d nbits: %4d nerrs: %3d ber: %4.3f\n", - Fs, Rs, df, spreadHz, EbNodB, num_frames, num_frames_rx, states.Tbits, states.Terrs, states.Terrs/states.Tbits); - ber = states.Terrs/states.Tbits; - else - ber = 0.5; - end - - if plots - figure(1); clf; - ideal=ones(length(f_log),1)*states.ftx; - plot((1:length(f_log)),ideal(:,1),'bk;ideal;') - hold on; plot((1:length(f_log)),ideal(:,2:states.M),'bk'); hold off; - hold on; - plot(f_log(:,1), 'linewidth', 2, 'b;peak;'); - plot(f_log(:,2:states.M), 'linewidth', 2, 'b'); - hold off; - xlabel('Time (frames)'); ylabel('Frequency (Hz)'); - figure(2); clf; plot(log_nerrs); title('Errors per frame'); - end - - per = 1 - num_frames_rx/num_frames; -end - - -% run BER v Eb/No curves over a range of frequency rate/change -function modem_run_curve(Fs, Rs, num_frames=100, dfmax=0.01) - EbNodB = 0:9; - m4fsk_ber_theory = [0.23 0.18 0.14 0.09772 0.06156 0.03395 0.01579 0.00591 0.00168 3.39E-4]; - figure(1); clf; semilogy(EbNodB, m4fsk_ber_theory, 'linewidth', 2, 'bk+-;theory;'); grid; - xlabel('Eb/No (dB)'); ylabel('BER'); - title(sprintf("Mask: Fs = %d Hz Rs = %d Hz", Fs, Rs)); hold on; - figure(2); clf; - xlabel('Eb/No (dB)'); ylabel('PER'); title(sprintf("Mask: Fs = %d Hz Rs = %d Hz", Fs, Rs)); - grid; axis([min(EbNodB) max(EbNodB) 0 1]); hold on; - - for df=-dfmax:dfmax:dfmax - ber_log = []; per_log = []; - for ne = 1:length(EbNodB) - [states ber per] = modem_run_test(EbNodB(ne), num_frames, Fs, Rs, df*Rs); - ber_log = [ber_log; ber]; per_log = [per_log; per]; - end - figure(1); semilogy(EbNodB, ber_log, 'linewidth', 2, sprintf("+-;df=% 3.2f Hz/s;",df*Rs)); - figure(2); plot(EbNodB, per_log, 'linewidth', 2, sprintf("+-;df=% 3.2f Hz/s;",df*Rs)); - end - - figure(1); hold off; print(sprintf("fsk_modem_ber_%d_%d.png",Fs,Rs), "-dpng") - figure(2); hold off; print(sprintf("fsk_modem_per_%d_%d.png",Fs,Rs), "-dpng") -end - -% run BER v Eb/No curve with some phase noise spreading the energy of the tones in frequency -function modem_run_curve_spread(Fs, Rs, num_frames=100) - EbNodB = 0:9; - m4fsk_ber_theory = [0.23 0.18 0.14 0.09772 0.06156 0.03395 0.01579 0.00591 0.00168 3.39E-4]; - figure(1); clf; semilogy(EbNodB, m4fsk_ber_theory, 'linewidth', 2, 'bk+-;theory;'); grid; - xlabel('Eb/No (dB)'); ylabel('BER'); - title(sprintf("Spread: Fs = %d Hz Rs = %d Hz", Fs, Rs)); hold on; - figure(2); clf; - xlabel('Eb/No (dB)'); ylabel('PER'); - title(sprintf("Spread: Fs = %d Hz Rs = %d Hz", Fs, Rs)); - grid; axis([min(EbNodB) max(EbNodB) 0 1]); hold on; - - spreadHz = [0.0 1 2 5]; - for ns = 1:length(spreadHz) - ber_log = []; per_log = []; - for ne = 1:length(EbNodB) - [states ber per] = modem_run_test(EbNodB(ne), num_frames, Fs, Rs, 0, 0, spreadHz(ns)); - ber_log = [ber_log; ber]; per_log = [per_log; per]; - end - figure(1); semilogy(EbNodB, ber_log, 'linewidth', 2, sprintf("+-;spread=% 3.2f Hz;",spreadHz(ns))); - figure(2); plot(EbNodB, per_log, 'linewidth', 2, sprintf("+-;spread=% 3.2f Hz;",spreadHz(ns))); - end - - figure(1); hold off; print(sprintf("fsk_modem_ber_spread_%d_%d.png",Fs,Rs), "-dpng") - figure(2); hold off; print(sprintf("fsk_modem_per_spread_%d_%d.png",Fs,Rs), "-dpng") -end - -% study code rate versus Rs and MDS -function code_rate_table - packet_duration_sec = 20; - k = 256; - noise_figure = 1; - bits_per_symbol = 2; - noise_bandwidth = 3000; - - code_rate=[1 0.8 0.5 1/3]; - raw_ber=[2E-3 0.04 0.08 0.16]; - EbNodB_4fsk=[8 4.5 3.5 1.5]; - - printf("Code Rate | Raw BER | 4FSK Eb/No | n,k | Rs | SNR | MDS |\n"); - printf("| --- | --- | --- | --- | --- | --- | --- |\n"); - for i=1:length(code_rate) - n = k/code_rate(i); - Rb = n/packet_duration_sec; - Rs = Rb/bits_per_symbol; - snr = EbNodB_4fsk(i) + 10*log10(Rb/noise_bandwidth); - mds = EbNodB_4fsk(i) + 10*log10(Rb) + noise_figure - 174; - printf("%3.2f | %4.3f | %2.1f | %d,%d | %4.1f | %4.1f | %5.1f |\n", - code_rate(i), raw_ber(i), EbNodB_4fsk(i), n, k, Rs, snr, mds); - end -end - -graphics_toolkit("gnuplot"); -more off; - -% same results every time -rand('state',1); -randn('state',1); - -% freq estimator tests (choose one) -#freq_run_single(3,10) -#freq_run_curve_peak_mask -#freq_run_curve_mask(8000,100) -#freq_run_curve_mask(24000,25) -#freq_run_curve_mask(8000,25) - -% complete modem tests (choose one) -#modem_run_curve(24000,25,100) -#modem_run_curve(8000,100,50,0.05) -#modem_run_curve_spread(8000,25,50) -#modem_run_curve(8000,100,20) -modem_run_test(2, 20, 8000, 25, 0, 1, 0, 270); - -% just print a table of code rates -#code_rate_table - diff --git a/octave/fsk_v_afsk.m b/octave/fsk_v_afsk.m deleted file mode 100644 index a692c18..0000000 --- a/octave/fsk_v_afsk.m +++ /dev/null @@ -1,232 +0,0 @@ -% fsk.m -% David Rowe Nov 2014 - -% Ideal non-coherent FSK and AFSK-over-analog-FM simulation. Can draw -% Eb/No curves or run single point simulations - -rand('state',1); -randn('state',1); -graphics_toolkit ("gnuplot"); - -fm; - -function sim_out = fsk_ber_test(sim_in) - Fs = 96000; - fmark = sim_in.fmark; - fspace = sim_in.fspace; - Rs = sim_in.Rs; - Ts = Fs/Rs; - emphasis = 50E-6; - verbose = sim_in.verbose; - - nsym = sim_in.nsym; - nsam = nsym*Ts; - EbNodB = sim_in.EbNodB; - - fm = sim_in.fm; - - if fm - fm_states.pre_emp = 0; - 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 - - % simulate over a range of Eb/No values - - for ne = 1:length(EbNodB) - Nerrs = Terrs = Tbits = 0; - - % 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); - variance = Fs/(Rs*EbNo); - - % Modulator ------------------------------- - - tx_bits = round(rand(1, nsym)); - tx = zeros(1,nsam); - tx_phase = 0; - - for i=1:nsym - for k=1:Ts - if tx_bits(i) == 1 - tx_phase += 2*pi*fmark/Fs; - else - tx_phase += 2*pi*fspace/Fs; - end - tx_phase = tx_phase - floor(tx_phase/(2*pi))*2*pi; - tx((i-1)*Ts+k) = exp(j*tx_phase); - end - end - - % Optional AFSK over FM modulator - - if sim_in.fm - % FM mod takes real input; +/- 1 for correct deviation - tx = analog_fm_mod(fm_states, real(tx)); - end - - % Channel --------------------------------- - - % We use complex (single sided) channel simulation, as it's convenient - % for the FM simulation. - - noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); - rx = tx + noise; - if verbose > 1 - printf("EbNo: %f Eb: %f var No: %f EbNo (meas): %f\n", - EbNo, var(tx)*Ts/Fs, var(noise)/Fs, (var(tx)*Ts/Fs)/(var(noise)/Fs)); - end - save fsk tx_bits rx - - % Optional AFSK over FM demodulator - - if sim_in.fm - % scaling factor for convenience to match pure FSK - rx_bb = 2*analog_fm_demod(fm_states, rx); - else - rx_bb = rx; - end - - % Demodulator ----------------------------- - - % non-coherent FSK demod - - mark_dc = rx_bb .* exp(-j*(0:nsam-1)*2*pi*fmark/Fs); - space_dc = rx_bb .* exp(-j*(0:nsam-1)*2*pi*fspace/Fs); - - rx_bits = zeros(1, nsym); - for i=1:nsym - st = (i-1)*Ts+1; - en = st+Ts-1; - mark_int(i) = sum(mark_dc(st:en)); - space_int(i) = sum(space_dc(st:en)); - rx_bits(i) = abs(mark_int(i)) > abs(space_int(i)); - end - - if fm - d = fm_states.nsym_delay; - error_positions = xor(rx_bits(1+d:nsym), tx_bits(1:(nsym-d))); - else - error_positions = xor(rx_bits, tx_bits); - end - Nerrs = sum(error_positions); - Terrs += Nerrs; - Tbits += length(error_positions); - - TERvec(ne) = Terrs; - BERvec(ne) = Terrs/Tbits; - - if verbose > 1 - figure(2) - clf - Rx = 10*log10(abs(fft(rx))); - plot(Rx(1:Fs/2)); - axis([1 Fs/2 0 50]); - - figure(3) - clf; - subplot(211) - plot(real(rx_bb(1:Ts*20))) - subplot(212) - Rx_bb = 10*log10(abs(fft(rx_bb))); - plot(Rx_bb(1:3000)); - axis([1 3000 0 50]); - - figure(4); - subplot(211) - stem(abs(mark_int(1:100))); - subplot(212) - stem(abs(space_int(1:100))); - end - - if verbose - printf("EbNo (db): %3.2f Terrs: %d BER: %3.2f \n", aEbNodB, Terrs, Terrs/Tbits); - end - end - - sim_out.TERvec = TERvec; - sim_out.BERvec = BERvec; -endfunction - - -function run_fsk_curves - sim_in.fmark = 1200; - sim_in.fspace = 2200; - sim_in.Rs = 1200; - sim_in.nsym = 12000; - sim_in.EbNodB = 0:2:20; - sim_in.fm = 0; - sim_in.verbose = 1; - - EbNo = 10 .^ (sim_in.EbNodB/10); - fsk_theory.BERvec = 0.5*exp(-EbNo/2); % non-coherent BFSK demod - fsk_sim = fsk_ber_test(sim_in); - - sim_in.fm = 1; - fsk_fm_sim = fsk_ber_test(sim_in); - - % BER v Eb/No curves - - figure(1); - clf; - semilogy(sim_in.EbNodB, fsk_theory.BERvec,'r;FSK theory;') - hold on; - semilogy(sim_in.EbNodB, fsk_sim.BERvec,'g;FSK sim;') - semilogy(sim_in.EbNodB, 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]) - 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(sim_in.Rs/1); - figure(2); - clf; - semilogy(sim_in.EbNodB+RsOnB_dB, fsk_theory.BERvec,'r;FSK theory;') - hold on; - semilogy(sim_in.EbNodB+RsOnB_dB, fsk_sim.BERvec,'g;FSK sim;') - semilogy(sim_in.EbNodB+RsOnB_dB, fsk_fm_sim.BERvec,'b;FSK over FM 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=1200 bit/s and 1 Hz noise bandwidth (dB)"); - ylabel("Bit Error Rate (BER)") -end - -function run_fsk_single - sim_in.fmark = 1000; - sim_in.fspace = 2000; - sim_in.Rs = 1000; - sim_in.nsym = 2000; - sim_in.EbNodB = 7; - sim_in.fm = 0; - sim_in.verbose = 1; - - fsk_sim = fsk_ber_test(sim_in); -endfunction - -# choose one of these functions below - -run_fsk_curves -#run_fsk_single - diff --git a/octave/fskdemodgui.py b/octave/fskdemodgui.py deleted file mode 100644 index f1a4f29..0000000 --- a/octave/fskdemodgui.py +++ /dev/null @@ -1,220 +0,0 @@ -#!/usr/bin/env python3 -# -# fsk_demod Statistics GUI -# Accepts the stats output from fsk_demod on stdin, and plots it. -# -# Mark Jessop 2016-03-13 <[email protected]> -# -# NOTE: This is intended to be run on a 'live' stream of samples, and hence expects -# updates at about 10Hz. Anything faster will fill up the input queue and be discarded. -# -# Call using: -# <producer>| ./fsk_demod --cu8 -s --stats=100 2 $SDR_RATE $BAUD_RATE - - 2> >(python fskdemodgui.py --wide) | <consumer> -# -# Dependencies: -# * Python (written for 2.7, only tested recently on 3+) -# * numpy -# * pyqtgraph -# * PyQt5 (Or some Qt5 backend compatible with pyqtgraph) -# -# -import sys, time, json, argparse -from threading import Thread -try: - from pyqtgraph.Qt import QtGui, QtCore -except ImportError: - print("Could not import PyQt5 - is it installed?") - sys.exit(1) - -try: - import numpy as np -except ImportError: - print("Could not import numpy - is it installed?") - sys.exit(1) - -try: - import pyqtgraph as pg -except ImportError: - print("Could not import pyqtgraph - is it installed?") - sys.exit(1) - -try: - # Python 2 - from Queue import Queue -except ImportError: - # Python 3 - from queue import Queue - -parser = argparse.ArgumentParser() -parser.add_argument("--wide", action="store_true", default=False, help="Alternate wide arrangement of widgets, for placement at bottom of 4:3 screen.") -args = parser.parse_args() - -# Some settings... -update_rate = 2 # Hz -history_size = 100 # 10 seconds at 10Hz... -history_scale = np.linspace((-1*history_size+1)/float(update_rate),0,history_size) - -# Input queue -in_queue = Queue(1) # 1-element FIFO... - -win = pg.GraphicsWindow() -win.setWindowTitle('FSK Demodulator Modem Statistics') - - -# Plot objects -ebno_plot = win.addPlot(title="Eb/No") -ppm_plot = win.addPlot(title="Sample Clock Offset") -if args.wide == False: - win.nextRow() -else: - win.resize(1024,200) -fest_plot =pg.PlotItem() # win.addPlot(title="Tone Frequency Estimation") -eye_plot = win.addPlot(title="Eye Diagram") -# Disable auto-ranging on eye plot and fix axes for a big speedup... -spec_plot = win.addPlot(title="Spectrum") -spec_plot.setYRange(0,40) -spec_plot.setLabel('left','SNR (dB)') -spec_plot.setLabel('bottom','FFT Bin') -# Configure plot labels and scales. -ebno_plot.setLabel('left','Eb/No (dB)') -ebno_plot.setLabel('bottom','Time (seconds)') -ebno_plot.setYRange(0,30) -ppm_plot.setLabel('left','Clock Offset (ppm)') -ppm_plot.setLabel('bottom','Time (seconds)') -fest_plot.setLabel('left','Frequency (Hz)') -fest_plot.setLabel('bottom','Time (seconds)') -eye_plot.disableAutoRange() -eye_plot.setYRange(0,1) -eye_plot.setXRange(0,15) -eye_xr = 15 - -# Data arrays... -ebno_data = np.zeros(history_size) -ppm_data = np.zeros(history_size) -fest_data = np.zeros((4,history_size)) - -# Curve objects, so we can update them... -spec_curve = spec_plot.plot([0]) -ebno_curve = ebno_plot.plot(x=history_scale,y=ebno_data) -ppm_curve = ppm_plot.plot(x=history_scale,y=ppm_data) -fest1_curve = fest_plot.plot(x=history_scale,y=fest_data[0,:],pen='r') # f1 = Red -fest2_curve = fest_plot.plot(x=history_scale,y=fest_data[1,:],pen='g') # f2 = Blue -fest3_curve = fest_plot.plot(x=history_scale,y=fest_data[2,:],pen='b') # f3 = Greem -fest4_curve = fest_plot.plot(x=history_scale,y=fest_data[3,:],pen='m') # f4 = Magenta - -# Plot update function. Reads from queue, processes and updates plots. -def update_plots(): - global timeout,timeout_counter,eye_plot,ebno_curve, ppm_curve, fest1_curve, fest2_curve, ebno_data, ppm_data, fest_data, in_queue, eye_xr, spec_curve - - try: - if in_queue.empty(): - return - in_data = in_queue.get_nowait() - in_data = json.loads(in_data) - except Exception as e: - - sys.stderr.write(str(e)) - return - - # Roll data arrays - ebno_data[:-1] = ebno_data[1:] - ppm_data[:-1] = ppm_data[1:] - fest_data = np.roll(fest_data,-1,axis=1) - - - # Try reading in the new data points from the dictionary. - try: - new_ebno = in_data['EbNodB'] - new_ppm = in_data['ppm'] - new_fest1 = in_data['f1_est'] - new_fest2 = in_data['f2_est'] - new_spec = in_data['samp_fft'] - except Exception as e: - print("ERROR reading dict: %s" % e) - - # Try reading in the other 2 tones. - try: - new_fest3 = in_data['f3_est'] - new_fest4 = in_data['f4_est'] - fest_data[2,-1] = new_fest3 - fest_data[3,-1] = new_fest4 - except: - # If we can't read these tones out of the dict, fill with NaN - fest_data[2,-1] = np.nan - fest_data[3,-1] = np.nan - - # Add in new data points - ebno_data[-1] = new_ebno - ppm_data[-1] = new_ppm - fest_data[0,-1] = new_fest1 - fest_data[1,-1] = new_fest2 - - - # Update plots - spec_data_log = 20*np.log10(np.array(new_spec)+0.01) - spec_curve.setData(spec_data_log) - spec_plot.setYRange(spec_data_log.max()-50,spec_data_log.max()+10) - ebno_curve.setData(x=history_scale,y=ebno_data) - ppm_curve.setData(x=history_scale,y=ppm_data) - fest1_curve.setData(x=history_scale,y=fest_data[0,:],pen='r') # f1 = Red - fest2_curve.setData(x=history_scale,y=fest_data[1,:],pen='g') # f2 = Blue - fest3_curve.setData(x=history_scale,y=fest_data[2,:],pen='b') # f3 = Green - fest4_curve.setData(x=history_scale,y=fest_data[3,:],pen='m') # f4 = Magenta - - #Now try reading in and plotting the eye diagram - try: - eye_data = np.array(in_data['eye_diagram']) - - #eye_plot.disableAutoRange() - eye_plot.clear() - col_index = 0 - for line in eye_data: - eye_plot.plot(line,pen=(col_index,eye_data.shape[0])) - col_index += 1 - #eye_plot.autoRange() - - #Quick autoranging for x-axis to allow for differing P and Ts values - if eye_xr != len(eye_data[0]) - 1: - eye_xr = len(eye_data[0]) - 1 - eye_plot.setXRange(0,len(eye_data[0])-1) - - except Exception as e: - pass - - -timer = pg.QtCore.QTimer() -timer.timeout.connect(update_plots) -timer.start(1000/update_rate) - - -# Thread to read from stdin and push into a queue to be processed. -def read_input(): - global in_queue - - while True: - in_line = sys.stdin.readline() - if type(in_line) == bytes: - in_line = in_line.decode() - - # Only push actual data into the queue... - # This stops sending heaps of empty strings into the queue when fsk_demod closes. - if in_line == "": - time.sleep(0.1) - continue - - if not in_queue.full(): - in_queue.put_nowait(in_line) - - -read_thread = Thread(target=read_input) -read_thread.daemon = True # Set as daemon, so when all other threads die, this one gets killed too. -read_thread.start() - -## Start Qt event loop unless running in interactive mode or using pyside. -if __name__ == '__main__': - import sys - if (sys.flags.interactive != 1) or not hasattr(QtCore, 'PYQT_VERSION'): - try: - QtGui.QApplication.instance().exec_() - except KeyboardInterrupt: - sys.exit(0) diff --git a/octave/gmsk.m b/octave/gmsk.m deleted file mode 100644 index a8c04a1..0000000 --- a/octave/gmsk.m +++ /dev/null @@ -1,1021 +0,0 @@ -% gmsk.m -% David Rowe Dec 2014 -% -% GMSK modem implementation and simulations to test - -% -% [X] plot eye diagram -% [X] BER curves with reas match to theoretical -% [X] fine timing estimator -% [X] test with fine timing error by resampling -% [X] phase/freq estimator -% + need initial acquisition and tracking -% [X] test with different freq offsets -% [X] coarse timing estimator (sync up to known test frames) -% [X] test with different coarse timing offsets -% [ ] file read/write interface -% [ ] refactor into tx/rx functions -% [X] modify for 1200 (or any) bit/s operation -% + ie GMSK filter coeff generation -% + or just re-sampling? e.g. ratio of Fs to Rs? -% [ ] way to measure input SNR to demod -% + Maybe based on test tone/carrier from the other side? -% + think about process ... total signal plus noise power? Increase power until S+N doubles? -% [X] generate curves for baseline modem and with sync algorithms -% [X] used coarse sync code to remove need for knowing delays -% [X] demod level indep -% + scaled OK +/- 20dB same BER -% [X] effect of DC signals from SDRs -% + simulated effect of general interferer at -1500Hz, at an amplitude of 4 -% (12 dB above GMSK signal), it started to affect BER e.g. 0.007 to 0.009 -% Line appeared about 30dB above top of GMSK signal. -% [ ] effect of quantisation noise - -% 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 low pass filter noise - -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,... -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]; - -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.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]; - -rand('state',1); -randn('state',1); -graphics_toolkit ("gnuplot"); -fm; -close all; - -% -% Functions that implement the GMSK modem ------------------------------------------------------ -% - -function gmsk_states = gmsk_init(gmsk_states, Rs) - - % general - - verbose = gmsk_states.verbose; - gmsk_states.Fs = 48000; - gmsk_states.Rs = Rs; - M = gmsk_states.M = gmsk_states.Fs/gmsk_states.Rs; - global gmsk_mod_coeff; - global gmsk_demod_coeff; - gmsk_states.mod_coeff = (Rs/4800)*resample(gmsk_mod_coeff, 4800, Rs); - - if verbose > 1 - figure; - plot(gmsk_mod_coeff,'r;original 4800;') - hold on; - plot(gmsk_states.mod_coeff,'g;interpolated;') - hold off; - title('GMSK pulse shaping filter') - end - - % set up FM modulator - - fm_states.Fs = gmsk_states.Fs; - fm_states.fc = 0; - fm_max = fm_states.fm_max = Rs/2; - fd = fm_states.fd = Rs/4; - 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); - -endfunction - - -function [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits) - M = gmsk_states.M; - nsym = length(tx_bits); - nsam = nsym*M; - verbose = gmsk_states.verbose; - - % NRZ sequence of symbols - - tx_symbols = zeros(1,nsam); - for i=1:nsym - tx_symbols(1+(i-1)*M:i*M) = -1 + 2*tx_bits(i); - end - - tx_filt = filter(gmsk_states.mod_coeff, 1, tx_symbols); - - if verbose > 1 - figure; - clf - plot(tx_filt(1:M*10)) - title('tx signal after filtering, before FM mod') - end - - tx = analog_fm_mod(gmsk_states.fm_states, tx_filt); -endfunction - - -function [rx_bits rx_int rx_filt] = gmsk_demod(gmsk_states, rx) - M = gmsk_states.M; - Rs = gmsk_states.Rs; - Fs = gmsk_states.Fs; - nsam = length(rx); - nsym = floor(nsam/M); - global gmsk_demod_coeff; - wd = 2*pi*gmsk_states.fm_states.fd/gmsk_states.Fs; - timing_angle_log = zeros(1,length(rx)); - rx_int = zeros(1,length(rx)); - - if gmsk_states.coherent_demod - - % See IEEE Trans on Comms, Muroyta et al, 1981, "GSM Modulation - % for Digital Radio Telephony" Fig 8: - - % matched filter - - rx_filt = filter(gmsk_states.mod_coeff, 1, rx); - - % Property of MSK that re and im arms are sequences of 2T - % long symbols, can be demodulated like QPSK with matched filter - % and integrate and dump. - - % integrate energy in symbols 2T long in re and im arms - % note this could be combined with matched filter - - rx_int = conv(rx_filt,ones(1,2*M)); - - % phase and fine frequency tracking and correction ------------------------ - - if gmsk_states.phase_track - - % DCO design from "Introduction To Phase-Lock Loop System Modeling", Wen Li - % http://www.ece.ualberta.ca/~ee401/parts/data/PLLIntro.pdf - - eta = 0.707; - wn = 2*pi*10*(Rs/4800); % (Rs/4800) -> found reducing the BW benifical with falling Rs - Ts = 1/Fs; - g1 = 1 - exp(-2*eta*wn*Ts); - g2 = 1 + exp(-2*eta*wn*Ts) - 2*exp(-eta*wn*Ts)*cos(wn*Ts*sqrt(1-eta*eta)); - Gpd = 2/pi; - Gvco = 1; - G1 = g1/(Gpd*Gvco); G2 = g2/(Gpd*Gvco); - %printf("g1: %e g2: %e G1: %e G2: %e\n", g1, g2, G1, G2); - - filt_prev = dco = lower = ph_err_filt = ph_err = 0; - dco_log = filt_log = zeros(1,nsam); - - % w is the ref sine wave at the timing clock frequency - % tw is the length of the window used to estimate timing - - k = 1; - tw = 200*M; - xr_log = []; xi_log = []; - w_log = []; - timing_clock_phase = 0; - timing_angle = 0; - timing_angle_log = zeros(1,nsam); - - for i=1:nsam - - % update sample timing estimate every tw samples - - if mod(i,tw) == 0 - l = i - tw+1; - xr = abs(real(rx_int(l:l+tw-1))); - xi = abs(imag(rx_int(l:l+tw-1))); - w = exp(j*(l:l+tw-1)*2*pi*(Rs/2)/Fs); - X = xr * w'; - timing_clock_phase = timing_angle = angle(X); - k++; - xr_log = [xr_log xr]; - xi_log = [xi_log xi]; - w_log = [w_log w]; - else - timing_clock_phase += (2*pi)/(2*M); - end - timing_angle_log(i) = timing_angle; - - rx_int(i) *= exp(-j*dco); - ph_err = sign(real(rx_int(i))*imag(rx_int(i)))*cos(timing_clock_phase); - lower = ph_err*G2 + lower; - filt = ph_err*G1 + lower; - dco = dco + filt; - filt_log(i) = filt; - dco_log(i) = dco; - end - - figure; - clf - subplot(211); - plot(filt_log); - title('PLL filter') - subplot(212); - plot(dco_log/pi); - title('PLL DCO phase'); - %axis([1 nsam -0.5 0.5]) - end - - % sample integrator output at correct timing instant - - timing_adj = timing_angle_log*2*M/(2*pi); - timing_adj_uw = unwrap(timing_angle_log)*2*M/(2*pi); - % Toff = floor(2*M+timing_adj); - Toff = floor(timing_adj_uw+0.5); - k = 1; - re_syms = im_syms = zeros(1,nsym/2); - - for i=2*M:2*M:nsam - if (i-Toff(i)+M) < nsam - re_syms(k) = real(rx_int(i-Toff(i))); - im_syms(k) = imag(rx_int(i-Toff(i)+M)); - end - %re_syms(k) = real(rx_int(i-10)); - %im_syms(k) = imag(rx_int(i+M-10)); - k++; - end - - figure - subplot(211) - stem(re_syms) - subplot(211) - stem(im_syms) - - figure; - clf - subplot(211) - plot(timing_adj); - title('Timing est'); - subplot(212) - plot(Toff); - title('Timing est unwrap'); - - % XORs/adders on the RHS of Muroyta et al Fig 8 (a) and (b). We - % simulate digital logic bit stream at clock rate Rs, even though - % we sample integrators at rate Rs/2. I can't explain how and why - % this logic works/is required. I think it can be worked out from - % comparing to MSK/OQPSK demod designs. - - l = length(re_syms); - l2 = 2*l; - re_bits = zeros(1,l2); - im_bits = zeros(1,l2); - clk_bits = zeros(1,l2); - for i=1:l-1 - re_bits(2*(i-1)+1) = re_syms(i) > 0; - re_bits(2*(i-1)+2) = re_syms(i) > 0; - im_bits(2*(i-1)+2) = im_syms(i) > 0; - im_bits(2*(i-1)+3) = im_syms(i) > 0; - clk_bits(2*(i-1)+1) = 0; - clk_bits(2*(i-1)+2) = 1; - end - - rx_bits = bitxor(bitxor(re_bits,im_bits), clk_bits); - rx_bits = rx_bits(2:length(rx_bits)-1); - else - % non-coherent demod - - % filter to get rid of most of noise before FM demod, but doesn't - % introduce any ISI - - fc = Rs/(Fs/2); - bin = firls(200,[0 fc*(1-0.05) fc*(1+0.05) 1],[1 1 0.01 0.01]); - rx_filt = filter(bin, 1, rx); - - % FM demod - - rx_diff = [ 1 rx_filt(2:nsam) .* conj(rx_filt(1:nsam-1))]; - rx_filt = (1/wd)*atan2(imag(rx_diff),real(rx_diff)); - - % low pass filter, trade off between ISI and removing noise - - rx_filt = filter(gmsk_demod_coeff, 1, rx_filt); - Toff = 7; - rx_bits = real(rx_filt(1+Toff:M:length(rx_filt)) > 0); - - end - -endfunction - - -% Initial frequency offset estimation. Look for line a centre -% frequency, which is the strongest component when ...101010... is -% used to modulate the GMSK signal. Note just searching for a single -% line will get false lock on random sine waves but that's OK for a -% PoC. It could be improved by checking for other lines, or -% demodulating the preamble and checking for bit errors. - -function [freq_offset_est ratio] = gmsk_est_freq_offset(gmsk_states, rx, verbose) - Fs = gmsk_states.Fs; - Rs = gmsk_states.Rs; - - % Suggest Rs/10 symbols of preamble (100ms), this works OK at - % Rs=4800 and Es/No = 6dB. The large, Fs sample FFT size is used - % for convenience (the bin resolution is 1 Hz), for real time we - % would decimate and use smaller FFT to save CPU and memory. - - ndft = Fs; - f = fft(rx .* hanning(length(rx))', ndft); - f = fftshift(f); - - start_bin = 1 + Fs/2-Rs/4; - stop_bin = start_bin + Rs/2; - [max_val max_bin] = max(abs(f(start_bin:stop_bin))); - - max_bin -= Rs/4 + 1; - if verbose > 1 - printf("ndft: %d start_bin: %d stop_bin: %d max_bin: %d\n", ndft, start_bin, stop_bin, max_bin); - end - - % calc ratio of line energy to total energy. For a valid preamble - % this was measured as about 0.20 to 0.25 depending on noise. - - sum_sq = sum(abs(f(start_bin:stop_bin)) .^ 2); - ratio = sqrt(max_val*max_val/sum_sq); - - % map max_bin to frequency offset - - freq_offset_est = max_bin; - - if verbose > 1 - printf("freq_offset_est: %f pk/rms ratio: %f \n", freq_offset_est, ratio); - figure; - clf - subplot(211) - plot(rx,'+') - title('rx signal on complex plane') - subplot(212) - plot(-Rs/4:Rs/4, 20*log10(abs(f(start_bin:stop_bin)))); - axis([-Rs/4 Rs/4 0 80]); - title('spectrum of rx signal'); - end - -endfunction - -% -% Functions for Testing the GMSK modem -------------------------------------------------------- -% - -function sim_out = gmsk_test(sim_in) - nsym = sim_in.nsym; - EbNodB = sim_in.EbNodB; - verbose = sim_in.verbose; - Rs = 4800; - - gmsk_states.verbose = verbose; - gmsk_states.coherent_demod = sim_in.coherent_demod; - gmsk_states.phase_track = 0; - gmsk_states = gmsk_init(gmsk_states, Rs); - M = gmsk_states.M; - Fs = gmsk_states.Fs; - Rs = gmsk_states.Rs; - Bfm = gmsk_states.fm_states.Bfm; - - for ne = 1:length(EbNodB) - aEbNodB = EbNodB(ne); - EbNo = 10^(aEbNodB/10); - variance = Fs/(Rs*EbNo); - - tx_bits = round(rand(1, nsym)); - %tx_bits = ones(1, nsym); - %tx_bits = zeros(1, nsym); - %tx_bits(1:2:nsym) = 0; - [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*exp(j*pi/2) + noise; - - [rx_bits rx_out rx_filt] = gmsk_demod(gmsk_states, rx(1:length(rx))); - - % search for frame location over a range - - Nerrs_min = nsym; Nbits_min = nsym; l = length(rx_bits); - for i=1:100; - Nerrs = sum(xor(rx_bits(i:l), tx_bits(1:l-i+1))); - if Nerrs < Nerrs_min - Nerrs_min = Nerrs; - Nbits_min = l; - end - end - - TERvec(ne) = Nerrs_min; - BERvec(ne) = Nerrs_min/Nbits_min; - - if verbose > 0 - printf("EbNo dB: %3.1f Nerrs: %d BER: %f BER Theory: %f\n", aEbNodB, Nerrs_min, BERvec(ne), 0.5*erfc(sqrt(0.75*EbNo))); - end - - if verbose > 1 - - if gmsk_states.coherent_demod == 0 - Toff = 0; dsam = M*30; - figure; - clf - eyesyms = 2; - plot(rx_filt(dsam+1+Toff:dsam+eyesyms*M+Toff)) - hold on; - for i=1:10 - st = dsam+1+Toff+i*eyesyms*M; - en = st + eyesyms*M; - plot(rx_filt(st:en)) - end - hold off; - %axis([dsam dsam+eyesyms*M -2 2]); - title('Eye Diagram'); - else - figure; - nplot = 16; - clf; - subplot(211) - plot(real(rx_filt(1:nplot*M))) - axis([1 nplot*M -1 1]) - title('Matched Filter'); - subplot(212) - plot(imag(rx_filt(1:nplot*M))) - axis([1 nplot*M -1 1]) - - figure; - nplot = 16; - clf; - subplot(211) - plot(real(rx_out(1:nplot*M))/(2*M)) - title('Integrator'); - axis([1 nplot*M -1 1]) - subplot(212) - plot(imag(rx_out(1:nplot*M)/(2*M))) - axis([1 nplot*M -1 1]) - end - - figure; - clf - subplot(211) - stem(tx_bits(1:20)) - title('Tx Bits') - subplot(212) - stem(rx_bits(1:20)) - title('Rx Bits') - - figure; - clf - subplot(211); - f = fft(rx); - Tx = 20*log10(abs(f)); - plot(Tx) - grid; - title('GMSK Demodulator Input Spectrum'); - axis([1 5000 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; - axis([1 5000 0 max(cs)]) - - printf("Bfm: %4.0fHz %3.0f%% power bandwidth %4.0fHz = %3.2f*Rb\n", Bfm, x*100, bw, bw/Rs); - - end - end - - sim_out.TERvec = TERvec; - sim_out.BERvec = BERvec; - sim_out.Rs = gmsk_states.Rs; -endfunction - - -function run_gmsk_single - sim_in.coherent_demod = 0; - sim_in.nsym = 4800; - sim_in.EbNodB = 10; - sim_in.verbose = 2; - - sim_out = gmsk_test(sim_in); -endfunction - - -% Generate a bunch of BER versus Eb/No curves for various demods - -function run_gmsk_curves - sim_in.coherent_demod = 1; - sim_in.nsym = 48000; - sim_in.EbNodB = 2:10; - sim_in.verbose = 1; - - gmsk_coh = gmsk_test(sim_in); - - sim_in.coherent_demod = 0; - gmsk_noncoh = gmsk_test(sim_in); - - Rs = gmsk_coh.Rs; - EbNo = 10 .^ (sim_in.EbNodB/10); - alpha = 0.75; % guess for BT=0.5 GMSK - gmsk_theory.BERvec = 0.5*erfc(sqrt(alpha*EbNo)); - - % BER v Eb/No curves - - figure; - clf; - semilogy(sim_in.EbNodB, gmsk_theory.BERvec,'r;GMSK theory;') - hold on; - semilogy(sim_in.EbNodB, gmsk_coh.BERvec,'g;GMSK sim coherent;') - semilogy(sim_in.EbNodB, gmsk_noncoh.BERvec,'b;GMSK sim non-coherent;') - 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; - clf; - semilogy(sim_in.EbNodB+RsOnB_dB, gmsk_theory.BERvec,'r;GMSK theory;') - hold on; - semilogy(sim_in.EbNodB+RsOnB_dB, gmsk_coh.BERvec,'g;GMSK sim coherent;') - semilogy(sim_in.EbNodB+RsOnB_dB, gmsk_noncoh.BERvec,'b;GMSK sim non-coherent;') - 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 - - -function [preamble_location freq_offset_est] = find_preamble(gmsk_states, M, npreamble, rx) - verbose = gmsk_states.verbose; - - % look through rx buffer and determine if there is a valid preamble. Use steps of half the - % preamble size in samples to try to bracket the pre-amble. - - preamble_step = npreamble*M/2; - ratio = 0; freq_offset_est = 0; preamble_location = 0; - ratio_log = []; - for i=1:preamble_step:length(rx)-preamble_step - [afreq_offset_est aratio] = gmsk_est_freq_offset(gmsk_states, rx(i:i+preamble_step-1), verbose); - ratio_log = [ratio_log aratio]; - if aratio > ratio - preamble_location = i; - ratio = aratio; - freq_offset_est = afreq_offset_est; - end - end - if verbose - printf("preamble location: %2.1f seconds est f_off: %5.1f Hz ratio: %3.2f\n", - preamble_location/gmsk_states.Fs, freq_offset_est, ratio); - figure; - plot(ratio_log); - title('Preamble ratio'); - end -endfunction - - -% attempt to perform "coarse sync" sync with the received frames, we -% check each frame for the best coarse sync position. Brute force -% approach, that would be changed for a real demod which has some -% sort of unique word. Start looking for valid frames 1 frame -% after start of pre-amble to give PLL time to lock - -function [total_errors total_bits Nerrs_log Nerrs_all_log errors_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits) - - Nerrs_log = zeros(1, nframes_rx); - Nerrs_all_log = zeros(1, nframes_rx); - total_errors = 0; - total_bits = 0; - framesize = length(tx_frame); - errors_log = []; - - for f=2:nframes_rx-1 - Nerrs_min = framesize; - for i=1:framesize; - st = (f-1)*framesize+i; en = st+framesize-1; - errors = xor(rx_bits(st:en), tx_frame); - Nerrs = sum(errors); - if Nerrs < Nerrs_min - Nerrs_min = Nerrs; - errors_min = errors; - end - end - Nerrs_all_log(f) = Nerrs_min; - if Nerrs_min/framesize < 0.1 - errors_log = [errors_log errors_min]; - Nerrs_log(f) = Nerrs_min; - total_errors += Nerrs_min; - total_bits += framesize; - end - end -endfunction - -function plot_spectrum(gmsk_states, rx, preamble_location, title_str) - Fs = gmsk_states.Fs; - st = preamble_location + gmsk_states.npreamble*gmsk_states.M; - sig = rx(st:st+Fs*0.5); - h = hanning(length(sig))'; - Rx=20*log10(abs(fftshift(fft(sig .* h, Fs)))); - figure; - plot(-Fs/2:Fs/2-1,Rx); - grid("minor"); - xlabel('Hz'); - ylabel('dB'); - topy = ceil(max(Rx)/10)*10; - axis([-4000 4000 topy-50 topy+10]) - title(title_str); -endfunction - -% Give the demod a hard time: frequency, phase, time offsets, sample clock difference - -function run_test_channel_impairments - Rs = 1200; - verbose = 1; - aEbNodB = 6; - phase_offset = pi/2; - freq_offset = -104; - timing_offset = 100E3; - sample_clock_offset_ppm = -500; - interferer_freq = -1500; - interferer_amp = 0; - nsym = 4800*2; - npreamble = 480; - - gmsk_states.npreamble = npreamble; - gmsk_states.verbose = verbose; - gmsk_states.coherent_demod = 1; - gmsk_states.phase_track = 1; - gmsk_states = gmsk_init(gmsk_states, Rs); - Fs = gmsk_states.Fs; - Rs = gmsk_states.Rs; - M = gmsk_states.M; - - % A frame consists of nsym random data bits. Some experimentation - % has shown they must be random-ish data (not say 11001100...) for - % timing estimator to work. However initial freq offset estimation - % is a lot easier with a 01010 type sequence, so we construct a - % frame with a pre-amble followed by frames of random data. - - framesize = 480; - nframes = floor(nsym/framesize); - tx_frame = round(rand(1, framesize)); - tx_bits = zeros(1,npreamble); - tx_bits(1:2:npreamble) = 1; - for i=1:nframes - tx_bits = [tx_bits tx_frame]; - end - - [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits); - - tx = resample(tx, 1E6, 1E6-sample_clock_offset_ppm); - tx = [zeros(1,timing_offset) tx]; - nsam = length(tx); - - if verbose > 1 - figure; - subplot(211) - st = timing_offset; en = st+M*10; - plot(real(tx(st:en))) - title('Real part of tx'); - subplot(212) - plot(imag(tx(st:en))) - title('Imag part of tx'); - end - - EbNo = 10^(aEbNodB/10); - variance = Fs/(Rs*EbNo); - noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); - w = (0:nsam-1)*2*pi*freq_offset/Fs + phase_offset; - interferer = interferer_amp*exp(j*interferer_freq*(2*pi/Fs)*(0:nsam-1)); - - rx = sqrt(2)*tx.*exp(j*w) + noise + interferer; - - % optional dump to file - - if 1 - fc = 1500; gain = 10000; - wc = 2*pi*fc/Fs; - w1 = exp(j*wc*(1:nsam)); - rx1 = gain*real(rx .* w1); - fout = fopen("rx_6dB.raw","wb"); - fwrite(fout, rx1, "short"); - fclose(fout); - end - - rx = rx1 .* conj(w1); - - [preamble_location freq_offset_est] = find_preamble(gmsk_states, M, npreamble, rx); - w_est = (0:nsam-1)*2*pi*freq_offset_est/Fs; - rx = rx.*exp(-j*w_est); - - plot_spectrum(gmsk_states, rx, preamble_location, "GMSK rx just after preamble"); - - % printf("ntx: %d nrx: %d ntx_bits: %d\n", length(tx), length(rx), length(tx_bits)); - - [rx_bits rx_out rx_filt] = gmsk_demod(gmsk_states, rx(preamble_location+framesize:nsam)); - nframes_rx = length(rx_bits)/framesize; - - % printf("ntx: %d nrx: %d ntx_bits: %d nrx_bits: %d\n", length(tx), length(rx), length(tx_bits), length(rx_bits)); - - [total_errors total_bits Nerrs_log Nerrs_all_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits); - - ber = total_errors/total_bits; - - printf("Eb/No: %3.1f f_off: %4.1f ph_off: %4.3f Nframes: %d Nbits: %d Nerrs: %d BER: %f\n", - aEbNodB, freq_offset, phase_offset, nframes_rx, total_bits, total_errors, ber); - - figure; - clf - subplot(211) - plot(Nerrs_log,'r;errors/frame counted for BER;'); - hold on; - plot(Nerrs_all_log,'g;all errors/frame;'); - hold off; - legend("boxoff"); - title('Bit Errors') - subplot(212) - stem(real(cumsum(Nerrs_log))) - title('Cumulative Bit Errors') - -endfunction - - -% Generates a Fs=48kHz raw file of 16 bit samples centred on 1500Hz, -% Suitable for transmitting with a SSB tx - -function gmsk_tx(tx_file_name) - rand('state',1); - Rs = 1200; - nsym = Rs*4; - framesize = 480; - npreamble = 480; - gain = 10000; - fc = 1500; - - gmsk_states.verbose = 0; - gmsk_states.coherent_demod = 1; - gmsk_states.phase_track = 1; - gmsk_states = gmsk_init(gmsk_states, Rs); - Fs = gmsk_states.Fs; - Rs = gmsk_states.Rs; - M = gmsk_states.M; - - % generate frame with preamble - - nframes = floor(nsym/framesize) - tx_frame = round(rand(1, framesize)); - tx_bits = zeros(1,npreamble); - tx_bits(1:2:npreamble) = 1; - for i=1:nframes - tx_bits = [tx_bits tx_frame]; - end - - [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits); - nsam = length(tx); - - wc = 2*pi*fc/Fs; - w = exp(j*wc*(1:nsam)); - tx = gain*real(tx .* w); - figure; - plot(tx(1:4000)) - fout = fopen(tx_file_name,"wb"); - fwrite(fout, tx, "short"); - fclose(fout); - -endfunction - - -% Reads a file of Fs=48kHz 16 bit samples centred on 1500Hz, and -% measures the BER. - -function gmsk_rx(rx_file_name, err_file_name) - rand('state',1); - - Rs = 1200; - framesize = 480; - npreamble = 480; - fc = 1500; - - gmsk_states.npreamble = npreamble; - gmsk_states.verbose = 1; - gmsk_states.coherent_demod = 1; - gmsk_states.phase_track = 1; - gmsk_states = gmsk_init(gmsk_states, Rs); - Fs = gmsk_states.Fs; - Rs = gmsk_states.Rs; - M = gmsk_states.M; - - tx_frame = round(rand(1, framesize)); - - % get real signal at fc offset and convert to baseband complex - % signal - - fin = fopen(rx_file_name,"rb"); - rx = fread(fin,"short")'; - fclose(fin); - rx = filter([1 -0.999],[1 -0.99],rx); - nsam = length(rx); - wc = 2*pi*fc/Fs; - w = exp(-j*wc*(1:nsam)); - rxbb = rx .* w; - - figure; - plot(rx); - - % find preamble - - [preamble_location freq_offset_est] = find_preamble(gmsk_states, M, npreamble, rxbb); - - % power of signal, averaged over window - % TODO: remove wave file header, scale of actual level - % filter so we measure only energy in our passband - % work out noise BW of filter. Use GMSK filter? - - [b a] = cheby2(6,40,[200 3000]/(Fs/2)); - %bpwr_lp = fir2([200,4000/(Fs/2)); - noise_bw = var(filter(b,a,randn(1,1E6))); - - rx_filt = filter(b, a, rx(1000:length(rx))); - npower_window = 200*M; - rx_power = conv(rx_filt.^2,ones(1,npower_window))/(npower_window); - rx_power_dB = 10*log10(rx_power); - figure; - subplot(211) - plot(rx_filt(1000:length(rx_filt))); - title('GMSK Power (narrow filter)'); - subplot(212) - plot(rx_power_dB); - axis([1 length(rx_power) max(rx_power_dB)-29 max(rx_power_dB)+1]) - grid("minor") - - % Work out where to sample N, and S+N - - noise_end = preamble_location - 2*npreamble*M; - noise_start = noise_end - Fs; - if noise_start < 1 - printf("Hmm, we really need >1 second of noise only before preamble to measure noise!\n"); - else - noise = mean(rx_power_dB(noise_start:noise_end)); - signal_noise_start = preamble_location + 2*npreamble*M; - signal_noise_end = signal_noise_start + Fs; - signal_noise = mean(rx_power_dB(signal_noise_start:signal_noise_end)); - hold on; - plot([noise_start noise_end],[noise noise],'color','r','linewidth',5); - plot([signal_noise_start signal_noise_end],[signal_noise signal_noise],'color','r','linewidth',5); - - % determine SNR - - noise_lin = 10 ^ (noise/10); - signal_noise_lin = 10 ^ (signal_noise/10); - signal_lin = signal_noise_lin - noise_lin; - signal = 10*log10(signal_lin); - snr = signal - noise; - fudge_factor = 3; % 3dB for single/double sided noise adjustment? Just a guess - CNo = snr + 10*log10(Fs*noise_bw) - fudge_factor; - EbNo = CNo - 10*log10(Rs); - - EbNo_lin = 10 .^ (EbNo/10); - alpha = 0.75; % guess for BT=0.5 GMSK - ber_theory = 0.5*erfc(sqrt(alpha*EbNo_lin)); - - printf("Estimated S: %3.1f N: %3.1f Nbw: %4.0f Hz SNR: %3.1f CNo: %3.1f EbNo: %3.1f BER theory: %f\n", - signal, noise, Fs*noise_bw, snr, CNo, EbNo, ber_theory); - - % FM signal is centred on 12 kHz and 16 kHz wide so lets also work out noise there - - [b a] = cheby2(6,40,[12000-8000 12000+8000]/(Fs/2)); - noise_bw_fm = var(filter(b,a,randn(1,1E6))); - - rx_filt_fm = filter(b, a, rx(1000:length(rx))); - rx_power_fm = conv(rx_filt_fm.^2,ones(1,npower_window))/(npower_window); - rx_power_dB_fm = 10*log10(rx_power_fm); - - noise = mean(rx_power_dB_fm(noise_start:noise_end))*ones(1, length(rx_power_fm)); - noise_lin = 10 .^ (noise/10); - - signal_lin = rx_power_fm - noise_lin; - signal = 10*log10(abs(signal_lin) + 1E-6); - snr = signal - noise; - - CNo = snr + 10*log10(Fs*noise_bw_fm) - fudge_factor; - - figure - plot(rx_power_dB_fm,'r;signal plus noise;'); - hold on; - plot(CNo,'g;C/No;'); - hold off; - top_fm = ceil(max(CNo)/10)*10; - axis([1 length(rx_power_dB_fm) 20 top_fm]) - grid("minor") - legend("boxoff"); - title('FM C/No'); - end - - % spectrum of a chunk of GMSK signal just after preamble - - plot_spectrum(gmsk_states, rx, preamble_location, "GMSK rx just after preamble"); - - % correct freq offset and demodulate - - w_est = (0:nsam-1)*2*pi*freq_offset_est/Fs; - rxbb = rxbb.*exp(-j*w_est); - st = preamble_location+npreamble*M; - %en = min(nsam,st + 4*framesize*M); - en = nsam; - gmsk_statres.verbose = 2; - [rx_bits rx_out rx_filt] = gmsk_demod(gmsk_states, rxbb(st:en)); - nframes_rx = length(rx_bits)/framesize; - - % count errors - - [total_errors total_bits Nerrs_log Nerrs_all_log errors_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits); - - ber = total_errors/total_bits; - - printf("Nframes: %d Nbits: %d Nerrs: %d BER: %f\n", - nframes_rx, total_bits, total_errors, ber); - - % Optionally save a file of bit errors so we can simulate the effect on Codec 2 - - if nargin == 2 - - % To simulate effects of these errors on Codec 2: - % $ ~/codec2-dev/octave$ ../build_linux/src/c2enc 1300 ../raw/hts1raw - | ../build_linux/src/insert_errors - - ssb7dbSNR.err 52 | ../build_linux/src/c2dec 1300 - - | play -t raw -r 8000 -s -2 - - % Note in this example I'm using the 1300 bit/s codec, it's sig more robust that 1200 bit/s, - % if we ran the GMSK modem at 1300 bit/s there would be a 10*log10(1300/1200) = 0.35dB SNR penalty - - fep=fopen(err_file_name,"wb"); fwrite(fep, errors_log, "short"); fclose(fep); - end - - figure; - clf - subplot(211) - plot(Nerrs_log,'r;errors/frame counted for BER;'); - hold on; - plot(Nerrs_all_log,'g;all errors/frame;'); - hold on; - title('Bit Errors') - legend("boxoff"); - subplot(212) - stem(real(cumsum(Nerrs_log))) - title('Cumulative Bit Errors') -endfunction - - -%run_gmsk_single -%run_gmsk_curves -%run_gmsk_init -%run_test_channel_impairments -gmsk_tx("test_gmsk.raw") -gmsk_rx("test_gmsk.raw") -%gmsk_rx("ssb25db.wav") -%gmsk_rx("~/Desktop/ssb_fm_gmsk_high.wav") -%gmsk_rx("~/Desktop/test_gmsk_28BER.raw") -%gmsk_rx("~/Desktop/gmsk_rec_reverse.wav") - diff --git a/octave/hp_filt.m b/octave/hp_filt.m deleted file mode 100644 index 1087bb9..0000000 --- a/octave/hp_filt.m +++ /dev/null @@ -1,12 +0,0 @@ -% hp_filt.m -% David Rowe 20 Feb 2012 - -function hp_filt(in_file, out_file) - fin = fopen(in_file,"rb"); - s = fread(fin,Inf,"short"); - b = fir1(256, 300/4000, "high"); - freqz(b); - s_hpf = filter(b,1,s); - fout = fopen(out_file,"wb"); - fwrite(fout, s_hpf, "short"); -endfunction |