diff options
| author | drowe67 <[email protected]> | 2023-07-07 14:22:31 +0930 |
|---|---|---|
| committer | David Rowe <[email protected]> | 2023-07-07 14:22:31 +0930 |
| commit | a291cfaf4ca3bb359852336f796aa5cc64568826 (patch) | |
| tree | 5d765b89334fb422e1b6b83a0ce5d1e75098c15e | |
| parent | 87025ac7d2a23192bf124ce8684626f821b099ea (diff) | |
some more Octave files rm-ed
41 files changed, 0 insertions, 5552 deletions
diff --git a/octave/hackrf_dc.m b/octave/hackrf_dc.m deleted file mode 100644 index a204e70..0000000 --- a/octave/hackrf_dc.m +++ /dev/null @@ -1,26 +0,0 @@ -% hackrf_dc.m -% -% David Rowe Nov 2015 -% -% Downconverts a HackRF IQ sample file to a lower sample rate -% -% To sample a -60dB signal: -% $ hackrf_transfer -r df1.iq -f 439200000 -n 10000000 -l 20 -g 40play file at 10.7MHz used: -% octave:25> d = hackrf_dc("df1.iq") - -function d = hackrf_dc(infilename) - Fs1 = 10E6; % input sample rate to HackRF - Fs2 = 96E3; % output sample rate - fc = 700E3; % offset to shift input by, HackRF doesn't like signals in the centre - - s1 = load_hackrf(infilename); - ls1 = length(s1); - ls1 = 20*Fs1; - t = 0:ls1-1; - - % shift down to baseband from Fc, not sure of rot90 rather than trasnpose operator ' - % to avoid unwanted complex conj - - s2 = rot90(s1(1:ls1)) .* exp(-j*2*pi*t*fc/Fs1); - d = resample(s2, Fs2, Fs1); -end diff --git a/octave/hackrf_twotone.m b/octave/hackrf_twotone.m deleted file mode 100644 index 1d5aa02..0000000 --- a/octave/hackrf_twotone.m +++ /dev/null @@ -1,24 +0,0 @@ -% hackrf_twotone.m -% -% David Rowe Nov 2015 -% -% Generates a two tone test signal that can be played out of HackRF -% -% To play file at 10.7MHz used: -% $ hackrf_transfer -t ../octave/twotone.iq -f 10000000 -a 0 -x 47 -% -% However 2nd harmonic at 21.4 was only -32dBC so not really useful for my application -% in testing an ADC - -Fs = 8E6; -fc = 2E6; -f1 = fc; -f2 = fc+1E3; -A = 127; -T = 2; - -N = T*Fs; -t = 0:N-1; -%s = A*exp(j*2*pi*t*f1/Fs) + A*exp(j*2*pi*t*f2/Fs); -s = A*exp(j*2*pi*t*f2/Fs); -save_hackrf("twotone.iq",s); diff --git a/octave/hackrf_uc.m b/octave/hackrf_uc.m deleted file mode 100644 index 0f55bd1..0000000 --- a/octave/hackrf_uc.m +++ /dev/null @@ -1,47 +0,0 @@ -% hackrf_uc.m -% -% David Rowe Nov 2015 -% -% Upconverts a real baseband sample file to a file suitable for input into a HackRF -% -% To play file at 10.7MHz used: -% octave:25> hackrf_uc("fsk_10M.iq","fsk_horus_rx_1200_96k.raw") -% $ hackrf_transfer -t ../octave/fsk_10M.iq -f 10000000 -a 1 -x 40 - -function hackrf_uc(outfilename, infilename) - pkg load signal; - Fs1 = 48E3; % input sample rate - Fs2 = 10E6; % output sample rate to HackRF - fc = 700E3-24E3; % offset to shift to, HackRF doesn't like signals in the centre - A = 100; % amplitude of signal after upc-nversion (max 127) - N = Fs1*20; - - fin = fopen(infilename,"rb"); - printf("1\n"); - s1 = fread(fin,"short"); - printf("1\n"); - fclose(fin); - printf("1\n"); - ls1 = length(s1); - printf("1\n"); - N = ls1; - % single sided freq shifts, we don't want DSB - printf("1\n"); - s1 = hilbert(s1(1:N)); - - % upsample to Fs2 - - M = Fs2/Fs1; - s2 = resample(s1(1:N),Fs2,Fs1); - ls2 = length(s2); - mx = max(abs(s2)); - t = 0:ls2-1; - printf("2\n"); - % shift up to Fc, note use of rot90 rather than trasnpose operator ' - % as we don't want complex conj, that would shift down in freq - - sout = rot90((A/mx)*s2) .* exp(j*2*pi*t*fc/Fs2); - - save_hackrf(outfilename,sout); - -end diff --git a/octave/horus_msg.txt b/octave/horus_msg.txt deleted file mode 100644 index 882a9b3..0000000 --- a/octave/horus_msg.txt +++ /dev/null @@ -1 +0,0 @@ - 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 diff --git a/octave/impulse_noise.m b/octave/impulse_noise.m deleted file mode 100644 index dee5c5d..0000000 --- a/octave/impulse_noise.m +++ /dev/null @@ -1,122 +0,0 @@ -% impulse_noise -% David Rowe May 2017 -% -% Experiments with impulsive noise and HF radio - -format; -more off; -rand('seed',1) - -% DFT function ------------------------------------------------ -% note k is on 0..K-1 format, unlike Octave fft() which is 1..K - -function H = calc_H(k, K, a, d) - L = length(d); - H = 0; - for i=1:L - H += a(i)*exp(-j*2*pi*k*d(i)/K); - end -endfunction - -% ----------------------------------------- -% PWM noise simulation -% ----------------------------------------- - -function pwm_noise - - Fs = 10E6; % sample rate of simulation - Fsig = 1E6; % frequency of our wanted signal - Fpwm = 255E3; % switcher PWM frequency - T = 1; % length of simulations in seconds - Nsam = T*Fs; - Nsamplot = 200; - Apwm = 0.1; - Asig = -40; % attenuation of wanted signal in dB - - % generate an impulse train with jitter to simulate switcher noise - - pwm = zeros(1,Fs); - Tpwm = floor(Fs/Fpwm); - pulse_positions_pwm = Tpwm*(1:T*Fpwm) + round(rand(1,T*Fpwm)); - - h_pwm = zeros(1,Nsam); - h_pwm(pulse_positions_pwm) = Apwm; - h_pwm = h_pwm(1:Nsam); - - % add in wanted signal and computer amplitude spectrum - - s = 10^(Asig/20)*cos(2*pi*Fsig*(1:Nsam)/Fs); - - h = h_pwm+s; - H = fft(h); - Hdb = 20*log10(abs(H)) - 20*log10(Nsam/2); - - figure(1); clf; - subplot(211) - plot(h(1:Nsamplot)); - subplot(212) - plot(Hdb(1:Nsam/2)); - axis([0 T*2E6 -120 0]); xlabel('Frequency Hz'); ylabel('Amplityude dBV'); grid; - - printf("pwm rms: %f signal rms: %f noise rms\n", std(h_pwm), std(s)); -endfunction - -% ----------------------------------------- -% Single pulse noise simulation -% ----------------------------------------- - -function pulse_noise - - % set up short pulse in wide window, consisting of two samples next - % to each other - - K = 1024; - a(1) = a(2) = 1; d(1) = 10; d(2) = d(1)+1; - h = zeros(1,K); - h(d(1)) = a(1); - h(d(2)) = a(2); - - % mag and phase spectrum, mag spectrum changes slowly - - figure(2); clf; - Hfft = fft(h); - subplot(311) - stem(h(1:100)); - axis([1 100 -0.2 1.2]); - subplot(312) - plot(abs(Hfft(1:K/2)),'+'); - title('Magnitude'); - subplot(313) - plot(angle(Hfft(1:K/2)),'+'); - title('Phase'); - - % simple test to estimate H(k+1) from H(k) -------------------- - - % brute force calculation - - k = 300; - H = zeros(1,K); - H(k-1) = calc_H(k-1, K, a, d); - H(k) = calc_H(k, K, a, d); - H(k+1) = calc_H(k+1, K, a, d); - - % calculation of k+1 from k using approximation that {d(i)} are - % close together compared to M, i.e it's a narrow pulse (assumes we - % can estimate d using other means) - - Hk1_ = exp(-j*2*pi*d(1)/K)*H(k); - - % plot zoomed in version around k to compare - - figure(3); clf; - plot(H(k-1:k+1),'b+','markersize', 10, 'linewidth', 2); - hold on; plot(Hk1_,'g+','markersize', 10, 'linewidth', 2); hold off; - title('H(k-1) .... H(k+1)'); - printf("H(k+1) match: %f dB\n", 20*log10(abs(H(k+1) - Hk1_))); -endfunction - -% Run various simulations here --------------------------------------------- - -%pwm_noise -pulse_noise - diff --git a/octave/load_comp.m b/octave/load_comp.m deleted file mode 100644 index b9fa686..0000000 --- a/octave/load_comp.m +++ /dev/null @@ -1,9 +0,0 @@ -% load_comp.m -% David Rowe Sep 2015 - -function s = load_comp(fn) - fs=fopen(fn,"rb"); - s = fread(fs,Inf,"float32"); - ls = length(s); - s = s(1:2:ls) + j*s(2:2:ls); -endfunction diff --git a/octave/load_hackrf.m b/octave/load_hackrf.m deleted file mode 100644 index 85bfb73..0000000 --- a/octave/load_hackrf.m +++ /dev/null @@ -1,11 +0,0 @@ -% load_hackrf.m -% -% David Rowe Oct 2015 - -function s = load_hackrf(fn) - fs = fopen(fn,"rb"); - iq = fread(fs,Inf,"schar"); - fclose(fs); - l = length(iq); - s = iq(1:2:l) + j*iq(2:2:l); -endfunction diff --git a/octave/load_rtlsdr.m b/octave/load_rtlsdr.m deleted file mode 100644 index 946cfb2..0000000 --- a/octave/load_rtlsdr.m +++ /dev/null @@ -1,11 +0,0 @@ -% load_rtlsdr.m -% -% David Rowe Oct 2015 - -function s = load_rtlsdr(fn) - fs = fopen(fn,"rb"); - iq = fread(fs,Inf,"uchar"); - fclose(fs); - l = length(iq); - s = iq(1:2:l) + j*iq(2:2:l); -endfunction diff --git a/octave/mancyfsk.m b/octave/mancyfsk.m deleted file mode 100644 index 1ee24ef..0000000 --- a/octave/mancyfsk.m +++ /dev/null @@ -1,500 +0,0 @@ -% mancyfsk.m -% David Rowe October 2015 -% -% Manchester encoded 2FSK & 4FSK simulation. -% -% Attempt to design a FSK waveform that can pass through legacy FM -% radios but still be optimally demodulated by SDRs. It doesn't have -% to be optimally demodulated by legacy radios. Trick is getting it -% to pass through 300-3000Hz audio filters in legacy radios. -% -% [X] code up modulator -% [X] manchester two bit symbols -% [X] plot spectrum -% [X] demodulate using analog FM and ideal demods -% [X] measure BER compared to ideal coherent FSK - -1; - -fm; % analog FM library - - -function states = legacyfsk_init(M,Rs) - Fs = states.Fs = 96000; - states.Rs = Rs; % symbol rate over channel - Ts = states.Ts = Fs/Rs; % symbol period in samples - states.M = M; % mFSK, either 2 or 4 - bpsym = state.Rb = log2(M); % bits per symbol over channel - rate = states.rate = 0.5; % Manchester code rate - nbits = 100; - nbits = states.nbits = 100; % number of payload data symbols/frame - nbits2 = states.nbits2 = nbits/rate; % number of symbols/frame over channel after manchester encoding - nsym = states.nsym = nbits2/log2(M); % number of symbols per frame - nsam = states.nsam = nsym*Ts; - - %printf(" Rs: %d M: %d bpsym: %d nbits: %d nbits2: %d nsym: %d nsam: %d\n", Rs, M, bpsym, nbits, nbits2, nsym, nsam); - - states.fc = states.Fs/4; - if states.M == 2 - states.f(1) = states.fc - Rs/2; - states.f(2) = states.fc + Rs/2; - else - states.f(1) = states.fc - 3*Rs/2; - states.f(2) = states.fc - Rs/2; - states.f(3) = states.fc + Rs/2; - states.f(4) = states.fc + 3*Rs/2; - end -endfunction - - -% test modulator function - -function tx = legacyfsk_mod(states, tx_bits) - Fs = states.Fs; - Ts = states.Ts; - Rs = states.Rs; - f = states.f; - M = states.M; - nsym = states.nsym; - - tx = zeros(Ts*length(tx_bits)/log2(M),1); - tx_phase = 0; - - step = log2(M); - k = 1; - for i=1:step:length(tx_bits) - if M == 2 - tone = tx_bits(i) + 1; - else - tone = (tx_bits(i:i+1) * [2 1]') + 1; - end - tx_phase_vec = tx_phase + (1:Ts)*2*pi*f(tone)/Fs; - tx((k-1)*Ts+1:k*Ts) = 2.0*cos(tx_phase_vec); k++; - tx_phase = tx_phase_vec(Ts) - floor(tx_phase_vec(Ts)/(2*pi))*2*pi; - end - -endfunction - - -function run_sim(sim_in) - - frames = sim_in.frames; - test_frame_mode = sim_in.test_frame_mode; - M = sim_in.M; - Rs = sim_in.Rs; - demod = sim_in.demod; - EbNodB = sim_in.EbNodB; - timing_offset = sim_in.timing_offset; - - % rx timing has been adjusted experimentally - - if Rs == 4800 - if demod == 1 - rx_timing = 4; - else - rx_timing = 0; - end - end - if Rs == 2400 - if demod == 1 - rx_timing = 40; - else - rx_timing = 0; - end - end - - % init fsk modem - - more off - rand('state',1); - randn('state',1); - states = legacyfsk_init(M,Rs); - Fs = states.Fs; - nbits = states.nbits; - nbits2 = states.nbits2; - Ts = states.Ts; - nsam = states.nsam; - rate = states.rate; - - % init analog FM modem - - fm_states.Fs = Fs; - fm_max = fm_states.fm_max = 3E3; - fd = fm_states.fd = 5E3; - fm_states.fc = states.fc; - - 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); - [b, a] = cheby1(4, 1, 300/Fs, 'high'); % 300Hz HPF to simulate FM radios - - % init sim states - - rx_bits_buf = zeros(1,2*nbits2); - Terrs = Tbits = 0; - state = 0; - nerr_log = []; - - % set up the channel noise. We have log(M)*rate payload bits/symbol - % we have log2(M) bits/symbol, and rate bits per payload symbol - % TODO: explain this better as Im confused! - - EbNo = 10^(EbNodB/10); - EsNo = EbNo*rate*log2(M); - variance = states.Fs/((states.Rs)*EsNo); - %printf("EbNodB: %3.1f EbNo: %3.2f EsNo: %3.2f\n", EbNodB, EbNo, EsNo); - - % set up the input bits - - if test_frame_mode == 1 - % test frame of bits, which we repeat for convenience when BER testing - test_frame = round(rand(1, nbits)); - 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, nbits*(frames+1))); - end - if test_frame_mode == 3 - % ...10101... sequence - tx_bits = zeros(1, nbits*(frames+1)); - tx_bits(1:2:length(tx_bits)) = 1; - end - - % Manchester Encoding ----------------------------------------------------------- - - % Manchester encoding, which removes DC term in baseband signal, - % making the waveform friendly to old-school legacy FM radios with - % voiceband filtering. The "code rate" is 0.5, which means we have - % encode one input bit into 2 output bits. The 2FSK encoder takes - % one input bit, the 4FSK encoder two input bits. - - tx_bits_encoded = zeros(1,length(tx_bits)*2); - fsk2_enc = [[1 0]; [0 1]]; - % -1.5 1.5 1.5 -1.5 -0.5 0.5 0.5 -0.5 - % 0 3 3 0 1 2 2 1 - fsk4_enc = [[0 0 1 1]; [1 1 0 0]; [0 1 1 0]; [1 0 0 1]]; - k=1; - if M == 2 - for i=1:2:length(tx_bits_encoded) - input_bit = tx_bits(k); k++; - tx_bits_encoded(i:i+1) = fsk2_enc(input_bit+1,:); - end - else - for i=1:4:length(tx_bits_encoded) - input_bits = tx_bits(k:k+1) * [2 1]'; k+=2; - tx_bits_encoded(i:i+3) = fsk4_enc(input_bits+1,:); - end - end - - % FSK Modulator -------------------------------------------------------------- - - % use ideal FSK modulator (note: need to try using analog FM modulator) - - tx = legacyfsk_mod(states, tx_bits_encoded); - noise = sqrt(variance)*randn(length(tx),1); - rx = tx + noise; - timing_offset_samples = round(timing_offset*Ts); - rx = [zeros(timing_offset_samples,1); rx]; - - % Demodulator ---------------------------------------------------------------------------- - - if demod == 1 - % use analog FM demodulator, aka a $40 Baofeng - - [rx_out rx_bb] = analog_fm_demod(fm_states, rx'); - if sim_in.hpf - rx_out_hp = filter(b,a,rx_out); - else - rx_out_hp = rx_out; - end - rx_filt = filter(ones(1,Ts),1,rx_out_hp); - rx_timing_sig = rx_filt; - - % TODO: for 4FSK determine amplitude/decn boundaries, choose closest to demod each symbol - - end - - if demod == 2 - - % optimal non-coherent demod at Rs - - rx_timing_sig = zeros(1,length(rx)); - for m=1:M - phi_vec = (1:length(rx))*2*pi*states.f(m)/Fs; - dc = rx' .* exp(-j*phi_vec); - rx_filt(m,:) = abs(filter(ones(1,Ts),1,dc)); - rx_timing_sig = rx_timing_sig + rx_filt(m,1:length(rx)); - end - end - - % Fine timing estimation ------------------------------------------------------ - - % Estimate fine timing using line at Rs/2 that Manchester encoding provides - % We need this to sync up to Manchester codewords. TODO plot signal and - % timing "line" we extract - - Np = length(rx_timing_sig); - w = 2*pi*(Rs)/Fs; - x = (rx_timing_sig .^ 2) * exp(-j*w*(0:Np-1))'; - norm_rx_timing = angle(x)/(2*pi) - 0.42; - %rx_timing = round(norm_rx_timing*Ts); - %printf("norm_rx_timing: %4.4f rx_timing: %d\n", norm_rx_timing, rx_timing); - - % Max likelihood decoding of Manchester encoded symbols. Search - % through all ML possibilities to extract bits. Use energy (filter - % output sq) - - % Manchester Decoding -------------------------------------------------------- - - if M == 2 - if demod == 1 - - % sample at optimum instant - - [tmp l] = size(rx_filt); - rx_filt_dec = rx_filt(:, Ts+rx_timing:Ts:l); - - [tmp l] = size(rx_filt_dec); - rx_bits = zeros(1,l); - k = 1; - for i=1:2:l-1 - ml = [rx_filt_dec(i)-rx_filt_dec(i+1) -rx_filt_dec(i)+rx_filt_dec(i+1)]; - [mx mx_ind] = max(ml); - rx_bits(k) = mx_ind-1; k++; - end - end - - if demod == 2 - - % sample at optimum instant - - [tmp l] = size(rx_filt); - rx_filt_dec = rx_filt(:, Ts+rx_timing:Ts:l); - - [tmp l] = size(rx_filt_dec); - rx_bits = zeros(1,l); - k = 1; - for i=1:2:l-1 - %ml = [rx_filt_dec(2,i)*rx_filt_dec(1,i+1) rx_filt_dec(1,i)*rx_filt_dec(2,i+1)]; - ml = [rx_filt_dec(2,i)+rx_filt_dec(1,i+1) rx_filt_dec(1,i)+rx_filt_dec(2,i+1)]; - [mx mx_ind] = max(ml); - rx_bits(k) = mx_ind-1; k++; - end - end - else % M == 4 - if demod == 1 - % TODO: 4FSK version of demod - rx_bits=tx_bits; - end - if demod == 2 - % sample at optimal instant - - [tmp l] = size(rx_filt); - rx_filt_dec = rx_filt(:, Ts+rx_timing:Ts:l); - [tmp l] = size(rx_filt_dec); - rx_bits = zeros(1,l); - - k = 1; - fsk4_dec = [[0 0]; [0 1]; [1 0]; [1 1]]; - for i=1:2:l-1 - %ml = [rx_filt_dec(1,i)*rx_filt_dec(4,i+1) rx_filt_dec(4,i)*rx_filt_dec(1,i+1) rx_filt_dec(2,i)*rx_filt_dec(3,i+1) rx_filt_dec(3,i)*rx_filt_dec(2,i+1)]; - ml = [(rx_filt_dec(1,i)+rx_filt_dec(4,i+1)) (rx_filt_dec(4,i)+rx_filt_dec(1,i+1)) (rx_filt_dec(2,i)+rx_filt_dec(3,i+1)) (rx_filt_dec(3,i)+rx_filt_dec(2,i+1))]; - [mx mx_ind] = max(ml); - rx_bits(k:k+1) = fsk4_dec(mx_ind,:); k+=2; - end - end - end - - % useful for getting decoding right - %tx_bits(1:20) - %rx_bits(1:20) - - % Frame sync and BER logic ------------------------------------------------------------- - - st = 1; - for f=1:frames - - % extract nin bits - - nin = nbits; - en = st + nin - 1; - - rx_bits_buf(1:nbits) = rx_bits_buf(nbits+1:2*nbits); - rx_bits_buf(nbits+1:2*nbits) = rx_bits(st:en); - - st += nin; - - % frame sync based on min BER - - if test_frame_mode == 1 - nerrs_min = nbits; - next_state = state; - if state == 0 - for i=1:nbits - error_positions = xor(rx_bits_buf(i:nbits+i-1), test_frame); - nerrs = sum(error_positions); - %printf("i: %d nerrs: %d nerrs_min: %d \n", i, nerrs, nerrs_min); - if nerrs < nerrs_min - nerrs_min = nerrs; - coarse_offset = i; - end - end - if nerrs_min < 3 - next_state = 1; - %printf("%d %d\n", coarse_offset, nerrs_min); - end - end - - if state == 1 - error_positions = xor(rx_bits_buf(coarse_offset:coarse_offset+nbits-1), test_frame); - nerrs = sum(error_positions); - Terrs += nerrs; - Tbits += nbits; - nerr_log = [nerr_log nerrs]; - end - - state = next_state; - - end - end - - if test_frame_mode == 1 - if sim_in.verbose - printf(" demod: %d frames: %d EbNodB: %3.1f Tbits: %d Terrs: %d BER %4.3f\n", demod, frames, EbNodB, Tbits, Terrs, Terrs/Tbits); - else - printf(" EbNodB: %3.1f BER %4.3f\n", EbNodB, Terrs/Tbits); - end - end - - % Bunch O'plots -------------------------------------------------------------- - - close all; - - st = 1; en=20; - - Tx=fft(tx, Fs); - TxdB = 20*log10(abs(Tx(1:Fs/2))); - figure(1) - clf; - plot(TxdB) - axis([1 Fs/2 (max(TxdB)-100) max(TxdB)]) - title('Tx Spectrum'); - - figure(2) - clf - if demod == 1 - subplot(211) - plot(rx_filt(st*Ts:en*Ts)); - title('After integrator'); - subplot(212) - plot(rx_filt_dec(st:en),'+'); - title('Decimated output'); - end - if demod == 2 - subplot(211); - plot(rx_filt(1,st*Ts:en*Ts)); - hold on; - plot(rx_filt(2,st*Ts:en*Ts),'g'); - if M == 4 - plot(rx_filt(3,st*Ts:en*Ts),'c'); - plot(rx_filt(4,st*Ts:en*Ts),'r'); - end - hold off; - title('Output of each filter'); - subplot(212); - plot(rx_filt_dec(1,st:en),'+'); - hold on; - plot(rx_filt_dec(2,st:en),'g+'); - if M == 4 - plot(rx_filt_dec(3,st:en),'c+'); - plot(rx_filt_dec(4,st:en),'r+'); - end - hold off; - title('Decimated output of each filter'); - end - - figure(3) - clf; - subplot(211) - plot(rx_timing_sig(st*Ts:en*Ts).^2) - title('rx-timing-sig') - subplot(212) - F = abs(fft(rx_timing_sig(1:Fs))); - plot(F(100:8000)) - title('FFT of rx-timing-sig') - - if demod == 1 - figure(4); - clf; - h = fft(rx_out, Fs); - hdB = 20*log10(abs(h)); - plot(hdB(1:4000)) - title('Spectrum of baseband modem signal after analog FM demod'); - axis([1 4000 (max(hdB)-40) max(hdB)]) - end - - if demod == 1 - figure(5) - clf; - subplot(211) - plot(rx_out(st*Ts:en*Ts)); - title('baseband modem signal after analog FM demod'); - subplot(212) - plot(rx_out_hp(st*Ts:en*Ts)); - title('baseband modem signal after 300Hz filter'); - end -end - - -% Run various permutations of simulation here --------------------------------------- - -function run_single - - sim_in.frames = 100; - sim_in.test_frame_mode = 1; - sim_in.M = 2; - sim_in.Rs = 2400; - sim_in.demod = 1; - sim_in.EbNodB = 15; - sim_in.timing_offset = 0.0; - sim_in.hpf = 1; - sim_in.verbose = 1; - - run_sim(sim_in); -endfunction - - -function run_lots - - % adjusted a few scenarios for about 2% BER so we can compare - - sim_in.frames = 100; - sim_in.test_frame_mode = 1; - sim_in.M = 2; - sim_in.Rs = 4800; - sim_in.demod = 1; - sim_in.EbNodB = 12; - sim_in.timing_offset = 0.0; - sim_in.hpf = 1; - sim_in.verbose = 0; - - printf("Rs=4800 2FSK ideal demod\n"); - sim_in.EbNodB = 8.5; sim_in.demod = 2; run_sim(sim_in); - printf("Rs=4800 2FSK analog FM demod, not too shabby and pushes 2400bit/s thru a $40 HT!\n"); - sim_in.EbNodB = 12; sim_in.demod = 1; run_sim(sim_in); - printf("Rs=2400 2FSK analog FM demod, needs more power for same BER! Che?\n"); - sim_in.Rs = 2400; sim_in.EbNodB = 15; run_sim(sim_in); - printf("Hmm, doesn't improve with no 300Hz HPF, maybe due to less deviation?\n"); - sim_in.hpf = 0; run_sim(sim_in); - printf("Rs=2400 4FSK ideal demod, nice low Eb/No!\n"); - sim_in.demod = 2; sim_in.M = 4; sim_in.Rs = 2400; sim_in.EbNodB = 6; run_sim(sim_in); -endfunction - -%run_single; -run_lots; diff --git a/octave/mfsk.m b/octave/mfsk.m deleted file mode 100644 index 0414b5e..0000000 --- a/octave/mfsk.m +++ /dev/null @@ -1,199 +0,0 @@ -% mfsk.m -% David Rowe Nov 2015 - -% Simulation to test m=2 and m=4 FSK demod - - -1; - -function sim_out = fsk_ber_test(sim_in) - Fs = 96000; - M = sim_in.M; - Rs = sim_in.Rs; - Ts = Fs/Rs; - verbose = sim_in.verbose; - - nbits = sim_in.nbits; - nsym = sim_in.nbits*2/M; - nsam = nsym*Ts; - EsNodB = sim_in.EbNodB + 10*log10(M/2); - - % printf("M: %d nbits: %d nsym: %d\n", M, nbits, nsym); - - if M == 2 - f(1) = -Rs/2; - f(2) = Rs/2; - end - if M == 4 - f(1) = -3*Rs/2; - f(2) = -Rs/2; - f(3) = Rs/2; - f(4) = 3*Rs/2; - end - - % simulate over a range of Eb/No values - - for ne = 1:length(EsNodB) - Nerrs = Terrs = Tbits = 0; - - aEsNodB = EsNodB(ne); - EsNo = 10^(aEsNodB/10); - variance = Fs/(Rs*EsNo); - - % Modulator ------------------------------- - - tx_bits = round(rand(1, nbits)); - tx = zeros(1,nsam); - tx_phase = 0; - - for i=1:nsym - if M == 2 - tone = tx_bits(i) + 1; - else - tone = (tx_bits(2*(i-1)+1:2*i) * [2 1]') + 1; - end - - tx_phase_vec = tx_phase + (1:Ts)*2*pi*f(tone)/Fs; - tx((i-1)*Ts+1:i*Ts) = exp(j*tx_phase_vec); - tx_phase = tx_phase_vec(Ts) - floor(tx_phase_vec(Ts)/(2*pi))*2*pi; - 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 - - % Demodulator ----------------------------- - - % non-coherent FSK demod - - rx_bb = rx; - dc = zeros(M,nsam); - for m=1:M - dc(m,:) = rx_bb .* exp(-j*(0:nsam-1)*2*pi*f(m)/Fs); - end - - rx_bits = zeros(1, nsym); - for i=1:nsym - st = (i-1)*Ts+1; - en = st+Ts-1; - for m=1:M - int(m,i) = abs(sum(dc(m,st:en))); - end - if m == 2 - rx_bits(i) = int(1,i) < int(2,i); - else - [max_amp tone] = max([int(1,i) int(2,i) int(3,i) int(4,i)]); - if tone == 1 - rx_bits(2*(i-1)+1:2*i) = [0 0]; - end - if tone == 2 - rx_bits(2*(i-1)+1:2*i) = [0 1]; - end - if tone == 3 - rx_bits(2*(i-1)+1:2*i) = [1 0]; - end - if tone == 4 - rx_bits(2*(i-1)+1:2*i) = [1 1]; - end - end - end - - error_positions = xor(rx_bits, tx_bits); - 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: %4.3f \n", aEsNodB - 10*log10(M/2), Terrs, Terrs/Tbits); - end - end - - sim_out.TERvec = TERvec; - sim_out.BERvec = BERvec; -endfunction - - -function run_fsk_curves - sim_in.M = 2; - sim_in.Rs = 1200; - sim_in.nbits = 12000; - sim_in.EbNodB = 0:2:20; - sim_in.verbose = 1; - - EbNo = 10 .^ (sim_in.EbNodB/10); - fsk_theory.BERvec = 0.5*exp(-EbNo/2); % non-coherent BFSK demod - fsk2_sim = fsk_ber_test(sim_in); - - sim_in.M = 4; - fsk4_sim = fsk_ber_test(sim_in); - - % BER v Eb/No curves - - figure(1); - clf; - semilogy(sim_in.EbNodB, fsk_theory.BERvec,'r;2FSK theory;') - hold on; - semilogy(sim_in.EbNodB, fsk2_sim.BERvec,'g;2FSK sim;') - semilogy(sim_in.EbNodB, fsk4_sim.BERvec,'b;4FSK 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)") - -end - - -function run_fsk_single - sim_in.M = 4; - sim_in.Rs = 1200; - sim_in.nbits = 5000; - sim_in.EbNodB = 8; - sim_in.verbose = 1; - - fsk_sim = fsk_ber_test(sim_in); -endfunction - - -rand('state',1); -randn('state',1); -graphics_toolkit ("gnuplot"); - -run_fsk_curves -%run_fsk_single - diff --git a/octave/newamp1_fbf.m b/octave/newamp1_fbf.m deleted file mode 100644 index 69111c6..0000000 --- a/octave/newamp1_fbf.m +++ /dev/null @@ -1,144 +0,0 @@ -% newamp1_fbf.m -% -% Copyright David Rowe 2016 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Interactive Octave script to explore frame by frame operation of newamp1 -% spectral amplitude modelling. -% -% Usage: -% Make sure codec2-dev is compiled with the -DDUMP option - see README for -% instructions. -% ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --dump hts1a -% $ cd ~/codec2-dev/octave -% octave:14> newamp1_fbf("../build_linux/src/hts1a",50) - - -function newamp1_fbf(samname, f=73, varargin) - more off; - - newamp_700c; melvq; - load train_120_1.txt; load train_120_2.txt; - train_120_vq(:,:,1)= train_120_1; train_120_vq(:,:,2)= train_120_2; m=5; - - Fs = 8000; K = 20; - - vq = 0; eq_en = 0; pf = 0; - - % load up text files dumped from c2sim --------------------------------------- - - 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); - [frames tmp] = size(model); - - % pre-process - [rate_K_surface sample_freqs_kHz] = resample_const_rate_f_mel(model(1:frames,:), K); - - % we need to know eq states on each frame - eq = zeros(frames,K); an_eq = zeros(1,K); - for ff=1:frames - mean_f = mean(rate_K_surface(ff,:)); - rate_K_vec_no_mean = rate_K_surface(ff,:) - mean_f; - [tmp an_eq] = front_eq(rate_K_vec_no_mean, an_eq); - eq(ff,:) = an_eq; - end - - % Keyboard loop -------------------------------------------------------------- - - k = ' '; - do - fg = 1; - s = [ Sn(2*f-1,:) Sn(2*f,:) ]; - figure(fg++); clf; plot(s); axis([1 length(s) -20000 20000]); - - Wo = model(f,1); L = model(f,2); Am = model(f,3:(L+2)); AmdB = 20*log10(Am); - Am_freqs_kHz = (1:L)*Wo*4/pi; - - % plots ---------------------------------- - - figure(fg++); clf; - l = sprintf(";rate %d AmdB;g+-", L); - plot((1:L)*Wo*4000/pi, AmdB, l); - axis([1 4000 -20 80]); - hold on; - stem(sample_freqs_kHz*1000, rate_K_surface(f,:), ";rate K;b+-"); - - % default - rate_K_vec_ = rate_K_surface(f,:); - - mean_f = mean(rate_K_surface(f,:)); - rate_K_vec_no_mean = rate_K_surface(f,:) - mean_f; - if eq_en - rate_K_vec_no_mean -= eq(f,:); - end - rate_K_vec_no_mean_ = rate_K_vec_no_mean; - if vq - [res rate_K_vec_no_mean_ ind] = mbest(train_120_vq, rate_K_vec_no_mean, m); - if pf - rate_K_vec_no_mean_ = post_filter(rate_K_vec_no_mean_, sample_freqs_kHz, 1.5); - end - rate_K_vec_ = rate_K_vec_no_mean_ + mean_f; - end - - % back to rate L - model_(f,:) = resample_rate_L(model(f,:), rate_K_vec_, sample_freqs_kHz); - Am_ = model_(f,3:(L+2)); AmdB_ = 20*log10(Am_); - varL = var(AmdB - AmdB_); - - plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-"); - l = sprintf(";error var %3.2f dB;bk+-", varL); - plot((1:L)*Wo*4000/pi, (AmdB - AmdB_), l); - hold off; - - figure(3); clf; - plot(sample_freqs_kHz*1000, 40+ rate_K_vec_no_mean, ";rate K no mean;g+-"); - axis([1 4000 -20 80]); hold on; - plot(sample_freqs_kHz*1000, 40 + rate_K_vec_no_mean_, ";rate K no mean bar;r+-"); - varK = var(rate_K_vec_no_mean - rate_K_vec_no_mean_); - l = sprintf(";error var %3.2f dB;bk+-", varK); - plot(sample_freqs_kHz*1000, rate_K_vec_no_mean - rate_K_vec_no_mean_, l); - - plot(sample_freqs_kHz*1000, eq(f,:), ";eq;b+-"); - hold off; - - % interactive menu ------------------------------------------ - - printf("\rframe: %d menu: n-next b-back q-quit v-vq[%d] p-pf[%d] e-eq[%d]", f, vq, pf, eq_en); - fflush(stdout); - k = kbhit(); - - if k == 'v' - if vq == 0; vq = 1; else vq = 0; end - endif - if k == 'p' - if pf == 0; pf = 1; else pf = 0; end - endif - if k == 'e' - if eq_en == 0; eq_en = 1; else eq_en = 0; end - endif - if k == 'n' - f = f + 1; - endif - if k == 'b' - f = f - 1; - endif - - until (k == 'q') - printf("\n"); - -endfunction - - -function ind = arg_exists(v, str) - ind = 0; - for i=1:length(v) - if !ind && strcmp(v{i}, str) - ind = i; - end - end -endfunction diff --git a/octave/nf_from_gr.m b/octave/nf_from_gr.m deleted file mode 100644 index bd36d02..0000000 --- a/octave/nf_from_gr.m +++ /dev/null @@ -1,129 +0,0 @@ -% nf_from_gr.m -% David Rowe Mar 2016 - -#{ - Calculate NF from GNU Radio output samples in - ...IQIQ... (32 bit float) sample files - - 1/ Take one sample with a -100dBm input carrier - 2/ Take another sample with no signal (just rx noise) - 3/ Set Fs, adjust st and en to use a chunk of spectrum without too - many birdies. - - Gotchas: - - 1/ Inspect Figure(1), the time domain plots. - 2/ Make sure plenty of ADC bits are being used with the noise-only sample, - we don't want ADC quantisation noise to dominate. Aim for about half - full scale. - 3/ Also watch out for clipping on either sample. - -#} - -1; - -function det_nf(p_filename, n_filename, title, Fs, st, en, Pin_dB, real_file=0) - - if real_file - % real samples files of 16 bit shorts - fs=fopen(p_filename,"rb"); - p = fread(fs,Inf,"short"); - fclose(fs); - fs=fopen(n_filename,"rb"); - pn = fread(fs,Inf,"short"); - fclose(fs); - else - % GNU radio complex file input - p = load_comp(p_filename); - pn = load_comp(n_filename); - end - - % skip any start up transients - - tst = floor(0.1*Fs); ten = st + Fs - 1; - P = fft(p(tst:ten)); - N = fft(pn(tst:ten)); - - PdB = 20*log10(abs(P)); - NdB = 20*log10(abs(N)); - - figure(1); clf; - subplot(211); plot(real(p(tst:tst+floor(Fs*0.1)))); - subplot(212); plot(real(pn(tst:tst+floor(Fs*0.1)))); - - figure(2); clf; - subplot(211); plot(st:en, PdB(st:en)); - subplot(212); plot(st:en, NdB(st:en)); - - #{ - ------------------------------------------------------------------------ - - From Wikipedia: The Noise Figure is the difference in decibels - (dB) between the noise output of the actual receiver to the noise - output of an “ideal” receiver - - An ideal receiver would have an output noise power of: - - Nout_dB = 10log10(B) -174 + G_dB - - The -174 dBm/Hz figure is the thermal noise density at 25C, for - every 1Hz of bandwidth your will get -174dBm of noise power. It's - the lower limit set by the laws of physics. G_dB is the Rx gain. The - 10log10(B) term takes into account the bandwidth of the Rx. A wider - bandwidth means more total noise power. - - So if you have a 1Hz bandwidth, and a gain of 100dB, you would - expect Nout_NdB = 0 -174 + 100 = -74dBm at the rx output with no - signal. If you have a 1000Hz bandwidth receiver you would have NdB_out - = 20 -174 + 100 = -44dBm of noise power at the output. - - To determine Noise Figure: - 1) Sample the Rx output first with a test signal and then with noise only. - 2) Find the Rx gain using the test signal. - 3) Find the noise output power, then using the gain we can find the noise - input power. - 4) Normalise the noise input power to 1Hz noise bandwidth and - compare to the thermal noise floor. - - ---------------------------------------------------------------------------- - #} - - % variance is the power of a sampled signal - - Pout_dB = 10*log10(var(P(st:en))); % Rx output power with test signal - G_dB = Pout_dB - Pin_dB; % Gain of Rx - Nout_dB = 10*log10(var(N(st:en))); % Rx output power with noise - Nin_dB = Nout_dB - G_dB; % Rx input power with noise - No_dB = Nin_dB - 10*log10(en-st); % Rx input power with noise in 1Hz bandwidth - NF_dB = No_dB + 174; % compare to thermal noise to get NF - printf("%10s: Pin: %4.1f Pout: %4.1f G: %4.1f NF: %3.1f dB\n", title, Pin_dB, Pout_dB, G_dB, NF_dB); -endfunction - - -% HackRF -------------------------- - -%p_filename = "~/Desktop/blogs/nf/hackrf_100dbm_4MHz.bin"; -%n_filename = "~/Desktop/blogs/nf/hackrf_nosignal_4MHz.bin"; -p_filename = "~/codec2-dev/build_linux/unittest/hackrf_100dbm_4MHz.bin"; -n_filename = "~/codec2-dev/build_linux/unittest/hackrf_nosignal_4MHz.bin"; -det_nf(p_filename, n_filename, "HackRF", 4E6, 180E3, 600E3, -100); - -#{ -% RTL-SDR -------------------------- - -p_filename = "~/Desktop/nf/neg100dBm_2MHz.bin"; -n_filename = "~/Desktop/nf/nosignal_2MHz.bin"; -det_nf(p_filename, n_filename, "RTL-SDR", 2E6, 100E3, 300E3, -100); - -% AirSpy ------------------------- - -p_filename = "~/Desktop/nf/airspy_100dbm_2.5MSPS.bin"; -n_filename = "~/Desktop/nf/airspy_nosig_2.5MSPS.bin"; -det_nf(p_filename, n_filename, "AirSpy", 2.5E6, 100E3, 300E3, -100); - -% Fun Cube Dongle Pro Plus ------------------------- - -p_filename = "~/Desktop/nf/fcdpp_100dbm_192khz.bin"; -n_filename = "~/Desktop/nf/fcdpp_nosig_192khz.bin"; -det_nf(p_filename, n_filename, "FunCube PP", 192E3, 25E3, 125E3, -100); -#} diff --git a/octave/nf_from_stdio.m b/octave/nf_from_stdio.m deleted file mode 100644 index e1d38f2..0000000 --- a/octave/nf_from_stdio.m +++ /dev/null @@ -1,133 +0,0 @@ -% nf_from_gr.m -% David Rowe Mar 2018 - -#{ - - Calculate NF in real time from 16 bit real samples from stdin - - 1/ Using gqrx: - - gqrx setup: - Configure I/O devices: - To switch on LNA bias for HackRF, in Configure I/O devices menu set: - Device String: hackrf,bias=1 - To switch on LNA bias for airspy run for a few seconds this before starting gqrx: - $ airspy_rx -r /dev/null -f 435 -b 1 - I used a sample rate of 250000 for Airspy R2, 3000000 for Airspy Mini - Input options...: start with set all gain sliders set to maximum - FFT Setting.....: freq Zoom to max - Receiver Options: On spectrum display, drag filter width until it's about 12k - Filter Shape Normal - Mode USB - Tune until tone is between 2 and 4 k - Press UDP button - - Then in a Linux Term: - - $ nc -ul 7355 | octave --no-gui -qf nf_from_stdio.m 48000 - - 2/ Using command line tools. Compile airspy tools and csdr from source: - - a) Airspy: - - $ airspy_rx -a 6000000 -l 14 -m 15 -v 15 -r - -f 434.998 -b 1 | \ - csdr convert_s16_f | csdr fir_decimate_cc 50 | csdr convert_f_s16 | \ - octave --no-gui -qf ~/codec2-dev/octave/nf_from_stdio.m 120000 complex - - Note: we tuned a few kHz down to put the test tone in the 2000 to 4000 Hz range. - - b) HackRF: - - Term 1: - - $ ~/codec2-dev/octave$ nc -ul 7355 | octave --no-gui -qf nf_from_stdio.m 80000 complex - - Term 2: - - $ hackrf_transfer -r - -f 434995000 -s 4000000 -a 1 -p 1 -l 40 -g 32 | \ - csdr convert_s8_f | csdr fir_decimate_cc 50 | csdr convert_f_s16 | \ - nc localhost -u 7355 - - Note: HackRF needed a bit of tuning to get test tone in 2000 to 4000 Hz range. This - can be tricky with the command line method, easier with gqrx. - - c) rtlsdr (assuming sig gen set to 144.5MHz, -100dBm) - - Term 1: - - $ ./rtl_sdr -g 50 -s 2400000 -f 144.498E6 - | csdr convert_u8_f | csdr fir_decimate_cc 50 | \ - csdr convert_f_s16 | octave --no-gui -qf ~/codec2/octave/nf_from_stdio.m 48000 complex - - TODO: - [ ] work out why noise power st bounces around so much, signal power seems stable - [ ] reduce CPU load, in particular of plotting -#} - -graphics_toolkit ("gnuplot") - -% command line arguments - -arg_list = argv (); -if nargin == 0 - printf("\nusage: %s FsHz [real|complex] [testToneLeveldBm]\n\n", program_name()); - exit(0); -end - -Fs = str2num(arg_list{1}); -shorts_per_sample = 1; - -if nargin == 2 - if strcmp(arg_list{2}, "real") - shorts_per_sample = 1; - end - if strcmp(arg_list{2}, "complex") - shorts_per_sample = 2; - end -end - -Pin_dB = -100; % level of input test tone -if nargin == 3 - Pin_dB = str2num(arg_list{3}); -end - -printf("Fs: %d shorts_per_sample: %d Pin_dB: %f\n", Fs, shorts_per_sample, Pin_dB); - -[s,c] = fread(stdin, shorts_per_sample*Fs, "short"); - -while c - if shorts_per_sample == 2 - s = s(1:2:end)+j*s(2:2:end); - end - S = fft(s.*hanning(Fs)); - SdB = 20*log10(abs(S)); - figure(1); plot(real(s)); axis([0 Fs -4E4 4E4]); - figure(2); plot(SdB); axis([0 12000 40 180]); - - % assume sine wave is between 2000 and 4000 Hz, and dominates energy in that - % region. Noise is between 5000 - 10000 Hz - - sig_st = 2000; sig_en = 5000; - noise_st = 6000; noise_en = 10000; - - % find peak and sum power a few bins either side, this ensure we don't capture - % too much noise as well - - [pk pk_pos] = max(abs(S)); - if pk_pos > 5 - Pout_dB1 = 10*log10(sum(abs(S(pk_pos-5:pk_pos+5)).^2)); % Rx output power with test signal - else - Pout_dB1 = 0; - end - - Pout_dB = 10*log10(sum(abs(S(sig_st:sig_en)).^2)); % Rx output power with test signal - G_dB = Pout_dB - Pin_dB; % Gain of Rx - Nout_dB = 10*log10(sum(abs(S(noise_st:noise_en)).^2)/(noise_en-noise_st)); % Rx output power with noise - Nin_dB = Nout_dB - G_dB; % Rx input power with noise - No_dB = Nin_dB; %- 10*log10(noise_en-noise_st); % Rx input power with noise in 1Hz bandwidth - NF_dB = No_dB + 174; % compare to thermal noise to get NF - printf("Pout: %4.1f %d %4.1f Nout: %4.1f G: %4.1f No: %4.1f NF: %3.1f dB\n", Pout_dB, pk_pos, Pout_dB1, Nout_dB, G_dB, No_dB, NF_dB); - - pause(2); - [s,c] = fread(stdin, shorts_per_sample*Fs, "short"); -endwhile - diff --git a/octave/oqpsk.m b/octave/oqpsk.m deleted file mode 100644 index f308617..0000000 --- a/octave/oqpsk.m +++ /dev/null @@ -1,521 +0,0 @@ -% oqpsk.m -% David Rowe Jan 2017 -% -% Unfiltered OQPSK modem implementation and simulations to test, -% derived from GMSK modem in gmsk.m -% -% Usage: see "choose one of these to run" at the end of this file. - -rand('state',1); -randn('state',1); -graphics_toolkit ("gnuplot"); -format -more off; - -% init modem states - -function oqpsk_states = oqpsk_init(oqpsk_states, Rs) - - % general - - verbose = oqpsk_states.verbose; - oqpsk_states.Fs = 4*Rs; - oqpsk_states.Rs = Rs; - oqpsk_states.bps = 2; % two bit/symbol for QPSK - - M = oqpsk_states.M = oqpsk_states.Fs/oqpsk_states.Rs; - assert(floor(M) == M, "oversampling factor M must be an integer"); - assert(floor(M/2) == M/2, "(oversampling factor M)/2 must be an integer to offset QPSK"); -endfunction - - -% Gray coded QPSK modulation function - -function symbol = qpsk_mod(two_bits) - two_bits_decimal = sum(two_bits .* [2 1]); - switch(two_bits_decimal) - case (0) symbol = 1; - case (1) symbol = j; - case (2) symbol = -j; - case (3) symbol = -1; - endswitch -endfunction - - -% Gray coded QPSK demodulation function - -function two_bits = qpsk_demod(symbol) - if isscalar(symbol) == 0 - printf("only works with scalars\n"); - return; - end - bit0 = real(symbol*exp(j*pi/4)) < 0; - bit1 = imag(symbol*exp(j*pi/4)) < 0; - two_bits = [bit1 bit0]; -endfunction - - -% Unfiltered OQPSK modulator - -function [tx tx_symb] = oqpsk_mod(oqpsk_states, tx_bits) - M = oqpsk_states.M; - bps = oqpsk_states.bps; - nsym = length(tx_bits)/bps; - nsam = nsym*M; - verbose = oqpsk_states.verbose; - - % Map bits to Gray coded QPSK symbols - - tx_symb = zeros(1,nsym); - - for i=1:nsym - tx_symb(i) = qpsk_mod(tx_bits(2*i-1:2*i))*exp(j*pi/4); - end - - % Oversample by M (sample and hold) to create unfiltered QPSK - - tx = zeros(1, nsam); - for i=1:nsym - tx((i-1)*M+1:(i*M)) = tx_symb(i); - end - - % delay Q arm by half of a symbol to make OQPSK - - tx = [real(tx) zeros(1,M/2)] + j*[zeros(1,M/2) imag(tx)]; -endfunction - - -#{ - - Unfiltered OQPSK demodulator function, with (optional) phase and - timing estimation. Adapted from Fig 8 of [1]. See also gmsk.m and - [2]. - - Note demodulator returns phase corrected symbols sampled at ideal - timing instant. These symbols may have a m*pi/2 phase ambiguity due - to properties of phase tracking loop. The caller is responsible for - determining this ambiguity and recovering the actual bits. - - [1] GMSK demodulator in IEEE Trans on Comms, Muroyta et al, 1981, - "GSM Modulation for Digital Radio Telephony". - - [2] GMSK Modem Simulation, http://www.rowetel.com/?p=3824 - -#} - - -function [rx_symb rx_int filt_log dco_log timing_adj Toff] = oqpsk_demod(oqpsk_states, rx) - M = oqpsk_states.M; - Rs = oqpsk_states.Rs; - Fs = oqpsk_states.Fs; - nsam = length(rx); - nsym = floor(nsam/M); - verbose = oqpsk_states.verbose; - - timing_angle_log = zeros(1,length(rx)); - rx_int = zeros(1,length(rx)); - dco_log = filt_log = zeros(1,nsam); - - % Unfiltered PSK - integrate energy in symbols M long in re and im arms - - rx_int = conv(rx,ones(1,M))/M; - - % phase and fine frequency tracking and correction ------------------------ - - if oqpsk_states.phase_est - - % 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 beneficial 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; - end - - if oqpsk_states.timing_est - % w is the ref sine wave at the timing clock frequency - % tw is the length of the window used to estimate timing - - tw = 200*M; - k = 1; - xr_log = []; xi_log = []; - w_log = []; - timing_clock_phase = 0; - timing_angle = 0; - timing_angle_log = zeros(1,nsam); - end - - % Sample by sample processing loop for timing and phase est. Note - % this operates at sample rate Fs, unlike many PSK modems that - % operate at the symbol rate Rs - - for i=1:nsam - - if oqpsk_states.timing_est - - % update sample timing estimate every tw samples, free wheel - % rest of the time - - 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/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)/M; - end - timing_angle_log(i) = timing_angle; - end - - if oqpsk_states.phase_est - - % PLL per-sample processing - - 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_log(i) = dco; - dco = dco + filt; - filt_log(i) = filt; - - end - - end - - % final adjustment of timing output to take into account slowly - % moving estimates due to sample clock offset. Unwrap ensures that - % when timing angle jumps from -pi to pi we move to the next symbol - % and frame sync isn't broken - - timing_adj = timing_angle_log*M/(2*pi); - timing_adj_uw = unwrap(timing_angle_log)*M/(2*pi); - % Toff = floor(2*M+timing_adj); - Toff = floor(timing_adj_uw+0.5); - - % sample integrator output at correct timing instant - - k = 1; - re_syms = im_syms = zeros(1,nsym); - rx_symb = []; - for i=M:M:nsam - if i-Toff(i)+M/2 <= nsam - re_syms(k) = real(rx_int(i-Toff(i))); - im_syms(k) = imag(rx_int(i-Toff(i)+M/2)); - %re_syms(k) = real(rx_int(i)); - %im_syms(k) = imag(rx_int(i)); - rx_symb = [rx_symb re_syms(k) + j*im_syms(k)]; - k++; - end - end - -endfunction - - -% Test modem over a range Eb/No points in an AWGN channel. Can -% simulate a variety of channel impairments and performs ambiguity -% resolution. - -function sim_out = oqpsk_test(sim_in) - bitspertestframe = sim_in.bitspertestframe; - nbits = sim_in.nbits; - EbNodB = sim_in.EbNodB; - verbose = sim_in.verbose; - Rs = 4800; - - oqpsk_states.verbose = verbose; - oqpsk_states.coherent_demod = sim_in.coherent_demod; - oqpsk_states.phase_est = sim_in.phase_est; - oqpsk_states.timing_est = sim_in.timing_est; - oqpsk_states = oqpsk_init(oqpsk_states, Rs); - M = oqpsk_states.M; - Fs = oqpsk_states.Fs; - Rs = oqpsk_states.Rs; - sample_clock_offset_ppm = sim_in.sample_clock_offset_ppm; - - tx_testframe = round(rand(1, bitspertestframe)); - ntestframes = floor(nbits/bitspertestframe); - tx_bits = []; - for i=1:ntestframes - tx_bits = [tx_bits tx_testframe]; - end - - for ne = 1:length(EbNodB) - aEbNodB = EbNodB(ne); - EbNo = 10^(aEbNodB/10); - variance = Fs/(Rs*EbNo*oqpsk_states.bps); - - [tx tx_symb] = oqpsk_mod(oqpsk_states, tx_bits); - if sample_clock_offset_ppm - tx = resample(tx, 1E6, 1E6-sample_clock_offset_ppm); - end - nsam = length(tx); - - phi = sim_in.phase_offset + 2*pi*sim_in.freq_offset*(1:nsam)/M; - - noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); - st = 1+sim_in.timing_offset; en = length(tx); - rx = tx(st:en).*exp(j*phi(st:en)) + noise(st:en); - - [rx_symb rx_int filt_log dco_log timing_adj Toff] = oqpsk_demod(oqpsk_states, rx); - - % OK so the phase and timing estimators get us close (e.g. a good - % scatter diagram), but no banana just yet. One problem is the - % PLL can lock up on mulitples of pi/2. Combinations of phase - % offsets can confuse the timing estimator. One tricky example is a - % phase offset of pi/2 which swaps I & Q, and with OQPSK (unlike - % MSK and friends) we can't easily tell which is I and which is Q - % after a phase rotation, e.g. could be IQIQIQI or QIQIQIQ - - % So we need to determine the ambiguities: - % a) could be m*pi/2 rotations of phase - % b) could be I and Q swapped by timing est - % c) time alignment of test frame - - nsymb = bitspertestframe/oqpsk_states.bps; - nrx_symb = length(rx_symb); - rx_bits = zeros(1, bitspertestframe); - atx_symb = tx_symb(1:nsymb); - - % Treat I and Q as separate sequences, each with their own unique - % word. In our case the UW is the whole test frame. Correlate rx - % sequence with tx sequence at each possible offset through the - % received symbols to find the test frames. Note we also - % correlate I of tx with Q of rx to trap any IQ swaps. - - % The sign of the I and Q correlation lets us sort out the pi/2 - % phase rotation issue. - - nerrs_tot = 0; nbits_tot = 0; - - max_corr = real(atx_symb) * real(atx_symb)'; - for offset=2:nrx_symb-nsymb+1 - corr_ii(offset) = real(atx_symb) * real(rx_symb(offset:offset+nsymb-1))'/max_corr; - corr_qq(offset) = imag(atx_symb) * imag(rx_symb(offset:offset+nsymb-1))'/max_corr; - corr_iq(offset) = real(atx_symb) * imag(rx_symb(offset:offset+nsymb-1))'/max_corr; - corr_qi(offset) = imag(atx_symb) * real(rx_symb(offset:offset+nsymb-1))'/max_corr; - %printf("offset: %2d ii: % 5f qq: % 5f iq: % 5f qi: % 5f\n", - %offset, corr_ii(offset), corr_qq(offset), corr_iq(offset), corr_qi(offset)); - - if abs(corr_ii(offset)) > 0.8 - - % no IQ swap, or time offset - - i_sign = sign(corr_ii(offset)); - q_sign = sign(corr_qq(offset)); - arx_symb = i_sign*real(rx_symb(offset:offset+nsymb-1)) + j*q_sign*imag(rx_symb(offset:offset+nsymb-1)); - - for i=1:nsymb - rx_bits(2*i-1:2*i) = qpsk_demod(arx_symb(i)*exp(-j*pi/4)); - end - nerrs = sum(xor(tx_testframe, rx_bits)); - if verbose > 2 - printf("offset: %5d swap: %d i_sign: % 2.1f q_sign: % 2.1f nerr: %d\n", - offset, 0, i_sign, q_sign, nerrs); - end - nerrs_tot += nerrs; - nbits_tot += bitspertestframe; - end - - if abs(corr_qi(offset)) > 0.8 - - % IQ swap, I part in Q part of symbol before - - i_sign = sign(corr_iq(offset-1)); - q_sign = sign(corr_qi(offset)); - arx_symb = i_sign*imag(rx_symb(offset-1:offset+nsymb-2)) + j*q_sign*real(rx_symb(offset:offset+nsymb-1)); - - for i=1:nsymb - rx_bits(2*i-1:2*i) = qpsk_demod(arx_symb(i)*exp(-j*pi/4)); - end - nerrs = sum(xor(tx_testframe, rx_bits)); - if verbose > 1 - printf("offset: %5d swap: %d i_sign: % 2.1f q_sign: % 2.1f nerr: %d\n", - offset, 1, i_sign, q_sign, nerrs); - end - nerrs_tot += nerrs; - nbits_tot += bitspertestframe; - end - end - - TERvec(ne) = nerrs_tot; - BERvec(ne) = nerrs_tot/nbits_tot; - - if verbose > 0 - printf("EbNo dB: %3.1f Nbits: %d Nerrs: %d BER: %4.3f BER Theory: %4.3f\n", - aEbNodB, nbits_tot, nerrs_tot, BERvec(ne), 0.5*erfc(sqrt(EbNo))); - end - - if find(sim_in.plots == 1) - figure(1); clf; - subplot(211) - stem(real(tx)) - title('Tx samples'); - ylabel('Inphase'); - subplot(212) - stem(imag(tx)) - ylabel('Quadrature'); - end - - if find(sim_in.plots == 2) - figure(2); clf; - f = fftshift(fft(rx)); - Tx = 20*log10(abs(f)); - plot((1:length(f))*Fs/length(f) - Fs/2, Tx) - grid; - title('OQPSK Demodulator Input Spectrum'); - end - - if find(sim_in.plots == 3) - figure(3); clf; - nplot = min(16, nbits/oqpsk_states.bps); - title('Rx Integrator'); - subplot(211) - stem(real(rx_int(1:nplot*M))) - axis([1 nplot*M -1 1]) - subplot(212) - stem(imag(rx_int(1:nplot*M))) - axis([1 nplot*M -1 1]) - end - - if find(sim_in.plots == 4) - figure(4); clf; - subplot(211); - plot(filt_log); - title('PLL filter') - subplot(212); - plot(dco_log); - title('PLL DCO phase'); - end - - if find(sim_in.plots == 5) - figure(5); clf; - subplot(211) - plot(timing_adj); - title('Timing est'); - subplot(212) - plot(Toff); - title('Timing est unwrap'); - end - - if find(sim_in.plots == 6) - figure(6); clf; - st = floor(0.5*nrx_symb); - plot(rx_symb(st:nrx_symb), '+'); - title('Scatter Diagram'); - axis([-1.5 1.5 -1.5 1.5]) - end - - if find(sim_in.plots == 7) - figure(7); clf; - subplot(211) - plot(corr_ii); - axis([1 length(corr_ii) -1.2 1.2]); - title('corr ii'); - subplot(212) - plot(corr_qi); - axis([1 length(corr_ii) -1.2 1.2]); - title('corr qi'); - end - - if find(sim_in.plots == 8) - figure(8); clf; - subplot(211); - stem(real(arx_symb)); - title('Rx symbols') - subplot(212); - stem(imag(arx_symb)); - end - - if find(sim_in.plots == 9) - figure(9); clf; - subplot(211) - stem(tx_testframe(1:min(20,length(rx_bits)))) - title('Tx Bits') - subplot(212) - stem(rx_bits(1:min(20,length(rx_bits)))) - title('Rx Bits') - end - end - - sim_out.TERvec = TERvec; - sim_out.BERvec = BERvec; - sim_out.Rs = oqpsk_states.Rs; -endfunction - - -function run_oqpsk_single - sim_in.coherent_demod = 1; - sim_in.phase_est = 1; - sim_in.timing_est = 1; - sim_in.bitspertestframe = 100; - sim_in.nbits = 10000; - sim_in.EbNodB = 4; - sim_in.verbose = 1; - sim_in.phase_offset = 3*pi/4; % in radians - sim_in.timing_offset = 4; % in samples 0..M-1 - sim_in.freq_offset = 0.001; % fraction of Symbol Rate - sim_in.plots = [1 2 4 5 6 7]; - sim_in.sample_clock_offset_ppm = 100; - - sim_out = oqpsk_test(sim_in); -endfunction - - -% Generate a bunch of BER versus Eb/No curves for various demods - -function run_oqpsk_curves - sim_in.coherent_demod = 1; - sim_in.EbNodB = 2:8; - sim_in.verbose = 1; - sim_in.phase_est = 1; - sim_in.timing_est = 1; - sim_in.bitspertestframe = 100; - sim_in.nbits = 50000; - sim_in.phase_offset = 3*pi/4; % in radians - sim_in.timing_offset = 4; % in samples 0..M-1 - sim_in.freq_offset = 0.001; % fraction of Symbol Rate - sim_in.plots = []; - sim_in.sample_clock_offset_ppm = 0; - - oqpsk_coh = oqpsk_test(sim_in); - - Rs = oqpsk_coh.Rs; - EbNo = 10 .^ (sim_in.EbNodB/10); - oqpsk_theory.BERvec = 0.5*erfc(sqrt(EbNo)); - - % BER v Eb/No curves - - figure; - clf; - semilogy(sim_in.EbNodB, oqpsk_theory.BERvec,'r+-;OQPSK theory;') - hold on; - semilogy(sim_in.EbNodB, oqpsk_coh.BERvec,'g+-;OQPSK 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)") -endfunction - - -% Choose one of these to run ------------------------------------------ - -run_oqpsk_single -%run_oqpsk_curves - diff --git a/octave/papr_test.m b/octave/papr_test.m deleted file mode 100644 index 3785e25..0000000 --- a/octave/papr_test.m +++ /dev/null @@ -1,407 +0,0 @@ -% papr_test.m -% -% Experiments with PAPR reduction using clipping/compression -% -% OFDM Tx -> compress -> filter -> normalise power -> channel -> OFDM Rx - -#{ - TODO: - [ ] option for normalised power after clipper - [ ] experiment to plot those curves -#} - -1; - -function symbol = qpsk_mod(two_bits) - two_bits_decimal = sum(two_bits .* [2 1]); - switch(two_bits_decimal) - case (0) symbol = 1; - case (1) symbol = j; - case (2) symbol = -j; - case (3) symbol = -1; - endswitch -endfunction - -function two_bits = qpsk_demod(symbol) - bit0 = real(symbol*exp(j*pi/4)) < 0; - bit1 = imag(symbol*exp(j*pi/4)) < 0; - two_bits = [bit1 bit0]; -endfunction - -function papr = calc_papr(tx) - papr = 10*log10(max(abs(tx).^2)/mean(abs(tx).^2)); -end - -% test PAPR calculation with a two tone signal of known PAPR (3dB) -function test_papr - f1=800; f2=1200; Fs=8000; n=(0:Fs-1); - tx=exp(j*2*pi*n*f1/Fs) + exp(j*2*pi*n*f2/Fs); - papr = calc_papr(tx); - assert(abs(papr-3.0) < 0.05, 'test_papr() failed!') -end - -% "Genie" OFDM modem simulation that assumes ideal sync - -function [ber papr] = run_sim(Nc, Nsym, EbNodB, channel='awgn', plot_en=0, filt_en=0, method="", threshold=1, norm_ebno=0) - rand('seed',1); - randn('seed',1); - - M = 160; % number of samples in each symbol - bps = 2; % two bits per symbol for QPSK - Ncp = 16; % cyclic prefix samples - Fs = 8000; - - phase_est = 1; % perform phase estimation/correction - timing = Ncp; - - if strcmp(method,"diversity") - % total power of tx symbol after combination the same. Scatter plot positions - % different but also twice as much noise (bandwidth) - Nd = 2; gain = 1/sqrt(2); - else - Nd = 1; gain = 1.0; - end - - if strcmp(channel,'multipath') - dopplerSpreadHz = 1; path_delay = Ncp/2; - Nsam = floor(Nsym*(M+Ncp)*1.1); - spread1 = doppler_spread(dopplerSpreadHz, Fs, Nsam); - spread2 = doppler_spread(dopplerSpreadHz, Fs, Nsam); - end - - papr_log = []; - for e=1:length(EbNodB) - % generate a 2D array of QPSK symbols - - Nphases = 2^bps; - tx_phases = pi/2*floor((rand(Nsym,Nc)*Nphases)); - if strcmp(method,"diversity") - % duplicate carriers but with opposite phase - tx_phases = [tx_phases (tx_phases-pi/2)]; - end - tx_sym = gain*exp(j*tx_phases); - - % carrier frequencies, centre about 0 - st = floor(Nc*Nd/2); - w = 2*pi/M*(-st:-st+Nc*Nd-1); - - % generate OFDM signal - - tx = []; - for s=1:Nsym - atx = zeros(1,M); - for c=1:Nc*Nd - atx += exp(j*(0:M-1)*w(c))*tx_sym(s,c); - end - % insert cyclic prefix and build up stream of time domain symbols - % note CP costs us 10*log10((Ncp+M)/M) in Eb, as energy in CP isn't used for demodulation - tx = [tx atx(end-Ncp+1:end) atx]; - end - Nsam = length(tx); - - if strcmp(channel,'multipath') - assert(length(spread1) >= Nsam); - assert(length(spread2) >= Nsam); - end - - % bunch of PAPR reduction options - tx_ = tx; - - % determine threshold based on CDF - cdf = empirical_cdf((1:Nc),abs(tx)); - if strcmp(method, "clip") || strcmp(method, "diversity") || strcmp(method, "compand") - if threshold < 1 - threshold_level = find(cdf >= threshold)(1); - else - threshold_level = 10*Nc; - end - - % printf("threshold: %f threshold_level: %f\n", threshold, threshold_level); - end - - if strcmp(method, "clip") || strcmp(method, "diversity") - ind = find(abs(tx) > threshold_level); - tx_(ind) = threshold_level*exp(j*angle(tx(ind))); - end - if strcmp(method, "compand") - # power law compander x = a*y^power, y = (x/a) ^ (1/power) - power=2; a=threshold_level/(threshold_level^power); - tx_mag = (abs(tx)/a) .^ (1/power); - tx_ = tx_mag.*exp(j*angle(tx)); - end - - if filt_en - Nfilt=80; - b = fir1(Nfilt,2*Nc*Nd/M); - tx_ = filter(b,1,[tx_ zeros(1,Nfilt/2)]); - tx_ = [tx_(Nfilt/2+1:end)]; - end - - rx = tx_; - - % multipath channel - - if phase_est - % estimate phase of each symbol before multipath simulation - - rx_phase1 = zeros(Nsym,Nc); - for s=1:Nsym - st = (s-1)*(M+Ncp)+1+timing; en = st+M-1; - for c=1:Nc*Nd - rx_phase1(s,c) = sum(exp(-j*(0:M-1)*w(c)) .* rx(st:en))/M; - end - end - end - - if strcmp(channel,'multipath') - rx = spread1(1:Nsam).*rx + spread2(1:Nsam).*[zeros(1,path_delay) rx(1:end-path_delay)]; - end - - % normalise power after multipath, so that Eb/No is set up - % correctly - - if norm_ebno == 0 - norm = sqrt(mean(abs(tx_).^2)/mean(abs(rx).^2)); - else - % normalise after clipper, this makes norm_pwr constant for all test - % conditions - norm = sqrt(mean(abs(tx).^2)/mean(abs(rx).^2)); - end - rx *= norm; - norm_pwr = 10*log10(mean(abs(rx).^2)); - - if phase_est - % auxiliary rx to get ideal phase ests on signal after multipath but before AWGN noise is added - - rx_phase = zeros(Nsym,Nc); - for s=1:Nsym - st = (s-1)*(M+Ncp)+1+timing; en = st+M-1; - for c=1:Nc*Nd - arx_sym = sum(exp(-j*(0:M-1)*w(c)) .* rx(st:en))/M; - rx_phase(s,c) = arx_sym * conj(rx_phase1(s,c)); - end - end - rx_phase = exp(j*arg(rx_phase)); - end - - % AWGN channel - - EsNodB = EbNodB(e) + 10*log10(bps); - variance = M/(10^(EsNodB/10)); - noise = sqrt(variance/2)*randn(1,Nsam) + j*sqrt(variance/2)*randn(1,Nsam); - rx += noise; - - % demodulate - rx_sym = zeros(Nsym,Nc); - for s=1:Nsym - st = (s-1)*(M+Ncp)+1+timing; en = st+M-1; - for c=1:Nc*Nd - rx_sym(s,c) = sum(exp(-j*(0:M-1)*w(c)) .* rx(st:en))/M; - if phase_est rx_sym(s,c) *= conj(rx_phase(s,c)); end - end - - if strcmp(method,"diversity") - for c=1:Nc - rx_sym(s,c) += rx_sym(s,c+Nc)*exp(j*pi/2); - end - end - end - - % count bit errors - - Tbits = Terrs = 0; ErrPerSym = zeros(1,Nsym); - for s=1:Nsym - Nerrs = 0; - for c=1:Nc - tx_bits = qpsk_demod(tx_sym(s,c)); - rx_bits = qpsk_demod(rx_sym(s,c)); - Tbits += bps; - Nerrs += sum(xor(tx_bits,rx_bits)); - end - ErrPerSym(s) = Nerrs; - Terrs += Nerrs; - end - - if plot_en - figure(1); clf; - plot(abs(tx(1:5*M))); hold on; plot(abs(tx_(1:5*M))); hold off; - axis([0 5*M 0 max(abs(tx))]); - figure(2); clf; [hh nn] = hist(abs(tx),25,1); - plotyy(nn,hh,1:Nc,cdf); title('PDF and CDF'); grid; - figure(3); clf; plot(real(rx_sym(:,1:Nc)), imag(rx_sym(:,1:Nc)), '+'); axis([-2 2 -2 2]); - figure(4); clf; Tx_ = 10*log10(abs(fft(tx_))); plot(fftshift(Tx_)); - mx = 10*ceil(max(Tx_)/10); axis([1 length(Tx_) mx-60 mx]); - figure(5); plot_specgram(real(rx.*exp(j*2*pi*(0:Nsam-1)/4))); - figure(6); clf; stem(ErrPerSym); - end - - papr1 = calc_papr(tx); - papr2 = calc_papr(tx_); - papr_log = [papr_log papr2]; - ber(e) = Terrs/Tbits; - printf("EbNodB: %4.1f %3.1f %4.1f PAPR: %5.2f %5.2f Tbits: %6d Terrs: %6d BER: %5.3f\n", - EbNodB(e), norm, norm_pwr, papr1, papr2, Tbits, Terrs, ber(e)) - end - - papr = mean(papr_log); -end - -% BER versus Eb/No curves ------------------------------------- - -% first pass at trying out a few different schemes -function curves_experiment1(Nc=8, channel='awgn', Nsym=1000, EbNodB=2:8) - - [ber1 papr1] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1); - [ber2 papr2] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.8); - [ber3 papr3] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.6); - [ber4 papr4] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "compand", threshold=0.6); - [ber5 papr5] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "diversity", threshold=0.6); - - figure(7); clf; - semilogy(EbNodB, ber1,sprintf('b+-;vanilla OFDM %3.1f;',papr1),'markersize', 10, 'linewidth', 2); hold on; - semilogy(EbNodB, ber2,sprintf('r+-;clip 0.8 %3.1f;',papr2),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber3,sprintf('g+-;clip 0.6 %3.1f;',papr3),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber4,sprintf('c+-;compand 0.6 %3.1f;',papr4),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber5,sprintf('bk+-;diversity 0.6 %3.1f;',papr5),'markersize', 10, 'linewidth', 2); - hold off; - axis([min(EbNodB) max(EbNodB) 1E-3 1E-1]); grid; - xlabel('Eb/No'); title(sprintf("%s Nc = %d", channel, Nc)) - fn = sprintf("papr_exp1_%s_BER_EbNo.png", channel); - print(fn,"-dpng"); - - figure(8); clf; - semilogy(EbNodB+papr1, ber1,sprintf('b+-;vanilla OFDM %3.1f;',papr1),'markersize', 10, 'linewidth', 2); hold on; - semilogy(EbNodB+papr2, ber2,sprintf('r+-;clip 0.8 %3.1f;',papr2),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr3, ber3,sprintf('g+-;clip 0.6 %3.1f;',papr3),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr4, ber4,sprintf('c+-;compand 0.6 %3.1f;',papr4),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr5, ber5,sprintf('bk+-;diversity 0.6 %3.1f;',papr5),'markersize', 10, 'linewidth', 2); - hold off; - xlabel('Peak Eb/No'); - axis([min(EbNodB)+papr2 max(EbNodB)+papr1 1E-3 1E-1]); grid; title(sprintf("%s Nc = %d", channel, Nc)) - fn = sprintf("papr_exp1_%s_BER_peakEbNo.png", channel); - print(fn,"-dpng"); -end - - -% vary threshold and plot BER v Eb/No curves -function curves_experiment2(Nc=8, channel='awgn', Nsym=1000, EbNodB=2:16) - - [ber1 papr1] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1); - [ber2 papr2] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.8); - [ber3 papr3] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.6); - [ber4 papr4] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.4); - [ber5 papr5] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.2); - [ber6 papr6] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "diversity", threshold=0.8); - - figure(7); clf; - semilogy(EbNodB, ber1,sprintf('b+-;vanilla OFDM %3.1f;',papr1),'markersize', 10, 'linewidth', 2); hold on; - semilogy(EbNodB, ber2,sprintf('r+-;clip 0.8 %3.1f;',papr2),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber3,sprintf('g+-;clip 0.6 %3.1f;',papr3),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber4,sprintf('c+-;clip 0.4 %3.1f;',papr4),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber5,sprintf('bk+-;clip 0.2 %3.1f;',papr5),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber6,sprintf('m+-;diversity 0.8 %3.1f;', papr6),'markersize', 10, 'linewidth', 2); - hold off; - axis([min(EbNodB) max(EbNodB) 1E-3 1E-1]); grid; - xlabel('Eb/No'); title(sprintf("%s Nc = %d", channel, Nc)) - fn = sprintf("papr_exp2_Nc%d_%s_BER_EbNo.png", Nc, channel); - print(fn,"-dpng"); - - figure(8); clf; - semilogy(EbNodB+papr1, ber1,sprintf('b+-;vanilla OFDM %3.1f;',papr1),'markersize', 10, 'linewidth', 2); hold on; - semilogy(EbNodB+papr2, ber2,sprintf('r+-;clip 0.8 %3.1f;',papr2),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr3, ber3,sprintf('g+-;clip 0.6 %3.1f;',papr3),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr4, ber4,sprintf('c+-;clip 0.4 %3.1f;',papr4),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr5, ber5,sprintf('bk+-;clip 0.2 %3.1f;',papr5),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr6, ber6,sprintf('m+-;diversity 0.8 %3.1f;', papr6),'markersize', 10, 'linewidth', 2); - hold off; - xlabel('Peak Eb/No'); - axis([min(EbNodB)+papr2 max(EbNodB)+papr1 1E-3 1E-1]); grid; title(sprintf("%s Nc = %d", channel, Nc)) - fn = sprintf("papr_exp2_Nc%d_%s_BER_peakEbNo.png", Nc, channel); - print(fn,"-dpng"); -end - -% PAPR against number of carriers Nc -function curves_experiment3(Nsym=3000) - - paper = zeros(1,32); - Nc = 2:2:32; - for i = 1:length(Nc) - aNc = Nc(i); - [aber apapr] = run_sim(aNc, Nsym, 100); - papr(aNc) = apapr; - end - - figure(9); clf; - plot(Nc, papr(Nc)); xlabel('Number of Carriers Nc'); ylabel('PAPR (dB)'); grid; - fn = sprintf("papr_exp3_Nc.png"); - print(fn,"-dpng"); -end - -% focus on diversity - vary threshold and plot BER v Eb/No curves -function curves_experiment4(Nc=8, channel='multipath', Nsym=3000, EbNodB=2:2:16) - - [ber1 papr1] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1); - [ber2 papr2] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "diversity", threshold=1); - [ber3 papr3] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "diversity", threshold=0.8); - [ber4 papr4] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "diversity", threshold=0.6); - - figure(7); clf; - semilogy(EbNodB, ber1,sprintf('b+-;vanilla OFDM %3.1f;',papr1),'markersize', 10, 'linewidth', 2); hold on; - semilogy(EbNodB, ber2,sprintf('r+-;diversity 1.0 %3.1f;',papr2),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber3,sprintf('g+-;diversity 0.8 %3.1f;',papr3),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber4,sprintf('c+-;diversity 0.6 %3.1f;',papr4),'markersize', 10, 'linewidth', 2); - hold off; - axis([min(EbNodB) max(EbNodB) 1E-3 1E-1]); grid; - xlabel('Eb/No'); title(sprintf("%s Nc = %d", channel, Nc)) - fn = sprintf("papr_exp4_Nc%d_%s_BER_EbNo.png", Nc, channel); - print(fn,"-dpng"); - - figure(8); clf; - semilogy(EbNodB+papr1, ber1,sprintf('b+-;vanilla OFDM %3.1f;',papr1),'markersize', 10, 'linewidth', 2); hold on; - semilogy(EbNodB+papr2, ber2,sprintf('r+-;diversity 1.0 %3.1f;',papr2),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr3, ber3,sprintf('g+-;diversity 0.8 %3.1f;',papr3),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB+papr4, ber4,sprintf('c+-;diversity 0.6 %3.1f;',papr4),'markersize', 10, 'linewidth', 2); - hold off; - xlabel('Peak Eb/No'); - axis([min(EbNodB)+papr4 max(EbNodB)+papr1 1E-3 1E-1]); grid; title(sprintf("%s Nc = %d", channel, Nc)) - fn = sprintf("papr_exp4_Nc%d_%s_BER_peakEbNo.png", Nc, channel); - print(fn,"-dpng"); -end - -% plot BER v Eb/No curves for clipping with normalised Eb/No after clipping -function curves_experiment5(Nc=8, channel='awgn', Nsym=1000, EbNodB=2:10) - - [ber1 papr1] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "", threshold=1, norm=1); - [ber2 papr2] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.8, norm=1); - [ber3 papr3] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.6, norm=1); - [ber4 papr4] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.4, norm=1); - [ber5 papr5] = run_sim(Nc, Nsym, EbNodB, channel, 0, filt_en=1, "clip", threshold=0.2, norm=1); - - figure(7); clf; - semilogy(EbNodB, ber1,sprintf('b+-;vanilla OFDM %3.1f;',papr1),'markersize', 10, 'linewidth', 2); hold on; - semilogy(EbNodB, ber2,sprintf('r+-;clip 0.8 %3.1f;',papr2),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber3,sprintf('g+-;clip 0.6 %3.1f;',papr3),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber4,sprintf('c+-;clip 0.4 %3.1f;',papr4),'markersize', 10, 'linewidth', 2); - semilogy(EbNodB, ber5,sprintf('bk+-;clip 0.2 %3.1f;',papr5),'markersize', 10, 'linewidth', 2); - hold off; - axis([min(EbNodB) max(EbNodB) 1E-3 1E-1]); grid; - xlabel('Eb/No'); title(sprintf("%s Nc = %d", channel, Nc)) - fn = sprintf("papr_exp5_Nc%d_%s_BER_EbNo.png", Nc, channel); - print(fn,"-dpng"); -end - -pkg load statistics; -more off; - -test_papr; - -% single point with lots of plots ----------- - -%run_sim(8, 1000, EbNo=100, channel='awgn', plot_en=1, filt_en=1); -%run_sim(8, 8, EbNo=100, channel='awgn', plot_en=1, filt_en=1, "diversity", threshold=0.8); -%run_sim(8, 1000, EbNo=10, channel='multipath', plot_en=1, filt_en=0, "diversity", threshold=5); -%curves_experiment2(Nc=16, 'awgn', Nsym=1000); -curves_experiment2(Nc=16,'multipath', Nsym=3000, EbNodB=2:2:16); -%curves_experiment3() -%curves_experiment4() -%curves_experiment5(Nc=16) diff --git a/octave/phase_noise.m b/octave/phase_noise.m deleted file mode 100644 index 37e933c..0000000 --- a/octave/phase_noise.m +++ /dev/null @@ -1,72 +0,0 @@ -% phase_noise.m -% David Nov 2019 - -% Close-in look at phase noise. Feed in a off-air sample file of a -% sine wave, extracts the phase noise contour and returns the Doppler -% spreading function that can be used to model the channel in -% simulations - -function spread_FsHz = phase_noise(file_name) - Fs = 8000; - s = load_raw(file_name); - % skip past wave header - s = [zeros(256,1); s(256:end)]; - S = abs(fft(s(1:Fs).*hanning(Fs))); - [mx mx_bin] = max(S); - ftone = mx_bin-1; - - figure(1); clf; - plot(20*log10(S(1:Fs/2))) - title('Input Spectrum'); - - % downshift to baseband and LPF. We just want the sinusoid with as little - % additive AWGN noise as possible - sbb = s' .* exp(-j*(1:length(s))*2*pi*ftone/Fs); - [b a] = cheby1(4, 1, 20/Fs); - sbb_lpf = filter(b,a,sbb); - - spread_fsHz = sbb_lpf; - - % estimate and remove fine freq offset, and HF phase noise - - st = Fs; en = 20*Fs; - phase = unwrap(angle(sbb_lpf(st:en))); - fine_freq = mean(phase(2:end) - phase(1:end-1)); - sbb_lpf_fine = sbb_lpf .* exp(-j*(1:length(sbb_lpf))*fine_freq); - phase = unwrap(angle(sbb_lpf_fine(st:en))); - - printf("length: %3.2fs freq: %5.1f\n", length(s)/Fs, ftone+fine_freq*Fs/(2*pi)); - - figure(2); clf; - plot3((st:en)/Fs, real(sbb_lpf_fine(st:en)),imag(sbb_lpf_fine(st:en))) - title('Polar phase trajectory'); - - figure(3); clf; - S2 = fftshift(fft(sbb_lpf_fine(Fs:Fs*11))); - [mx mx_bin] = max(abs(S2)); - S2dB = 20*log10(abs(S2)); - mxdB = 10*ceil(max(S2dB)/10); - x = -10:0.1:10; - plot(x,S2dB(mx_bin-100:mx_bin+100)); - axis([-10 10 mxdB-40 mxdB]) - title('Close in Phase Noise Spectrum'); - xlabel('Freq (Hz)'); - grid; - - figure(5); clf; - t = (st:en)/Fs; - plot(t, phase,'b;phase;'); - title('Unwrapped Phase'); - xlabel('Time (sec)') - ylabel('Phase (radians)') - - figure(6); clf; - beta = 0.00001; - rate_of_change_Hz = filter(beta, [1 -(1-beta)],phase(2:end) - phase(1:end-1))*Fs/pi; - plot(t(2:end), rate_of_change_Hz) - title('Rate of change of phase (Hz)'); - xlabel('Time (sec)') - ylabel('Freq (Hz)') - - spread_FsHz = sbb_lpf_fine/std(sbb_lpf_fine); -end diff --git a/octave/pitch_test.m b/octave/pitch_test.m deleted file mode 100644 index 3fe0d1a..0000000 --- a/octave/pitch_test.m +++ /dev/null @@ -1,39 +0,0 @@ -% pitch_test.m -% David Rowe Sep 2009 -% Constructs a sequence to test the pitch estimator - -function pitch_test(samname) - M=320; - F=200; - - fs=fopen(samname,"wb"); - - f0 = 100; - for f=1:200 - Wo=2*pi*f0/8000; - P=2*pi/Wo; - L = floor(pi/Wo); - A = 10000/L; - phi = zeros(1,L); - s = zeros(1,M); - - for m=1:L - s = s + A*cos(m*Wo*(0:(M-1)) + phi(m)); - endfor - - figure(1); - clf; - plot(s); - - fwrite(fs,s,"short"); - - f0 = f0 + 5; - if (f0 > 400) - f0 = 100; - endif - endfor - - fclose(fs); - -endfunction - diff --git a/octave/pl.m b/octave/pl.m deleted file mode 100644 index 0d54788..0000000 --- a/octave/pl.m +++ /dev/null @@ -1,45 +0,0 @@ -% Copyright David Rowe 2009 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Plots a raw speech sample file, you can optionally specify the start and end -% samples and create a large and small PNGs - -function pl(samname1, start_sam, end_sam, pngname) - - fs=fopen(samname1,"rb"); - s=fread(fs,Inf,"short"); - - st = 1; - en = length(s); - if (nargin >= 2) - st = start_sam; - endif - if (nargin >= 3) - en = end_sam; - endif - - figure(1); - clf; - plot(s(st:en)); - axis([1 en-st 1.1*min(s) 1.1*max(s)]); - - if (nargin == 4) - - % small image - - __gnuplot_set__ terminal png size 420,300 - ss = sprintf("__gnuplot_set__ output \"%s.png\"", pngname); - eval(ss) - replot; - - % larger image - - __gnuplot_set__ terminal png size 800,600 - ss = sprintf("__gnuplot_set__ output \"%s_large.png\"", pngname); - eval(ss) - replot; - - endif - -endfunction diff --git a/octave/pl2.m b/octave/pl2.m deleted file mode 100644 index a25c6a2..0000000 --- a/octave/pl2.m +++ /dev/null @@ -1,44 +0,0 @@ -% Copyright David Rowe 2009 -% This program is distributed under the terms of the GNU General Public License -% Version 2 - -function pl2(samname1, samname2, start_sam, end_sam, offset) - - fs1=fopen(samname1,"rb"); - s1=fread(fs1,Inf,"short"); - fs2=fopen(samname2,"rb"); - s2=fread(fs2,Inf,"short"); - - st1 = st2 = 1; - en1 = en2 = length(s1); - if (nargin >= 3) - st1 = st2 = start_sam; - endif - if (nargin >= 4) - en1 = en2 = end_sam; - endif - - if (nargin == 5) - st2 += offset - en2 += offset - endif - - figure(1); - clf; - subplot(211); - l1 = strcat("r;",samname1,";"); - plot(s1(st1:en1), l1); grid minor; - axis([1 en1-st1 min(s1(st1:en1)) max(s1(st1:en1))]); - subplot(212); - l2 = strcat("r;",samname2,";"); - plot(s2(st2:en2),l2); grid minor; - axis([1 en2-st2 min(s1(st2:en2)) max(s1(st2:en2))]); - - figure(2) - plot(s1(st1:en1)-s2(st2:en2)); grid minor; - - f=fopen("diff.raw","wb"); - d = s1(st1:en1)-s2(st2:en2); - fwrite(f,d,"short"); - -endfunction diff --git a/octave/plamp.m b/octave/plamp.m deleted file mode 100644 index 1d575b8..0000000 --- a/octave/plamp.m +++ /dev/null @@ -1,135 +0,0 @@ -% Copyright David Rowe 2009 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Plot ampltiude modelling information from dump files. - -function plamp(samname, f, samname2) - - % switch some stuff off to unclutter display - - plot_sw = 0; - - sn_name = strcat(samname,"_sn.txt"); - Sn = load(sn_name); - - sw_name = strcat(samname,"_sw.txt"); - Sw = load(sw_name); - - sw__name = strcat(samname,"_sw_.txt"); - if (file_in_path(".",sw__name)) - Sw_ = load(sw__name); - endif - - ew_name = strcat(samname,"_ew.txt"); - if (file_in_path(".",ew_name)) - Ew = load(ew_name); - endif - - rk_name = strcat(samname,"_rk.txt"); - if (file_in_path(".",rk_name)) - Rk = load(rk_name); - endif - - model_name = strcat(samname,"_model.txt"); - model = load(model_name); - - modelq_name = strcat(samname,"_qmodel.txt"); - if (file_in_path(".",modelq_name)) - modelq = load(modelq_name); - endif - - pw_name = strcat(samname,"_pw.txt"); - if (file_in_path(".",pw_name)) - Pw = load(pw_name); - endif - - lsp_name = strcat(samname,"_lsp.txt"); - if (file_in_path(".",lsp_name)) - lsp = load(lsp_name); - endif - - phase_name = strcat(samname,"_phase.txt"); - if (file_in_path(".",phase_name)) - phase = load(phase_name); - endif - - phase_name_ = strcat(samname,"_phase_.txt"); - if (file_in_path(".",phase_name_)) - phase_ = load(phase_name_); - endif - - snr_name = strcat(samname,"_snr.txt"); - if (file_in_path(".",snr_name)) - snr = load(snr_name); - endif - - % optional second file, for exploring post filter - - if nargin == 3 - model2_name = strcat(samname2,"_model.txt"); - model2 = load(model2_name); - sn2_name = strcat(samname2,"_sn.txt"); - Sn2 = load(sn2_name); - - sw_name2 = strcat(samname2,"_sw.txt"); - Sw2 = load(sw_name2); - end - - k = ' '; - do - figure(1); - clf; - s = [ Sn(2*f-1,:) Sn(2*f,:) ]; - plot(s,'b'); - if (nargin == 3) - s2 = [ Sn2(2*f-1,:) Sn2(2*f,:) ]; - hold on; plot(s2,'r'); hold off; - end - axis([1 length(s) -30000 30000]); - - figure(2); - Wo = model(f,1); - L = model(f,2); - Am = model(f,3:(L+2)); - plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;+-b"); - axis([1 4000 -10 80]); - hold on; - if plot_sw; plot((0:255)*4000/256, Sw(f,:),";Sw;b"); end - - if (nargin == 3) - Wo2 = model2(f,1); - L2 = model2(f,2); - Am2 = model2(f,3:(L2+2)); - plot((1:L2)*Wo2*4000/pi, 20*log10(Am2),";Am2;+-r" ); - if plot_sw; plot((0:255)*4000/256, Sw2(f,:),";Sw2;r"); end - endif - - hold off; grid minor; - - % interactive menu - - printf("\rframe: %d menu: n-next b-back p-png s-plot_sw q-quit", f); - fflush(stdout); - k = kbhit(); - if k == 'n'; f = f + 1; endif - if k == 'b'; f = f - 1; endif - if k == 's' - if plot_sw; plot_sw = 0; else; plot_sw = 1; end - endif - % optional print to PNG - - if (k == 'p') - figure(1); - pngname = sprintf("%s_%d_sn.png",samname,f); - print(pngname, '-dpng', "-S800,600") - - figure(2); - pngname = sprintf("%s_%d_sw.png",samname,f); - print(pngname, '-dpng', "-S800,600") - endif - - until (k == 'q') - printf("\n"); - -endfunction diff --git a/octave/plinterp.m b/octave/plinterp.m deleted file mode 100644 index 794a085..0000000 --- a/octave/plinterp.m +++ /dev/null @@ -1,11 +0,0 @@ -load ../unittest/tinterp_prev.txt; -load ../unittest/tinterp_interp.txt; -load ../unittest/tinterp_next.txt; - -clf; -plot(tinterp_prev(:,1), 20.0*log10(tinterp_prev(:,2)),";prev;") -hold on; -plot(tinterp_interp(:,1), 20.0*log10(tinterp_interp(:,2)),'g+-;interp;') -plot(tinterp_next(:,1), 20.0*log10(tinterp_next(:,2)),'ro-;next;') -hold off; -axis([0 pi 0 80]) diff --git a/octave/pllpcpf.m b/octave/pllpcpf.m deleted file mode 100644 index 924e045..0000000 --- a/octave/pllpcpf.m +++ /dev/null @@ -1,150 +0,0 @@ -% Copyright David Rowe 2012 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Plot amplitude modelling information from dump files to test and develop -% LPC post filter. - -function pllpcpf(samname, f) - - % switch some stuff off to unclutter display - - plot_Am = 0; - plot_Amq = 0; - plot_err = 0; - plot_lsp = 0; - plot_snr = 0; - plot_vsnr = 0; - plot_sw = 0; - plot_pw = 1; - plot_pwb = 1; - plot_rw = 1; - - sn_name = strcat(samname,"_sn.txt"); - Sn = load(sn_name); - - sw_name = strcat(samname,"_sw.txt"); - Sw = load(sw_name); - - sw__name = strcat(samname,"_sw_.txt"); - if (file_in_path(".",sw__name)) - Sw_ = load(sw__name); - endif - - model_name = strcat(samname,"_model.txt"); - model = load(model_name); - - modelq_name = strcat(samname,"_qmodel.txt"); - if (file_in_path(".",modelq_name)) - modelq = load(modelq_name); - endif - - % Pw (LPC synth filter spectrum) before post filter - - pwb_name = strcat(samname,"_pwb.txt"); - if (file_in_path(".",pwb_name)) - Pwb = load(pwb_name); - endif - - % Rw (Post filter spectrum) - - rw_name = strcat(samname,"_rw.txt"); - if (file_in_path(".",rw_name)) - Rw = load(rw_name); - endif - - % Pw (LPC synth filter spectrum) after post filter - - pw_name = strcat(samname,"_pw.txt"); - if (file_in_path(".",pw_name)) - Pw = load(pw_name); - endif - - - Ew_on = 1; - k = ' '; - do - figure(1); - clf; - s = [ Sn(2*f-1,:) Sn(2*f,:) ]; - size(s); - plot(s); - axis([1 length(s) -20000 20000]); - - figure(2); - clf; - Wo = model(f,1); - L = model(f,2); - Am = model(f,3:(L+2)); - if plot_Am - plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;r"); - end - axis([1 4000 -10 80]); - hold on; - if plot_sw - plot((0:255)*4000/256, Sw(f,:),";Sw;"); - end - - if (file_in_path(".",modelq_name)) - - Amq = modelq(f,3:(L+2)); - if plot_Amq - plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" ); - end - - if (file_in_path(".",pwb_name) && plot_pwb) - plot((0:255)*4000/256, 10*log10(Pwb(f,:)),";Pwb;r"); - endif - - if (file_in_path(".",rw_name) && plot_rw) - plot((0:255)*4000/256, 10*log10(Rw(f,:)),";Rw;b"); - endif - - if (file_in_path(".",pw_name) && plot_pw) - plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;g."); - endif - - signal = Am * Am'; - noise = (Am-Amq) * (Am-Amq)'; - snr1 = 10*log10(signal/noise); - Am_err_label = sprintf(";Am error SNR %4.2f dB;m",snr1); - if plot_err - plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label); - end - endif - - - hold off; - - % interactive menu - - printf("\rframe: %d menu: n-next b-back p-png q-quit", f); - fflush(stdout); - k = kbhit(); - if (k == 'n') - f = f + 1; - endif - if (k == 'b') - f = f - 1; - endif - - % optional print to PNG - - if (k == 'p') - figure(1); - pngname = sprintf("%s_%d_sn.png",samname,f); - print(pngname, '-dpng', "-S500,500") - pngname = sprintf("%s_%d_sn_large.png",samname,f); - print(pngname, '-dpng', "-S800,600") - - figure(2); - pngname = sprintf("%s_%d_sw.png",samname,f); - print(pngname, '-dpng', "-S500,500") - pngname = sprintf("%s_%d_sw_large.png",samname,f); - print(pngname, '-dpng', "-S1200,800") - endif - - until (k == 'q') - printf("\n"); - -endfunction diff --git a/octave/pllsp.m b/octave/pllsp.m deleted file mode 100644 index 0606d3c..0000000 --- a/octave/pllsp.m +++ /dev/null @@ -1,46 +0,0 @@ -% Copyright David Rowe 2010 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Plots a bunch of information related to LSP quantisation: -% - speech file -% - LSPs before and after quantisation -% - SNR for each frame -% -% Note: there is a 160 sample (two frame delay) from the when a sample -% enters the input buffer until it is at the centre of the analysis window - -function pllsp(rawfile, - dumpfile_prefix_lpc_only, - dumpfile_prefix_lsp, - start_f, end_f) - - fs=fopen(rawfile,"rb"); - s=fread(fs,Inf,"short"); - - lpc_snr_name = strcat(dumpfile_prefix_lpc_only,"_lpc_snr.txt"); - lpc10_snr = load(lpc_snr_name); - lpc_snr_name = strcat(dumpfile_prefix_lsp,"_lpc_snr.txt"); - lsp_snr = load(lpc_snr_name); - - lsp_name = strcat(dumpfile_prefix_lsp,"_lsp.txt"); - lsps = load(lsp_name); - [m,n]=size(lsps); - lsp = lsps(1:2:m,:); - lsp_ = lsps(2:2:m,:); - - figure(1); - clf; - subplot(211); - sp = s((start_f-2)*80:(end_f-2)*80); - plot(sp); - - subplot(212); - plot(lpc10_snr((start_f+1):end_f)-lsp_snr((start_f+1):end_f)); - - figure(2); - plot((4000/pi)*lsp((start_f+1):end_f,:)); - hold on; - plot((4000/pi)*lsp_((start_f+1):end_f,:),'+-'); - hold off; -endfunction diff --git a/octave/pllspdt.m b/octave/pllspdt.m deleted file mode 100644 index c711aa4..0000000 --- a/octave/pllspdt.m +++ /dev/null @@ -1,27 +0,0 @@ -% pllspdt.m -% Copyright David Rowe 2010 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Test script to plot differences in LSps between frames - -function pllspdt(rawfile,dumpfile_prefix_lsp,lspn, start_f, end_f) - - fs=fopen(rawfile,"rb"); - s=fread(fs,Inf,"short"); - - lsp_name = strcat(dumpfile_prefix_lsp,"_lsp.txt"); - lsps = load(lsp_name); - [m,n]=size(lsps); - lsp = lsps(1:2:m,:); - lsp_ = lsps(2:2:m,:); - lspdt = lsp(2:m/2,:) - lsp(1:m/2-1,:); - - figure(1); - clf; - sp = s((start_f-2)*80:(end_f-2)*80); - plot(sp); - - figure(2); - plot((4000/pi)*lspdt((start_f+1):end_f,lspn)); -endfunction diff --git a/octave/plnlp.m b/octave/plnlp.m deleted file mode 100644 index 01b4931..0000000 --- a/octave/plnlp.m +++ /dev/null @@ -1,134 +0,0 @@ -% Copyright David Rowe 2009 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Plot NLP states from dump files. - -function plnlp(samname, f) - - sn_name = strcat(samname,"_sn.txt"); - Sn = load(sn_name); - - sw_name = strcat(samname,"_sw.txt"); - Sw = load(sw_name); - - fw_name = strcat(samname,"_fw.txt"); - if (file_in_path(".",fw_name)) - fw = load(fw_name); - endif - - e_name = strcat(samname,"_e.txt"); - if (file_in_path(".",e_name)) - e = load(e_name); - endif - - p_name = strcat(samname,".p"); - if (file_in_path(".",p_name)) - p = load(p_name); - endif - - sq_name = strcat(samname,"_sq.txt"); - if (file_in_path(".",sq_name)) - sq = load(sq_name); - endif - - dec_name = strcat(samname,"_dec.txt"); - if (file_in_path(".",dec_name)) - dec = load(dec_name); - endif - - do - figure(1); - clf; - s = [ Sn(2*f-1,:) Sn(2*f,:) ]; - plot(s, ";Sn;"); - grid - axis([1 length(s) -20000 20000]); - - figure(2); - plot((0:255)*4000/256, Sw(f,:),";Sw;"); - grid - axis([1 4000 -10 80]); - hold on; - - f0 = 8000/p(f); - Wo = 2*pi/p(f); - L = floor(pi/Wo); - f0_label = sprintf("b;P=%3.1f F0=%3.0f;",p(f),f0); - for m=1:L-1 - plot([ m*Wo*4000/pi m*Wo*4000/pi], [10 60], 'b'); - endfor - plot([ L*Wo*4000/pi L*Wo*4000/pi], [10 60], f0_label); - - hold off; - - if (file_in_path(".",fw_name)) - figure(3); - if (file_in_path(".",e_name)) - subplot(211); - endif - plot((0:255)*800/256, fw(f,:)/max(fw(f,:)), ";Fw;"); - axis([1 400 0 1]); - if (file_in_path(".",e_name)) - subplot(212); - e_concat = [ e(2*f-1,:) e(2*f,:) ]; - plot(e_concat(1:400)/max(e_concat(1:400)), "+;MBE E(f);"); - axis([1 400 0 1]); - endif - endif - - if (file_in_path(".",sq_name)) - figure(4); - sq_concat = [ sq(2*f-1,:) sq(2*f,:) ]; - axis - plot(sq_concat, ";sq;"); - endif - - if (file_in_path(".",dec_name)) - figure(5); - plot(dec(f,:), ";dec;"); - endif - - figure(2); - - % interactive menu - - printf("\rframe: %d menu: n-next b-back p-png q-quit ", f); - fflush(stdout); - k = kbhit(); - if (k == 'n') - f = f + 1; - endif - if (k == 'b') - f = f - 1; - endif - - % optional print to PNG - - if (k == 'p') - - pngname = sprintf("%s_%d",samname,f); - - % small image - - __gnuplot_set__ terminal png size 420,300 - ss = sprintf("__gnuplot_set__ output \"%s.png\"", pngname); - eval(ss) - replot; - - % larger image - - __gnuplot_set__ terminal png size 800,600 - ss = sprintf("__gnuplot_set__ output \"%s_large.png\"", pngname); - eval(ss) - replot; - - % for some reason I need this to stop large plot getting wiped - __gnuplot_set__ output "/dev/null" - - endif - - until (k == 'q') - printf("\n"); - -endfunction diff --git a/octave/plphase.m b/octave/plphase.m deleted file mode 100644 index c12422e..0000000 --- a/octave/plphase.m +++ /dev/null @@ -1,198 +0,0 @@ -% Copyright David Rowe 2009 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Plot phase modelling information from dump files. - -function plphase(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); - - sw__name = strcat(samname,"_sw_.txt"); - if (file_in_path(".",sw__name)) - Sw_ = load(sw__name); - endif - - pw_name = strcat(samname,"_pw.txt"); - if (file_in_path(".",pw_name)) - Pw = load(pw_name); - endif - - ak_name = strcat(samname,"_ak.txt"); - if (file_in_path(".",ak_name)) - ak = load(ak_name); - endif - - phase_name = strcat(samname,"_phase.txt"); - if (file_in_path(".",phase_name)) - phase = load(phase_name); - endif - - phase_name_ = strcat(samname,"_phase_.txt"); - if (file_in_path(".",phase_name_)) - phase_ = load(phase_name_); - endif - - snr_name = strcat(samname,"_snr.txt"); - if (file_in_path(".",snr_name)) - snr = load(snr_name); - endif - - sn_name_ = strcat(samname,".raw"); - if (file_in_path(".",sn_name_)) - fs_ = fopen(sn_name_,"rb"); - sn_ = fread(fs_,Inf,"short"); - endif - - k = ' '; - do - figure(1); - clf; - s = [ Sn(2*f-1,:) Sn(2*f,:) ]; - plot(s); - grid; - axis([1 length(s) -20000 20000]); - if (k == 'p') - pngname = sprintf("%s_%d_sn",samname,f); - png(pngname); - endif - - figure(2); - Wo = model(f,1); - L = model(f,2); - Am = model(f,3:(L+2)); - plot((1:L)*Wo*4000/pi, 20*log10(Am),"r;Am;"); - axis([1 4000 -10 80]); - hold on; - plot((0:255)*4000/256, Sw(f,:),";Sw;"); - grid; - - if (file_in_path(".",sw__name)) - plot((0:255)*4000/256, Sw_(f,:),"g;Sw_;"); - endif - - if (file_in_path(".",pw_name)) - plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;"); - endif - - if (file_in_path(".",snr_name)) - snr_label = sprintf(";phase SNR %4.2f dB;",snr(f)); - plot(1,1,snr_label); - endif - - % phase model - determine SNR and error spectrum for phase model 1 - - if (file_in_path(".",phase_name_)) - orig = Am.*exp(j*phase(f,1:L)); - synth = Am.*exp(j*phase_(f,1:L)); - signal = orig * orig'; - noise = (orig-synth) * (orig-synth)'; - snr_phase = 10*log10(signal/noise); - - phase_err_label = sprintf("g;phase_err SNR %4.2f dB;",snr_phase); - plot((1:L)*Wo*4000/pi, 20*log10(orig-synth), phase_err_label); - endif - - hold off; - if (k == 'p') - pngname = sprintf("%s_%d_sw",samname,f); - png(pngname); - endif - - if (file_in_path(".",phase_name)) - figure(3); - plot((1:L)*Wo*4000/pi, phase(f,1:L)*180/pi, "-o;phase;"); - axis; - if (file_in_path(".", phase_name_)) - hold on; - plot((1:L)*Wo*4000/pi, phase_(f,1:L)*180/pi, "g;phase after;"); - grid - hold off; - endif - if (k == 'p') - pngname = sprintf("%s_%d_phase",samname,f); - png(pngname); - endif - endif - - % synthesised speech - - if (file_in_path(".",sn_name_)) - figure(4); - s_ = sn_((f-3)*80+1:(f+1)*80); - plot(s_); - axis([1 length(s_) -20000 20000]); - if (k == 'p') - pngname = sprintf("%s_%d_sn_",samname,f) - png(pngname); - endif - endif - - if (file_in_path(".",ak_name)) - figure(5); - axis; - akw = ak(f,:); - weight = 1.0 .^ (0:length(akw)-1); - akw = akw .* weight; - H = 1./fft(akw,8000); - subplot(211); - plot(20*log10(abs(H(1:4000))),";LPC mag spec;"); - grid; - subplot(212); - plot(angle(H(1:4000))*180/pi,";LPC phase spec;"); - grid; - if (k == 'p') - % stops multimode errors from gnuplot, I know not why... - figure(2); - figure(5); - - pngname = sprintf("%s_%d_lpc",samname,f); - png(pngname); - endif - endif - - - % autocorrelation function to research voicing est - - %M = length(s); - %sw = s .* hanning(M)'; - %for k=0:159 - % R(k+1) = sw(1:320-k) * sw(1+k:320)'; - %endfor - %figure(4); - %R_label = sprintf(";R(k) %3.2f;",max(R(20:159))/R(1)); - %plot(R/R(1),R_label); - %grid - - figure(2); - - % interactive menu - - printf("\rframe: %d menu: n-next b-back p-png q-quit ", f); - fflush(stdout); - k = kbhit(); - if (k == 'n') - f = f + 1; - endif - if (k == 'b') - f = f - 1; - endif - - % optional print to PNG - - if (k == 'p') - pngname = sprintf("%s_%d",samname,f); - png(pngname); - endif - - until (k == 'q') - printf("\n"); - -endfunction diff --git a/octave/plpitch.m b/octave/plpitch.m deleted file mode 100644 index 69ad533..0000000 --- a/octave/plpitch.m +++ /dev/null @@ -1,36 +0,0 @@ -% Copyright David Rowe 2009 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% plpitch.m -% Plots two pitch tracks on top of each other, used for comparing pitch -% estimators - -function plpitch(pitch1_name, pitch2_name, start_fr, end_fr) - - pitch1 = load(pitch1_name); - pitch2 = load(pitch2_name); - - st = 1; - en = length(pitch1); - if (nargin >= 3) - st = start_fr; - endif - if (nargin >= 4) - en = end_fr; - endif - - figure(1); - clf; - l1 = strcat("r;",pitch1_name,";") - l1 - st - en - plot(pitch1(st:en), l1); - axis([1 en-st 20 160]); - l2 = strcat("g;",pitch2_name,";"); - hold on; - plot(pitch2(st:en),l2); - hold off; -endfunction - diff --git a/octave/plvoicing.m b/octave/plvoicing.m deleted file mode 100644 index a531747..0000000 --- a/octave/plvoicing.m +++ /dev/null @@ -1,89 +0,0 @@ -% Copyright David Rowe 2009 -% This program is distributed under the terms of the GNU General Public License -% Version 2 -% -% Plot voicing information from sample and dump files. -% -% samfilename is the raw source file, e.g. "../raw/hts1a.raw" -% samname is the dumpfile prefix, e.g. "../src/hts1a" -% -% There is a 160 sample (two frame delay) from the when a sample -% enters the input buffer until it is at the centre of the analysis window - -function plvoicing(samfilename, samname, start_f, end_f, pngname) - - fs=fopen(samfilename,"rb"); - s=fread(fs,Inf,"short"); - - snr_name = strcat(samname,"_snr.txt"); - snr = load(snr_name); - model_name = strcat(samname,"_model.txt"); - model = load(model_name); - - Wo = model((start_f+1):end_f,1); - F0 = Wo*4000/pi; - dF0 = F0(1:length(Wo)-1) - F0(2:length(Wo)); - - % work out LP and HP energy - - for f=(start_f+1):end_f - L = model(f,2); - Am = model(f,3:(L+2)); - L2 = floor(L/2); - elow = Am(1:L2) * Am(1:L2)'; - ehigh = Am(L2:L) * Am(L2:L)'; - erat(f-(start_f+1)+1) = 10*log10(elow/ehigh); - endfor - - figure(1); - clf; - sp = s((start_f-2)*80:(end_f-2)*80); - plot(sp); - hold on; - vhigh = snr((start_f+1):end_f) > 7; - vlow = snr((start_f+1):end_f) > 4; - - % test correction based on erat - - vlowadj = vlow; - - for f=1:length(erat)-1 - if (vlow(f) == 0) - if (erat(f) > 10) - vlowadj(f) = 1; - endif - endif - if (vlow(f) == 1) - if (erat(f) < -10) - vlowadj(f) = 0; - endif - if (abs(dF0(f)) > 15) - vlowadj(f) = 0; - endif - endif - endfor - - x = 1:(end_f-start_f); - plot(x*80,snr((start_f+1):end_f)*1000,';SNRdB x 1000;g+'); - plot(x*80,-8000 + vhigh*2000,';7dB thresh;g'); - plot(x*80,-11000 + vlowadj*2000,';vlow with corr;g'); - plot(x*80,erat*1000,';elow/ehigh in dB;r'); - plot(x*80,-14000 + vlow*2000,';4dB thresh;r'); - hold off; - grid - if (nargin == 5) - print(pngname, "-dpng", "-S500,500") - endif - - figure(2) - Wo = model((start_f+1):end_f,1); - F0 = Wo*4000/pi; - dF0 = F0(1:length(Wo)-1) - F0(2:length(Wo)); - %plot(dF0,'+--') - %hold on; - %plot([ 1 length(dF0) ], [10 10] ,'r') - %plot([ 1 length(dF0) ], [-10 -10] ,'r') - %axis([1 length(dF0) -50 50]) - %hold off; - plot(F0,'+--') -endfunction diff --git a/octave/power_from_stdio.m b/octave/power_from_stdio.m deleted file mode 100644 index a4dbddc..0000000 --- a/octave/power_from_stdio.m +++ /dev/null @@ -1,27 +0,0 @@ -% power_from_gr.m -% David Rowe June 2018 -% -% Measure power of signal from stdio, used for SNR tests from analog radios - -#{ - $ rec -t raw -r 8000 -s -2 -c 1 - -q | octave --no-gui -qf power_from_stdio.m -#} - -graphics_toolkit ("gnuplot") - -Fs = 48000; % sample rate in Hz -shorts_per_sample = 1; % real samples - -[s,c] = fread(stdin, shorts_per_sample*Fs, "short"); - -while c - S = fft(s.*hanning(Fs)); - SdB = 20*log10(abs(S)); - figure(1); plot(real(s)); axis([0 Fs -3E4 3E4]); - figure(2); plot(SdB); axis([0 12000 40 160]); - - printf("power: %f dB\n", 10*log10(var(s))); - %pause(2); - [s,c] = fread(stdin, shorts_per_sample*Fs, "short"); -endwhile - diff --git a/octave/pulse.m b/octave/pulse.m deleted file mode 100644 index 223389e..0000000 --- a/octave/pulse.m +++ /dev/null @@ -1,37 +0,0 @@ -% pulse.m -% David Rowe August 2009 -% -% Experiments with human pulse perception for sinusoidal codecs - -function pulse(samname) - - A = 1000; - K = 16000; - N = 80; - frames = K/N; - s = zeros(1,K); - - for f=1:frames - % lets try placing np random pulses in every frame - - P = 20 + (160-20)*rand(1,1); - Wo = 2*pi/P; - L = floor(pi/Wo); - sf = zeros(1,N); - for m=1:L/2:L - pos = floor(rand(1,1)*N)+1; - %pos = 50; - for l=m:m+L/2-1 - sf = sf + A*cos(l*Wo*((f-1)*N+1:f*N) - pos*l*Wo); - endfor - endfor - s((f-1)*N+1:f*N) = sf; - endfor - - plot(s(1:250)); - - fs=fopen(samname,"wb"); - fwrite(fs,s,"short"); - fclose(fs); -endfunction - diff --git a/octave/save_array_c_header.m b/octave/save_array_c_header.m deleted file mode 100644 index 1fd87bb..0000000 --- a/octave/save_array_c_header.m +++ /dev/null @@ -1,14 +0,0 @@ -% save_array_c_header.m -% -% David Rowe Sep 2015 - -function save_array_c_header(array, array_name, filename) - f=fopen(filename,"wt"); - fprintf(f,"/* Generated by save_array_c_header.m Octave function */\n\n"); - fprintf(f,"const float %s[]={\n", array_name); - for m=1:length(array)-1 - fprintf(f," % .16f,\n",array(m)); - endfor - fprintf(f," % .16f\n};\n",array(length(array))); - fclose(f); -endfunction diff --git a/octave/save_comp.m b/octave/save_comp.m deleted file mode 100644 index b8e663f..0000000 --- a/octave/save_comp.m +++ /dev/null @@ -1,12 +0,0 @@ -% save_comp.m -% David Rowe Aug 2020 - -function save_comp(fn, iq) - l = length(iq); - s = zeros(1,2*l); - s(1:2:2*l) = real(iq); - s(2:2:2*l) = imag(iq); - fs=fopen(fn,"wb"); - s = fwrite(fs,s,"float32"); - fclose(fs); -endfunction diff --git a/octave/save_f32.m b/octave/save_f32.m deleted file mode 100644 index 62214d4..0000000 --- a/octave/save_f32.m +++ /dev/null @@ -1,12 +0,0 @@ -% save_f32.m -% David Rowe Sep 2021 -% -% save a matrix to .f32 binary files in row-major order - -function save_f32(fn, m) - f=fopen(fn,"wb"); - [r c] = size(m); - mlinear = reshape(m', 1, r*c); - fwrite(f, mlinear, 'float32'); - fclose(f); -endfunction diff --git a/octave/save_hackrf.m b/octave/save_hackrf.m deleted file mode 100644 index 5dda704..0000000 --- a/octave/save_hackrf.m +++ /dev/null @@ -1,13 +0,0 @@ -% save_hackrf.m -% -% David Rowe Aug 2020 - -function save_hackrf(fn,iq) - l = length(iq); - s = zeros(1,2*l); - s(1:2:2*l) = real(iq); - s(2:2:2*l) = imag(iq); - fs = fopen(fn,"wb"); - fwrite(fs,s,"schar"); - fclose(fs); -endfunction diff --git a/octave/save_raw.m b/octave/save_raw.m deleted file mode 100644 index 7f17277..0000000 --- a/octave/save_raw.m +++ /dev/null @@ -1,7 +0,0 @@ -% save_raw.m -% David Rowe 9 Feb 2015 - -function s = save_raw(fn,s) - fs=fopen(fn,"wb"); - fwrite(fs,s,"short"); -endfunction diff --git a/octave/snr_curves_plot.m b/octave/snr_curves_plot.m deleted file mode 100644 index 8a2e234..0000000 --- a/octave/snr_curves_plot.m +++ /dev/null @@ -1,262 +0,0 @@ -% snr_curves_plot.m -% -% Companion script for unittest/raw_data_curves - -1; - -function state_vec = set_graphics_state_print() - textfontsize = get(0,"defaulttextfontsize"); - linewidth = get(0,"defaultlinelinewidth"); - markersize = get(0, "defaultlinemarkersize"); - set(0, "defaulttextfontsize", 16); - set(0, "defaultaxesfontsize", 16); - set(0, "defaultlinelinewidth", 1); - state_vec = [textfontsize linewidth markersize]; -endfunction - -function set_graphics_state_screen(state_vec) - textfontsize = state_vec(1); - linewidth = state_vec(2); - markersize = state_vec(3); - set(0, "defaulttextfontsize", textfontsize); - set(0, "defaultaxesfontsize", textfontsize); - set(0, "defaultlinelinewidth", linewidth); - set(0, "defaultlinemarkersize", markersize); -endfunction - -function [snr_ch per] = snr_scatter(source, mode, channel, colour) - suffix = sprintf("_%s_%s_%s",source, mode, channel); - snr = load(sprintf("snr%s.txt",suffix)); - offset = load(sprintf("offset%s.txt",suffix)); - snr -= offset; - snr_x = []; snrest_y = []; - for i=1:length(snr) - fn = sprintf('snrest%s_%d.txt',suffix,i); - if exist(fn,'file') == 2 - snrest=load(fn); - if i == length(snr) - plot(snr(i)*ones(1,length(snrest)), snrest, sprintf('%s;%s %s;',colour,source,mode)); - else - plot(snr(i)*ones(1,length(snrest)), snrest, sprintf('%s',colour)); - end - snr_x = [snr_x snr(i)]; snrest_y = [snrest_y mean(snrest)]; - end - end - plot(snr_x, snrest_y, sprintf('%s', colour)); -endfunction - -function [snr_ch per] = per_snr(mode, colour) - snrch = load(sprintf("snrch_%s.txt",mode)); - snroffset = load(sprintf("snroffset_%s.txt",mode)); - snrch -= snroffset; - per = load(sprintf("per_%s.txt",mode)); - plot(snrch, per, sprintf('%so-;%s;', colour, mode)); -endfunction - -function snrest_snr_screen(source, channel) - clf; hold on; - snr_scatter(source, 'datac0', channel,'b+-') - snr_scatter(source, 'datac1', channel,'g+-') - snr_scatter(source, 'datac3', channel,'r+-') - snr_scatter(source, 'datac4', channel,'c+-') - snr_scatter(source, 'datac13', channel,'m+-') - xlabel('SNR (dB)'); ylabel('SNRest (dB)'); grid('minor'); - axis([-12 12 -12 12]); - a = axis; - plot([a(1) a(2)],[a(1) a(2)],'bk-'); - hold off; grid; - if strcmp(source,'ctx') - title(sprintf('SNR estimate versus SNR (%s) (no compression)', channel)); - else - title(sprintf('SNR estimate versus SNR (%s) (with compression)', channel)); - end - legend('location','northwest'); -endfunction - -function snrest_snr_print(source, channel) - state_vec = set_graphics_state_print(); - snrest_snr_screen(source, channel); - print(sprintf("snrest_snr_%s.png", source), "-dpng", "-S1000,800"); - set_graphics_state_screen(state_vec); -endfunction - -function ber_per_v_snr(source, mode, channel, colour) - suffix = sprintf("_%s_%s_%s.txt",source, mode, channel); - snr = load(sprintf("snr%s",suffix)); - offset = load(sprintf("offset%s",suffix)); - snr -= offset; - ber = load(sprintf("ber%s",suffix)) + 1E-6; - per = load(sprintf("per%s",suffix)) + 1E-6; - semilogy(snr, ber, sprintf('%s;%s %s ber;', colour, source, mode)); - semilogy(snr, per, sprintf('%s;%s %s per;', colour, source, mode),'linewidth',3,'markersize',10); -endfunction - -function per_v_snr(source, mode, channel, colour) - suffix = sprintf("_%s_%s_%s.txt",source, mode, channel); - snr = load(sprintf("snr%s",suffix)); - offset = load(sprintf("offset%s",suffix)); - snr -= offset; - per = load(sprintf("per%s",suffix)) + 1E-6; - if strcmp(channel,"awgn") - semilogy(snr, per, sprintf('%s;%s %s;', colour, mode, channel)); - else - semilogy(snr, per, sprintf('%s;%s %s;', colour, mode, channel),'linewidth',3,'markersize',10); - end -endfunction - -function thruput_v_snr(source, mode, channel, colour) - suffix = sprintf("_%s_%s_%s.txt",source, mode, channel); - snr = load(sprintf("snr%s",suffix)); - offset = load(sprintf("offset%s",suffix)); - snr -= offset; - per = load(sprintf("per%s",suffix)) + 1E-6; - if strcmp(mode,"datac0") Rb=291; end; - if strcmp(mode,"datac1") Rb=980; end; - if strcmp(mode,"datac3") Rb=321; end; - if strcmp(mode,"datac4") Rb=87; end; - if strcmp(mode,"datac13") Rb=65; end; - if strcmp(channel,"awgn") - plot(snr, Rb*(1-per), sprintf('%s;%s %s;', colour, mode, channel)); - else - plot(snr, Rb*(1-per), sprintf('%s;%s %s;', colour, mode, channel),'linewidth',3,'markersize',10); - end -endfunction - -function octave_ch_noise_screen(channel) - clf; hold on; - ber_per_v_snr('oct','datac0',channel,'bo-') - ber_per_v_snr('ch' ,'datac0',channel,'bx-') - ber_per_v_snr('oct','datac1',channel,'go-') - ber_per_v_snr('ch' ,'datac1',channel,'gx-') - ber_per_v_snr('oct','datac3',channel,'ro-') - ber_per_v_snr('ch' ,'datac3',channel,'rx-') - xlabel('SNR (dB)'); grid; - hold off; - if strcmp(channel,"awgn") - axis([-6 8 1E-3 1]); - else - axis([-2 12 1E-3 1]); - end - title(sprintf('Comparsion of Measuring SNR from Octave and ch tool (%s)', channel)); -endfunction - -function octave_ch_noise_print(channel) - state_vec = set_graphics_state_print(); - octave_ch_noise_screen(channel); - print(sprintf("octave_ch_noise_%s.png", channel), "-dpng","-S1000,800"); - set_graphics_state_screen(state_vec); -endfunction - -function octave_c_tx_screen(channel) - clf; hold on; - ber_per_v_snr('oct','datac0',channel,'bo-') - ber_per_v_snr('ctx','datac0',channel,'bx-') - ber_per_v_snr('oct','datac1',channel,'go-') - ber_per_v_snr('ctx','datac1',channel,'gx-') - ber_per_v_snr('oct','datac3',channel,'ro-') - ber_per_v_snr('ctx','datac3',channel,'rx-') - xlabel('SNR (dB)'); grid; - hold off; - if strcmp(channel,"awgn") - axis([-6 8 1E-3 1]); - else - axis([-2 12 1E-3 1]); - end - title(sprintf('Comparsion of Octave Tx and C Tx (no compression) (%s)', channel)); -endfunction - -function octave_c_tx_print(channel) - state_vec = set_graphics_state_print(); - octave_c_tx_screen(channel); - print(sprintf("octave_c_tx_%s.png", channel), "-dpng","-S1000,800"); - set_graphics_state_screen(state_vec); -endfunction - -function octave_c_tx_comp_screen(channel) - clf; hold on; - ber_per_v_snr('oct','datac0',channel,'bo-') - ber_per_v_snr('ctxc','datac0',channel,'bx-') - ber_per_v_snr('oct','datac1',channel,'go-') - ber_per_v_snr('ctxc','datac1',channel,'gx-') - ber_per_v_snr('oct','datac3',channel,'ro-') - ber_per_v_snr('ctxc','datac3',channel,'rx-') - xlabel('SNR (dB)'); grid; - hold off; - if strcmp(channel,"awgn") - axis([-6 8 1E-3 1]); - else - axis([-2 12 1E-3 1]); - end - title(sprintf('Comparsion of Octave Tx and C Tx (with compression) (%s)', channel)); -endfunction - -function octave_c_tx_comp_print(channel) - state_vec = set_graphics_state_print(); - octave_c_tx_comp_screen(channel); - print(sprintf("octave_c_tx_comp_%s.png", channel), "-dpng","-S1000,800"); - set_graphics_state_screen(state_vec); -endfunction - -% composite AWGN and MPP for compressed -function c_tx_comp_screen - clf; hold on; - per_v_snr('ctxc','datac0','awgn','bo-') - per_v_snr('ctxc','datac1','awgn','go-') - per_v_snr('ctxc','datac3','awgn','ro-') - per_v_snr('ctxc','datac4','awgn','co-') - per_v_snr('ctxc','datac13','awgn','mo-') - per_v_snr('ctxc','datac0','mpp','bx-') - per_v_snr('ctxc','datac1','mpp','gx-') - per_v_snr('ctxc','datac3','mpp','rx-') - per_v_snr('ctxc','datac4','mpp','cx-') - per_v_snr('ctxc','datac13','mpp','mx-') - xlabel('SNR (dB)'); ylabel('PER'); grid; - hold off; - axis([-10 10 1E-3 1]); - title('PER of C Raw Data Modes (with compression)'); -endfunction - -function c_tx_comp_print; - state_vec = set_graphics_state_print(); - c_tx_comp_screen; - print("c_tx_comp.png", "-dpng","-S1000,800"); - set_graphics_state_screen(state_vec); -endfunction - -function c_tx_comp_thruput_screen - clf; hold on; - thruput_v_snr('ctxc','datac0','awgn','bo-') - thruput_v_snr('ctxc','datac1','awgn','go-') - thruput_v_snr('ctxc','datac3','awgn','ro-') - thruput_v_snr('ctxc','datac4','awgn','co-') - thruput_v_snr('ctxc','datac13','awgn','mo-') - thruput_v_snr('ctxc','datac0','mpp','bx-') - thruput_v_snr('ctxc','datac1','mpp','gx-') - thruput_v_snr('ctxc','datac3','mpp','rx-') - thruput_v_snr('ctxc','datac4','mpp','cx-') - thruput_v_snr('ctxc','datac13','mpp','mx-') - xlabel('SNR (dB)'); ylabel('bits/s'); grid; - hold off; - axis([-10 10 0 1000]); - title(' Throughput for C Tx (with compression)'); - legend('location','west'); -endfunction - -function c_tx_comp_thruput_print; - state_vec = set_graphics_state_print; - c_tx_comp_thruput_screen; - print("c_tx_comp_thruput.png", "-dpng","-S1000,800"); - set_graphics_state_screen(state_vec); -endfunction - -#{ -figure(1); octave_ch_noise_screen; -figure(2); octave_c_tx_screen; -figure(3); octave_c_tx_comp_screen -figure(4); snrest_snr_screen; - -figure(5); octave_ch_noise_print; -figure(6); octave_c_tx_print; -figure(7); octave_c_tx_comp_print; -figure(8); snrest_snr_print; -#} diff --git a/octave/tdetphase.m b/octave/tdetphase.m deleted file mode 100644 index 8d7f96c..0000000 --- a/octave/tdetphase.m +++ /dev/null @@ -1,84 +0,0 @@ -% tdetphase.m -% David Rowe August 2017 -% -% Testing Hilbert Transform recover of phase from magnitude spectra - -newamp; -Fs = 8000; - -w = 2*pi*500/Fs; gamma = 0.95 -ak = [1 -2*gamma*cos(w) gamma*gamma]; -Nfft = 512; - -% Test 1 - compare phase from freqz for 2nd order system (all pole filter) -% - uses internal test of determine_phase() - -h = freqz(1,ak,Nfft/2); - -% note dummy_model not used, as determine_phase() is used in test mode - -L = 20; Wo = pi/(L+1); -dummy_model = [Wo L ones(1,L)]; -phase = determine_phase(dummy_model, 1, Nfft, ak); - -fg = 1; -figure(fg++); clf; -subplot(211); plot(20*log10(abs(h))); title('test 1'); -subplot(212); plot(angle(h)); hold on; plot(phase(1:Nfft/2),'g'); hold off; - -% Test 2 - feed in harmonic magnitudes - -F0 = 100; Wo = 2*pi*F0/Fs; L = floor(pi/Wo); -Am = zeros(1,L); -for m=1:L - b = round(m*Wo*Nfft/(2*pi)); - Am(m) = abs(h(b)); -end -AmdB = 20*log10(Am); -model = [Wo L Am]; -[phase Gdbfk s] = determine_phase(model, 1, Nfft); - -fftx = (1:Nfft/2)*(Fs/Nfft); -harmx = (1:L)*Wo*Fs/(2*pi); - -figure(fg++); clf; -subplot(211); plot(fftx, Gdbfk(1:Nfft/2)); -subplot(212); plot(s(1:Nfft/2)) - -figure(fg++); clf; -subplot(211); plot(fftx, 20*log10(abs(h))); - hold on; plot(harmx, AmdB, 'g+'); plot(fftx, Gdbfk(1:Nfft/2), 'r'); hold off; -subplot(212); plot(fftx, angle(h)); hold on; plot(fftx, phase(1:Nfft/2),'g'); hold off; - -% Test 3 - Use real harmonic amplitudes - -model = load("../build_linux/src/hts1a_model.txt"); -phase_orig = load("../build_linux/src/hts1a_phase.txt"); - -f = 184; -Wo = model(f,1); L = model(f,2); Am = model(f,3:L+2); AmdB = 20*log10(Am); -[phase Gdbfk s] = determine_phase(model, f, Nfft); - -fftx = (1:Nfft/2)*(Fs/Nfft); -harmx = (1:L)*Wo*Fs/(2*pi); - -figure(fg++); clf; -subplot(211); plot(fftx, Gdbfk(1:Nfft/2)); -subplot(212); plot(s(1:Nfft/2)) - -figure(fg++); clf; -subplot(211); plot(harmx, AmdB, 'g+'); - hold on; plot(fftx, Gdbfk(1:Nfft/2), 'r'); hold off; -subplot(212); plot(fftx, phase(1:Nfft/2),'g'); - -% synthesise using phases - -N = 320; -s = s_phase = zeros(1,N); -for m=1:L/4 - s = s + Am(m)*cos(m*Wo*(1:N) + phase_orig(f,m)); - b = round(m*Wo*Nfft/(2*pi)); - s_phase = s_phase + Am(m)*cos(m*Wo*(1:N) + phase(b)); -end -figure(fg++); clf; -subplot(211); plot(s); subplot(212); plot(s_phase,'g'); diff --git a/octave/test_dqpsk2.m b/octave/test_dqpsk2.m deleted file mode 100644 index d421d62..0000000 --- a/octave/test_dqpsk2.m +++ /dev/null @@ -1,465 +0,0 @@ -% test_dqpsk2.m -% David Rowe April 2014 -% -% DQPSK modem simulation inclduing filtering to test modulating modem -% tx power based on speech energy. Unlike test_dpsk runs at sample -% rate Fs. - -1; - -% main test function - -function sim_out = ber_test(sim_in) - Fs = 8000; - - verbose = sim_in.verbose; - framesize = sim_in.framesize; - Ntrials = sim_in.Ntrials; - Esvec = sim_in.Esvec; - phase_offset = sim_in.phase_offset; - w_offset = sim_in.w_offset; - plot_scatter = sim_in.plot_scatter; - Rs = sim_in.Rs; - hf_sim = sim_in.hf_sim; - Nhfdelay = floor(sim_in.hf_delay_ms*Fs/1000); - Nc = sim_in.Nc; - symbol_amp = sim_in.symbol_amp; - - bps = 2; - Nsymb = framesize/bps; - for k=1:Nc - prev_sym_tx(k) = qpsk_mod([0 0]); - prev_sym_rx(k) = qpsk_mod([0 0]); - end - - % design root nyquist (root raised cosine) filter and init tx and rx filter states - - alpha = 0.5; T=1/Fs; Nfiltsym=7; M=Fs/Rs; - if floor(Fs/Rs) != Fs/Rs - printf("oversampling ratio must be an integer\n"); - return; - end - hrn = gen_rn_coeffs(alpha, T, Rs, Nfiltsym, M); - Nfilter = length(hrn); - - % convert "spreading" samples from 1kHz carrier at Fs to complex - % baseband, generated by passing a 1kHz sine wave through PathSim - % with the ccir-poor model, enabling one path at a time. - - Fc = 1000; - fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb"); - spread1k = fread(fspread, "int16")/10000; - fclose(fspread); - fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb"); - spread1k_2ms = fread(fspread, "int16")/10000; - fclose(fspread); - - % down convert to complex baseband - spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))'); - spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))'); - - % remove -2000 Hz image - b = fir1(50, 5/Fs); - spread = filter(b,1,spreadbb); - spread_2ms = filter(b,1,spreadbb_2ms); - - % discard first 1000 samples as these were near 0, probably as - % PathSim states were ramping up. Transpose for convenience - - spread = transpose(spread(1000:length(spread))); - spread_2ms = transpose(spread_2ms(1000:length(spread_2ms))); - - % Determine "gain" of HF channel model, so we can normalise - % carrier power during HF channel sim to calibrate SNR. I imagine - % different implementations of ccir-poor would do this in - % different ways, leading to different BER results. Oh Well! - - hf_gain = 1.0/sqrt(var(spread)+var(spread_2ms)); - - % Start Simulation ---------------------------------------------------------------- - - for ne = 1:length(Esvec) - EsNodB = Esvec(ne); - EsNo = 10^(EsNodB/10); - - variance = Fs/(Rs*EsNo); - if verbose > 1 - printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance); - end - - Terrs = 0; Tbits = 0; - - tx_symb_log = []; - rx_symb_log = []; - noise_log = []; - sim_out.errors_log = []; - sim_out.tx_baseband_log = []; - sim_out.rx_filt_log = []; - symbol_amp_index = 1; - - % init filter memories and LOs - - tx_filter_memory = zeros(Nc, Nfilter); - rx_filter_memory = zeros(Nc, Nfilter); - s_delay_line_filt = zeros(Nc, Nfiltsym); - phase_tx = ones(1,Nc); - phase_rx = ones(1,Nc); - Fcentre = 1500; Fsep = (1+alpha)*Rs; - freq = Fcentre + Fsep*((-Nc/2+0.5):(Nc/2-0.5)); - freq = exp(j*freq*2*pi/Fs); - - % init HF channel - - sc = 1; hf_n = 1; - hf_sim_delay_line = zeros(1,M+Nhfdelay); - freq_sample_hz = Fcentre + ((Fsep*(-Nc/2)):50:(Fsep*(Nc/2))); - freq_sample_rads = (2*pi/Fs)*freq_sample_hz; - hf_model = ones(Ntrials*Nsymb/Nc, length(freq_sample_rads)); % defaults for plotting surface - - % bunch of outputs we log for graphing - - sim_out.errors_log = []; - sim_out.Nerrs = []; - sim_out.snr_log = []; - sim_out.hf_model_pwr = []; - sim_out.tx_fdm_log = []; - C_log = []; - - for nn = 1: Ntrials - - tx_bits = round( rand( 1, framesize ) ); - - % modulate -------------------------------------------- - - s = zeros(1, Nsymb); - for i=1:Nc:Nsymb - for k=1:Nc - tx_symb = qpsk_mod(tx_bits(2*(i-1+k-1)+1:2*(i+k-1))); - s_qpsk(i+k-1) = tx_symb; - tx_symb *= prev_sym_tx(k); - prev_sym_tx(k) = tx_symb; - s(i+k-1) = symbol_amp(symbol_amp_index)*tx_symb; - end - end - symbol_amp_index++; - s_ch = s; - - % Now we start processing frame Nc symbols at a time to model parallel carriers - - tx_fdm_sym_log = []; - for i=1:Nc:Nsymb - - % Delay tx symbols to match delay due to filters. qpsk - % (rather than dqpsk) symbols used for convenience as - % it's easy to shift symbols than pairs of bits - - s_delay_line_filt(:,1:Nfiltsym-1) = s_delay_line_filt(:,2:Nfiltsym); - s_delay_line_filt(:,Nfiltsym) = s_qpsk(i:i+Nc-1); - s_qpsk(i:i+Nc-1) = s_delay_line_filt(:,1); - for k=1:Nc - tx_bits(2*(i-1+k-1)+1:2*(i+k-1)) = qpsk_demod(s_qpsk(i+k-1)); - end - - % tx filter - - tx_baseband = zeros(Nc,M); - - % tx filter each symbol, generate M filtered output samples for each symbol. - % Efficient polyphase filter techniques used as tx_filter_memory is sparse - - tx_filter_memory(:,Nfilter) = s(i:i+Nc-1); - - for k=1:M - tx_baseband(:,k) = M*tx_filter_memory(:,M:M:Nfilter) * hrn(M-k+1:M:Nfilter)'; - end - tx_filter_memory(:,1:Nfilter-M) = tx_filter_memory(:,M+1:Nfilter); - tx_filter_memory(:,Nfilter-M+1:Nfilter) = zeros(Nc,M); - - sim_out.tx_baseband_log = [sim_out.tx_baseband_log tx_baseband]; - - % upconvert - - tx_fdm = zeros(1,M); - - for c=1:Nc - for k=1:M - phase_tx(c) = phase_tx(c) * freq(c); - tx_fdm(k) = tx_fdm(k) + tx_baseband(c,k)*phase_tx(c); - end - end - - sim_out.tx_fdm_log = [sim_out.tx_fdm_log tx_fdm]; - - % HF channel - - if hf_sim - hf_sim_delay_line(1:Nhfdelay) = hf_sim_delay_line(M+1:M+Nhfdelay); - hf_sim_delay_line(Nhfdelay+1:M+Nhfdelay) = tx_fdm; - - tx_fdm = tx_fdm.*spread(sc:sc+M-1) + hf_sim_delay_line(1:M).*spread_2ms(sc:sc+M-1); - tx_fdm *= hf_gain; - - % sample HF channel spectrum in middle of this symbol for plotting - - hf_model(hf_n,:) = hf_gain*(spread(sc+M/2) + exp(-j*freq_sample_rads*Nhfdelay)*spread_2ms(sc+M/2)); - - sc += M; - hf_n++; - end - - tx_fdm_sym_log = [tx_fdm_sym_log tx_fdm ]; - - % AWGN noise and phase/freq offset channel simulation - % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im - - noise = sqrt(variance*0.5)*(randn(1,M) + j*randn(1,M)); - noise_log = [noise_log noise]; - - % apply frequency and phase offset and noise - - for k=1:M - rx_fdm(k) = tx_fdm(k)*exp(j*phase_offset) + noise(k); - phase_offset += w_offset; - end - - % downconvert - - rx_baseband = zeros(Nc,M); - for c=1:Nc - for k=1:M - phase_rx(c) = phase_rx(c) * freq(c); - rx_baseband(c,k) = rx_fdm(k)*phase_rx(c)'; - end - end - - % rx filter - - rx_filter_memory(:,Nfilter-M+1:Nfilter) = rx_baseband; - rx_filt = rx_filter_memory * hrn'; - rx_filter_memory(:,1:Nfilter-M) = rx_filter_memory(:,1+M:Nfilter); - sim_out.rx_filt_log = [sim_out.rx_filt_log rx_filt]; - - s_ch(i:i+Nc-1) = rx_filt; - end - - % est HF model power for entire code frame (which could be several symbols) - - if hf_sim - frame_hf_model = reshape(hf_model(hf_n-Nsymb/Nc:hf_n-1,:),1,(Nsymb/Nc)*length(freq_sample_hz)); - sim_out.hf_model_pwr = [sim_out.hf_model_pwr mean(abs(frame_hf_model).^2)]; - else - sim_out.hf_model_pwr = [sim_out.hf_model_pwr 1]; - end - - % "genie" SNR estimate - - snr = (tx_fdm_sym_log*tx_fdm_sym_log')/(M*variance); - sim_out.snr_log = [sim_out.snr_log snr]; - - % de-modulate - - rx_bits = zeros(1, framesize); - for i=1:Nc:Nsymb - for k=1:Nc - rx_symb = s_ch(i+k-1); - tmp = rx_symb; - rx_symb *= conj(prev_sym_rx(k)/abs(prev_sym_rx(k))); - prev_sym_rx(k) = tmp; - rx_bits((2*(i-1+k-1)+1):(2*(i+k-1))) = qpsk_demod(rx_symb); - rx_symb_log = [rx_symb_log rx_symb]; - end - end - - % ignore data until we have enough frames to fill filter memory - % then count errors - - if nn > ceil(Nfiltsym/(Nsymb/Nc)) - error_positions = xor(rx_bits, tx_bits); - sim_out.errors_log = [sim_out.errors_log error_positions]; - Nerrs = sum(error_positions); - sim_out.Nerrs = [sim_out.Nerrs Nerrs]; - Terrs += Nerrs; - Tbits += length(tx_bits); - end - - end - - TERvec(ne) = Terrs; - BERvec(ne) = Terrs/Tbits; - - if verbose - printf("EsNo (dB): %f Terrs: %d BER %f ", EsNodB, Terrs, Terrs/Tbits); - printf("\n"); - end - if verbose > 1 - printf("Terrs: %d BER %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs, - Terrs/Tbits, var(sim_out.tx_fdm_log), var(noise_log), - var(sim_out.tx_fdm_log)/(Nc*Rs), var(noise_log)/Fs, (var(sim_out.tx_fdm_log)/(Nc*Rs))/(var(noise_log)/Fs)); - end - end - - Ebvec = Esvec - 10*log10(bps); - - sim_out.BERvec = BERvec; - sim_out.Ebvec = Ebvec; - sim_out.TERvec = TERvec; - - if plot_scatter - figure(2); - clf; - scat = rx_symb_log(Nfiltsym*Nc:length(rx_symb_log)) .* exp(j*pi/4); - plot(real(scat), imag(scat),'+'); - title('Scatter plot'); - - figure(3); - clf; - y = 1:Rs*2; - EsNodBSurface = 20*log10(abs(hf_model(y,:))) + EsNodB; - mesh(1:length(freq_sample_hz),y,EsNodBSurface); - grid - title('HF Channel Es/No'); - end - -endfunction - -% Gray coded QPSK modulation function - -function symbol = qpsk_mod(two_bits) - two_bits_decimal = sum(two_bits .* [2 1]); - switch(two_bits_decimal) - case (0) symbol = 1; - case (1) symbol = j; - case (2) symbol = -j; - case (3) symbol = -1; - endswitch -endfunction - -% Gray coded QPSK demodulation function - -function two_bits = qpsk_demod(symbol) - if isscalar(symbol) == 0 - printf("only works with scalars\n"); - return; - end - bit0 = real(symbol*exp(j*pi/4)) < 0; - bit1 = imag(symbol*exp(j*pi/4)) < 0; - two_bits = [bit1 bit0]; -endfunction - -function sim_in = standard_init - sim_in.verbose = 1; - sim_in.plot_scatter = 0; - - sim_in.Esvec = 5:15; - sim_in.Ntrials = 100; - sim_in.framesize = 64; - sim_in.Rs = 100; - sim_in.Nc = 8; - - sim_in.phase_offset = 0; - sim_in.w_offset = 0; - sim_in.phase_noise_amp = 0; - - sim_in.hf_delay_ms = 2; - sim_in.hf_sim = 0; - sim_in.hf_phase_only = 0; - sim_in.hf_mag_only = 0; -endfunction - -function awgn_hf_ber_curves() - sim_in = standard_init(); - - Ebvec = sim_in.Esvec - 10*log10(2); - BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); - - dpsk_awgn = ber_test(sim_in); - sim_in.hf_sim = 1; - dpsk_hf = ber_test(sim_in); - - figure(1); - clf; - semilogy(Ebvec, BER_theory,'r;QPSK theory;') - hold on; - semilogy(dpsk_awgn.Ebvec, dpsk_awgn.BERvec,'g;DQPSK;') - semilogy(dpsk_hf.Ebvec, dpsk_hf.BERvec,'g;DQPSK HF;') - hold off; - xlabel('Eb/N0') - ylabel('BER') - grid("minor") - axis([min(Ebvec) max(Ebvec) 1E-3 1]) -end - -sim_in = standard_init(); - -% energy file sampled every 10ms - -load ../src/ve9qrp.txt -pdB=10*log10(ve9qrp); -for i=1:length(pdB) - if pdB(i) < 0 - pdB(i) = 0; - end -end - -% Down sample to 40ms rate used for 1300 bit/s codec, every 4th sample is transmitted - -pdB = pdB(4:4:length(pdB)); - -% Use linear mapping function in dB domain to map to symbol power - -%power_map_x = [ 0 20 24 40 50 ]; -%power_map_y = [--6 -6 0 6 6]; -power_map_x = [ 0 50 ]; -power_map_y = [ -15 12]; -mapped_pdB = interp1(power_map_x, power_map_y, pdB); - -sim_in.symbol_amp = 10 .^ (mapped_pdB/20); -%sim_in.symbol_amp = ones(1,length(pdB)); -sim_in.plot_scatter = 1; -sim_in.verbose = 2; -sim_in.hf_delay_ms = 2; -sim_in.hf_sim = 1; -sim_in.Esvec = 10; -sim_in.Ntrials = 400; - -dqpsk_pwr_hf = ber_test(sim_in); - -% note: need way to test that power is aligned with speech - -figure(4) -clf; -plot((1:sim_in.Ntrials)*80*4, pdB(1:sim_in.Ntrials)); -hold on; -plot((1:sim_in.Ntrials)*80*4, mapped_pdB(1:sim_in.Ntrials),'r'); -hold off; - -figure(5) -clf; -s = load_raw("../raw/ve9qrp.raw"); -M=320; M_on_2 = M/2; % processing delay between input speech and centre of analysis window -subplot(211) -plot(M_on_2:(M_on_2-1+sim_in.Ntrials*M),s(1:sim_in.Ntrials*M)) -hold on; -plot((1:sim_in.Ntrials)*M, 5000*sim_in.symbol_amp(1:sim_in.Ntrials),'r'); -hold off; -axis([1 sim_in.Ntrials*M -3E4 3E4]); -subplot(212) -plot(real(dqpsk_pwr_hf.tx_fdm_log)); - - -figure(6) -clf; -plot((1:sim_in.Ntrials)*M, 20*log10(sim_in.symbol_amp(1:sim_in.Ntrials)),'b;Es (dB);'); -hold on; -plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.hf_model_pwr),'g;Fading (dB);'); -plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.snr_log),'r;Es/No (dB);'); - -ber = dqpsk_pwr_hf.Nerrs/sim_in.framesize; -ber_clip = ber; -ber_clip(find(ber > 0.2)) = 0.2; -plot((1:length(ber_clip))*M, -20+100*ber_clip,'k;BER (0-20%);'); -hold off; -axis([1 sim_in.Ntrials*M -20 20]) - -fep=fopen("dqpsk_errors_pwr.bin","wb"); fwrite(fep, dqpsk_pwr_hf.errors_log, "short"); fclose(fep); -fber=fopen("ber.bin","wb"); fwrite(fber, ber, "float"); fclose(fber); diff --git a/octave/trellis.m b/octave/trellis.m deleted file mode 100644 index 01649ec..0000000 --- a/octave/trellis.m +++ /dev/null @@ -1,594 +0,0 @@ -% trellis.m -% David Rowe July 2021 -% -% Testing trellis decoding of Codec 2 Vector Quantiser (VQ) -% information. Uses soft decision information, probablility of state -% transitions, and left over redundancy to correct errors on VQ -% reception. -% -% VQ indexes are transmitted as codewords mapped to +-1 -% -% y = c + n -% -% where c is the transmitted codeword, y is the received codeword, -% and n is Gaussian noise. -% -% This script generates the test data files: -% -% cd codec2/build_linux -% ../script/train_trellis.sh -% -% Results so far (August 2021): -% -% 1/ 2dB improvement with nstages=3, dec=1 -% 2/ No useful improvement with nstages=3, dec=4. This is required for a practical codec to -% get a useful bit rate. - -1; - -% converts a decimal value to a soft dec binary value -function c = dec2sd(dec, nbits) - - % convert to binary - - c = zeros(1,nbits); - for j=0:nbits-1 - mask = 2.^j; - if bitand(dec,mask) - c(nbits-j) = 1; - end - end - - % map to +/- 1 - - c = -1 + 2*c; -endfunction - - -% y is vector of received soft decision values (e.g +/-1 + noise) -function [txp indexes] = ln_tx_codeword_prob_given_rx_codeword_y(y, nstates, C) - nbits = length(y); - np = 2.^nbits; - - % Find log probability of all possible transmitted codewords - txp = C * y'; - - % return most probable codewords (number of states to search) - [txp indexes] = sort(txp,"descend"); - txp = txp(1:nstates); - indexes = indexes(1:nstates) - 1; -endfunction - -% A matrix of all possible tx codewords C, one per row -function C = precompute_C(nbits) - np = 2.^nbits; - - C = zeros(np, nbits); - for r=0:np-1 - C(r+1,:) = dec2sd(r,nbits); - end - -endfunction - - -% work out transition probability matrix, given lists of current and next -% candidate codewords - -function tp = calculate_tp(vq, sd_table, h_table, indexes_current, indexes_next, verbose) - ntxcw = length(indexes_current); - tp = zeros(ntxcw, ntxcw); - for txcw_current=1:ntxcw - index_current = indexes_current(txcw_current); - for txcw_next=1:ntxcw - index_next = indexes_next(txcw_next); - dist = vq(index_current+1,:) - vq(index_next+1,:); - sd = mean(dist.^2); - p = prob_from_hist(sd_table, h_table, sd); - if bitand(verbose, 0x2) - printf("index_current: %d index_next: %d sd: %f p: %f\n", index_current, index_next, sd, p); - end - tp(txcw_current, txcw_next) = log(p); - end - end -endfunction - - -% y is the sequence received soft decision codewords, each row is one -% codeword in time. sd_table and h_table map SD to -% probability. Returns the most likely transmitted VQ index ind in the -% middle of the codeword sequence y. We search the most likely ntxcw -% tx codewords out of 2^nbits possibilities. - -function ind = find_most_likely_index(y, vq, C, sd_table, h_table, nstages, ntxcw, verbose) - [ncodewords nbits] = size(y); - - % populate the nodes of the trellis with the most likely transmitted codewords - txp = zeros(nstages, ntxcw); indexes = zeros(nstages, ntxcw); - for s=1:nstages - [atxp aindexes] = ln_tx_codeword_prob_given_rx_codeword_y(y(s,:), ntxcw, C); - txp(s,:) = atxp; - indexes(s,:) = aindexes; - end - - if verbose - printf("rx_codewords:\n"); - for r=1:ncodewords - for c=1:nbits - printf("%7.2f", y(r,c)); - end - printf("\n"); - end - - printf("\nProbability of each tx codeword index/binary/ln(prob):\n"); - printf(" "); - for s=1:nstages - printf("Time n%+d ", s - (floor(nstages/2)+1)); - end - printf("\n"); - - for i=1:ntxcw - printf("%d ", i); - for s=1:nstages - ind = indexes(s,i); - printf("%4d %12s %5.2f ", ind, dec2bin(ind,nbits), txp(s, i)); - end - printf("\n"); - end - printf("\n"); - end - - % Determine transition probability matrix for each stage, this - % changes between stages as lists of candidate tx codewords - % changes - - tp = zeros(nstages, ntxcw, ntxcw); - for s=1:nstages-1 - if verbose printf("Calc tp(%d,:,:)\n", s), end - tp(s,:,:) = calculate_tp(vq, sd_table, h_table, indexes(s,:), indexes(s+1,:), verbose); - end - - if verbose - printf("Evaluation of all possible paths:\n"); - printf(" "); - for s=1:nstages - printf(" n%+d", s - (floor(nstages/2)+1)); - end - printf(" indexes"); - printf(" "); - - for s=1:nstages - printf(" txp(%d)", s-1); - if s < nstages - printf(" tp(%d,%d) ", s-1,s); - end - end - printf(" prob max_prob\n"); - end - - % OK lets search all possible paths and find most probable - - n = ones(1,nstages); % current node at each stage through trellis, describes current path - max_prob = -100; - do - - if bitand(verbose, 0x4) - printf(" "); - for s=1:nstages - printf("%4d", n(s)-1); - end - printf(" "); - for s=1:nstages - printf("%4d ", indexes(s,n(s))); - end - end - - % find the probability of current path - prob = 0; - for s=1:nstages - prob += txp(s, n(s)); - if bitand(verbose,0x4) - printf("%8.2f ", txp(s, n(s))); - end - if s < nstages - prob += tp(s, n(s), n(s+1)); - if bitand(verbose,0x4) - printf("%8.2f ", tp(s, n(s), n(s+1))); - end - end - end - - if (prob > max_prob) - max_prob = prob; - max_n = n; - end - - if bitand(verbose,0x4) - printf("%9.2f %9.2f\n", prob, max_prob); - end - - % next path - - s = nstages; - n(s)++; - while (s && (n(s) == (ntxcw+1))) - n(s) = 1; - s--; - if s > 0 - n(s)++; - end - end - until (sum(n) == nstages) - - middle = floor(nstages/2)+1; - ind = indexes(middle, max_n(middle)); - if verbose - printf("\nMost likely path through nodes... "); - for s=1:nstages - printf("%4d ", max_n(s)-1); - end - printf("\nMost likely path through indexes: "); - for s=1:nstages - printf("%4d ", indexes(s,max_n(s))); - end - printf("\nMost likely VQ index at time n..: %4d\n", ind); - end -endfunction - - -% Given a normalised histogram, estimate probability from SD -function p = prob_from_hist(sd_table, h_table, sd) - p = interp1 (sd_table, h_table, sd, "extrap", "nearest"); -endfunction - - -% Calculate a normalised histogram of the SD of adjacent frames from -% a file of output vectors from the VQ. -function [sd_table h_table] = vq_hist(vq_output_fn, dec=1) - K=20; K_st=2+1; K_en=16+1; - vq_out = load_f32(vq_output_fn, K); - [r c]= size(vq_out); - diff = vq_out(dec+1:end,K_st:K_en) - vq_out(1:end-dec,K_st:K_en); - % Octave efficient way to determine MSE or each row of matrix - sd_adj = meansq(diff'); - [h_table sd_table] = hist(sd_adj,100,1); - h_table = max(h_table, 1E-5); -endfunction - - -% vector quantise a sequence of target input vectors, returning the VQ indexes and -% quantised vectors target_ -function [indexes target_] = vector_quantiser(vq, target, verbose=1) - [vq_size K] = size(vq); - [ntarget tmp] = size(target); - target_ = zeros(ntarget,K); - indexes = zeros(1,ntarget); - for i=1:ntarget - best_e = 1E32; - for ind=1:vq_size - e = sum((vq(ind,:)-target(i,:)).^2); - if verbose printf("i: %d ind: %d e: %f\n", i, ind, e), end; - if e < best_e - best_e = e; - best_ind = ind; - end - end - if verbose printf("best_e: %f best_ind: %d\n", best_e, best_ind), end; - target_(i,:) = vq(best_ind,:); indexes(i) = best_ind; - end -endfunction - - -% faster version of vector quantiser -function [indexes target_] = vector_quantiser_fast(vq, target, verbose=1) - [vq_size K] = size(vq); - [ntarget tmp] = size(target); - target_ = zeros(ntarget,K); - indexes = zeros(1,ntarget); - - % pre-compute energy of each VQ vector - vqsq = zeros(vq_size,1); - for i=1:vq_size - vqsq(i) = vq(i,:)*vq(i,:)'; - end - - % use efficient matrix multiplies to search for best match to target - for i=1:ntarget - best_e = 1E32; - e = vqsq - 2*(vq * target(i,:)'); - [best_e best_ind] = min(e); - if verbose printf("best_e: %f best_ind: %d\n", best_e, best_ind), end; - target_(i,:) = vq(best_ind,:); indexes(i) = best_ind; - end -endfunction - - -% VQ a target sequence of frames then run a test using vanilla uncoded/trellis decoder -function results = run_test(target, vq, sd_table, h_table, ntxcw, nstages, EbNo, verbose) - [frames tmp] = size(target); - [vq_length tmp] = size(vq); - nbits = log2(vq_length); - nerrors = 0; - nerrors_vanilla = 0; - tbits = 0; - nframes = 0; - nper = 0; - nper_vanilla = 0; - - C = precompute_C(nbits); - - % Vector Quantise target vectors sequence - [tx_indexes target_ ] = vector_quantiser_fast(vq, target, verbose); - % use convention of indexes starting from 0 - tx_indexes -= 1; - % mean SD of VQ with no errors - diff = target - target_; - mse_noerrors = mean(diff(:).^2); - - % construct tx symbol codewords from VQ indexes - tx_codewords = zeros(frames, nbits); - for f=1:frames - tx_codewords(f,:) = dec2sd(tx_indexes(f), nbits); - end - - rx_codewords = tx_codewords + randn(frames, nbits)*sqrt(1/(2*EbNo)); - rx_indexes = zeros(1,frames); - rx_indexes_vanilla = ones(1,frames); - - ns2 = floor(nstages/2); - for f=ns2+1:frames-ns2 - %if f==10 verbose = 1+0x2, else verbose = 0;, end - if verbose - printf("f: %d tx_indexes: ", f); - for i=f-ns2:f+ns2 - printf("%d ", tx_indexes(i)); - end - printf("\n"); - end - tx_bits = tx_codewords(f,:) > 0; - if verbose - printf("tx_bits: "); - for i=1:nbits - printf("%d",tx_bits(i)); - end - printf("\n"); - end - rx_bits_vanilla = rx_codewords(f,:) > 0; - rx_indexes(f) = find_most_likely_index(rx_codewords(f-ns2:f+ns2,:)*EbNo, - vq, C, sd_table, h_table, nstages, ntxcw, verbose); - rx_bits = dec2sd(rx_indexes(f), nbits) > 0; - rx_indexes_vanilla(f) = sum(rx_bits_vanilla .* 2.^(nbits-1:-1:0)); - errors = sum(xor(tx_bits, rx_bits)); - nerrors += errors; - if errors nper++;, end - errors = sum(xor(tx_bits, rx_bits_vanilla)); - nerrors_vanilla += errors; - if errors nper_vanilla++;, end - if verbose - printf("[%d] %d %d\n", f, nerrors, nerrors_vanilla); - end - tbits += nbits; - nframes++; - end - - EbNodB = 10*log10(EbNo); - target = target(ns2+1:frames-ns2,:); - target_vanilla_ = vq(rx_indexes_vanilla(ns2+1:frames-ns2)+1,:); - target_ = vq(rx_indexes(ns2+1:frames-ns2)+1,:); - diff_vanilla = target - target_vanilla_; - mse_vanilla = mean(diff_vanilla(:).^2); - diff = target - target_; - mse = mean(diff(:).^2); - printf("Eb/No: %3.2f dB nframes: %2d nerrors %3d %3d BER: %4.3f %4.3f PER: %3.2f %3.2f mse: %3.2f %3.2f %3.2f\n", - EbNodB, nframes, nerrors, nerrors_vanilla, nerrors/tbits, nerrors_vanilla/tbits, - nper/nframes, nper_vanilla/nframes, - mse_noerrors, mse, mse_vanilla); - results.ber = nerrors/tbits; - results.ber_vanilla = nerrors_vanilla/tbits; - results.per = nper/nframes; - results.per_vanilla = nper_vanilla/nframes; - results.mse_noerrors = mse_noerrors; - results.mse = mse; - results.mse_vanilla = mse_vanilla; - results.tx_indexes = tx_indexes; - results.rx_indexes = rx_indexes; - results.rx_indexes_vanilla = rx_indexes_vanilla; -endfunction - -% Simulations --------------------------------------------------------------------- - -% top level function to set up and run a test -function [results target_] = test_trellis(target_fn, nframes=100, dec=1, ntxcw=8, nstages=3, EbNodB=3, verbose=0) - K = 20; K_st=2+1; K_en=16+1; - vq_fn = "../build_linux/vq_stage1_bs004.f32"; - vq_output_fn = "../build_linux/all_speech_8k_test.f32"; - - % load VQ - vq = load_f32(vq_fn, K); - [vq_size tmp] = size(vq); - vqsub = vq(:,K_st:K_en); - - % load file of VQ-ed vectors to train up SD PDF estimator - [sd_table h_table] = vq_hist(vq_output_fn, dec); - - % load sequence of target vectors we wish to VQ - target = load_f32(target_fn, K); - - % limit test to the first nframes vectors - if nframes != -1 - last = nframes; - else - last = length(target); - end - target = target(1:dec:last,K_st:K_en); - - % run a test - EbNo=10^(EbNodB/10); - results = run_test(target, vqsub, sd_table, h_table, ntxcw, nstages, EbNo, verbose); - if verbose - for f=2:nframes-1 - printf("f: %03d tx_index: %04d rx_index: %04d\n", f, results.tx_indexes(f), results.rx_indexes(f)); - end - end - - % return full band vq-ed vectors - target_ = zeros(last,K); - target_(1:dec:last,:) = vq(results.rx_indexes+1,:); - - % use linear interpolation to restore original frame rate - for f=1:dec:last-dec - prev = f; next = f + dec; - for g=prev+1:next-1 - cnext = (g-prev)/dec; cprev = 1 - cnext; - target_(g,:) = cprev*target_(prev,:) + cnext*target_(next,:); - %printf("f: %d g: %d cprev: %f cnext: %f\n", f, g, cprev, cnext); - end - end -endfunction - -% Plot histograms of SD at different decimations in time -function vq_hist_dec(vq_output_fn) - figure(1); clf; - [sd_table h_table] = vq_hist(vq_output_fn, dec=1); - plot(sd_table, h_table, "b;dec=1;"); - hold on; - [sd_table h_table] = vq_hist(vq_output_fn, dec=2); - plot(sd_table, h_table, "r;dec=2;"); - [sd_table h_table] = vq_hist(vq_output_fn, dec=3); - plot(sd_table, h_table, "g;dec=3;"); - [sd_table h_table] = vq_hist(vq_output_fn, dec=4); - plot(sd_table, h_table, "c;dec=4;"); - hold off; - axis([0 300 0 0.5]) - xlabel("SD dB*dB"); title('Histogram of SD(n,n+1)'); -endfunction - -% Automated tests for vanilla and fast VQ search functions -function test_vq(vq_fn) - K=20; - vq = load_f32(vq_fn, K); - vq_size = 100; - target = vq(1:vq_size,:); - indexes = vector_quantiser(target,target, verbose=0); - assert(indexes == 1:vq_size); - printf("Vanilla OK!\n"); - indexes = vector_quantiser_fast(target,target, verbose=0); - assert(indexes == 1:vq_size); - printf("Fast OK!\n"); -endfunction - -% Test trellis decoding a single vector in a sequence of 3 -function ind = run_test_single(tx_codewords, ntxcw, var, verbose) - nstages = 3; - nbits = 2; - - rx_codewords = tx_codewords + randn(nstages, nbits)*var; - vq = [0 0 0 1; - 0 0 1 0; - 0 1 0 0; - 1 0 0 0]; - sd_table = [0 1 2 4]; - h_table = [0.5 0.25 0.15 0.1]; - C = precompute_C(nbits); - ind = find_most_likely_index(rx_codewords, vq, C, sd_table, h_table, nstages, ntxcw, verbose); -endfunction - -% Series of single point sanity checks -function test_single - printf("Single vector decode tests....\n"); - ind = run_test_single([-1 -1; -1 -1; -1 -1], ntxcw=1, var=0, verbose=0); - assert(ind == 0); - printf("00 with no noise OK!\n"); - - ind = run_test_single([-1 1; 1 1; -1 1], ntxcw=1, var=0, verbose=0); - assert(ind == 3); - printf("11 with no noise OK!\n"); - - ind = run_test_single([-1 -1; -1 1; -1 -1], ntxcw=4, var=1, verbose=0); - assert(ind == 1); - printf("01 with noise OK!\n"); -endfunction - -% BPSK simulation to check noise injection -function test_bpsk_ber - nbits = 12; - frames = 10000; - tx_codewords = zeros(frames,nbits); - tx_bits = zeros(frames,nbits); - for f=1:frames - tx_codewords(f,:) = dec2sd(f, nbits); - tx_bits(f,:) = tx_codewords(f,:) > 0; - end - - EbNodB = 5; - EbNo = 10^(EbNodB/10); - rx_codewords = tx_codewords + randn(frames, nbits)*sqrt(1/(2*EbNo)); - rx_bits = rx_codewords > 0; - nerrors = sum(xor(tx_bits, rx_bits)(:)); - tbits = frames*nbits; - printf("EbNo: %4.2f dB tbits: %d errs: %d BER: %4.3f %4.3f\n", EbNodB, tbits, nerrors, nerrors/tbits, 0.5*erfc(sqrt(EbNo))); -endfunction - -% generate sets of curves -function [EbNodB rms_sd] = run_curves(frames=100, dec=1, nstages=5) - results_log = []; - EbNodB = [0 1 2 3 4 5]; - target_fn = "../build_linux/all_speech_8k_lim.f32"; - - for i=1:length(EbNodB) - results = test_trellis(target_fn, frames, dec, ntxcw=8, nstages, EbNodB(i), verbose=0); - results_log = [results_log results]; - end - for i=1:length(results_log) - ber(i) = results_log(i).ber; - ber_vanilla(i) = results_log(i).ber_vanilla; - per(i) = results_log(i).per; - per_vanilla(i) = results_log(i).per_vanilla; - rms_sd_noerrors(i) = sqrt(results_log(i).mse_noerrors); - rms_sd(i) = sqrt(results_log(i).mse); - rms_sd_vanilla(i) = sqrt(results_log(i).mse_vanilla); - end - - figure(1); clf; semilogy(EbNodB, ber_vanilla, "r+-;uncoded;"); hold on; - semilogy(EbNodB, ber, "g+-;trellis;"); hold off; - grid('minor'); title(sprintf("BER dec=%d nstages=%d",dec,nstages)); - print("-dpng", sprintf("trellis_dec_%d_ber.png",dec)); - - figure(2); clf; semilogy(EbNodB, per_vanilla, "r+-;uncoded;"); hold on; - semilogy(EbNodB, per, "g+-;trellis;"); - grid('minor'); title(sprintf("PER dec=%d nstages=%d",dec,nstages)); - print("-dpng", sprintf("trellis_dec_%d_per.png",dec)); - - figure(3); clf; plot(EbNodB, rms_sd_noerrors, "b+-;no errors;"); hold on; - plot(EbNodB, rms_sd_vanilla, "r+-;uncoded;"); - plot(EbNodB, rms_sd, "g+-;trellis;"); hold off; - grid('minor'); title(sprintf("RMS SD dec=%d nstages=%d",dec,nstages)); - print("-dpng", sprintf("trellis_dec_%d_rms_sd.png",dec)); -endfunction - -function vq_file(vq_fn, dec, EbNodB, in_fn, out_fn) - [results target_] = test_trellis(in_fn, nframes=-1, dec, ntxcw=8, nstages=3, EbNodB, verbose=0); - save_f32(out_fn, target_); -endfunction - -% ------------------------------------------------------------------- - -more off; -randn('state',1); - -% uncomment one of the below to run a test or simulation - -% These two tests show where we are at: -%test_trellis(target_fn, nframes=600, dec=1, ntxcw=8, nstages=3, EbNodB=3, verbose=0); -%test_trellis(target_fn, nframes=600, dec=4, ntxcw=8, nstages=3, EbNodB=3, verbose=0); - -%run_curves(600,1) -%run_curves(600,2) -%run_curves(600,4) -%[EbNodB rms_sd] = run_curves(30*100,3,3) - -%test_trellis(target_fn, nframes=200, dec=1, ntxcw=1, nstages=3, EbNodB=3, verbose=0); -%test_trellis(target_fn, nframes=100, dec=2, ntxcw=8, nstages=3, EbNodB=3, verbose=0); -%test_vq("../build_linux/vq_stage1.f32"); -%vq_hist_dec("../build_linux/all_speech_8k_test.f32"); -%test_single -%test_bpsk_ber diff --git a/octave/vq_700c_eq.m b/octave/vq_700c_eq.m deleted file mode 100644 index 3e57b9a..0000000 --- a/octave/vq_700c_eq.m +++ /dev/null @@ -1,366 +0,0 @@ -% vq_700c.m -% David Rowe May 2019 -% -% Researching Codec 2 700C VQ equaliser ideas -% See also scripts/train_700c_quant.sh, tnewamp1.m - -melvq; -newamp_700c; - -% general purpose plot function for looking at averages of K-band -% sequences in scripts dir and VQs: -% vq_700c_plots({"hts2a.f32" "vk5qi.f32" "train_120_1.txt"}) - -function vq_700c_plots(fn_array) - K = 20; rate_K_sample_freqs_kHz = mel_sample_freqs_kHz(K); - freq_Hz = rate_K_sample_freqs_kHz * 1000; - - figure(1); clf; hold on; axis([200 4000 40 90]); title('Max Hold'); - figure(2); clf; hold on; axis([200 4000 0 40]); title('Average'); - - for i=1:length(fn_array) - [dir name ext] = fileparts(fn_array{i}); - if strcmp(ext, ".f32") - % f32 feature file - fn = sprintf("../build_linux/%s%s", name, ext) - bands = load_f32(fn , K); - else - % text file (e.g. existing VQ) - bands = load(fn_array{i}); - end - % for max hold: break into segments of Nsec, find max, average maximums - % this avoids very rare global peaks setting the max - Nsec = 10; Tframe = 0.01; frames_per_seg = Nsec/Tframe - Nsegs = floor(length(bands)/frames_per_seg) - max_holds = zeros(Nsegs, K); - if Nsegs == 0 - max_holds = max(bands) - else - for s=1:Nsegs - st = (s-1)*frames_per_seg+1; en = st + frames_per_seg - 1; - max_holds(s,:) = max(bands(st:en,:)); - end - max_holds = mean(max_holds); - end - figure(1); plot(freq_Hz, max_holds, '+-', 'linewidth', 2); - figure(2); plot(freq_Hz, mean(bands), '+-', 'linewidth', 2); - end - figure(1); legend(fn_array); grid; xlabel('Freq (Hz)'); ylabel('Amp dB'); - figure(2); legend(fn_array); grid; xlabel('Freq (Hz)'); ylabel('Amp dB'); -endfunction - - -% limit mean of each vector to between lower_lim and upper_lim -function vout = limit_vec(vin, lower_lim, upper_lim) - m = mean(vin'); - vout = zeros(size(vin)); - for i=1:length(vin) - vec_no_mean = vin(i,:) - m(i); - if m(i) < lower_lim - m(i) = lower_lim; - end - if m(i) > upper_lim - m(i) = upper_lim; - end - vout(i,:) = vec_no_mean + m(i); - end -endfunction - - -% single stage vq a target matrix -function errors = vq_targets(vq, targets) - errors = []; - for i=1:length(targets) - [mse_list index_list] = search_vq(vq, targets(i,:), 1); - error = targets(i,:) - vq(index_list(1),:); - errors = [errors; error]; - end -endfunction - - -% single stage vq a target matrix with adaptive EQ, this didn't work - -function [errors eqs] = vq_targets_adap_eq(vq, targets, eqs) - errors = []; gain=0.02; - eq = eqs(end,:); - for i=1:length(targets) - t = targets(i,:) - eq; - mean(t) - %t -= mean(t); - [mse_list index_list] = search_vq(vq, t, 1); - error = t - vq(index_list(1),:); - eq = (1-gain)*eq + gain*error; - errors = [errors; error]; eqs = [eqs; eq]; - end -endfunction - - -% single stage vq a target matrix with block adaptive EQ, this works -% well with nblock == 10 - -function [errors eq] = vq_targets_block_eq(vq, targets, eq, nblock) - errors = []; n = 0; [tmp K] = size(vq); error_eq = zeros(1,K); gain=0.20; - for i=1:length(targets) - t = targets(i,:) - eq; - [mse_list index_list] = search_vq(vq, t, 1); - error = t - vq(index_list(1),:); - error_eq += error; - errors = [errors; error]; - n++; - if n == nblock - eq = 0.99*eq + gain*error_eq/nblock; - n = 0; error_eq = zeros(1,K); - end - end -endfunction - - -% two stage mbest VQ a target matrix - -function [errors targets_] = vq_targets2(vq1, vq2, targets) - vqset(:,:,1)= vq1; vqset(:,:,2)=vq2; m=5; - [errors targets_] = mbest(vqset, targets, m); -endfunction - - -% two stage mbest VQ a target matrix, with adap_eq - -function [errors targets_ eq] = vq_targets2_adap_eq(vq1, vq2, targets, eq) - vqset(:,:,1)= vq1; vqset(:,:,2)=vq2; m=5; gain=0.02; - errors = []; targets_ = []; - for i=1:length(targets) - t = targets(i,:)-eq; - t -= mean(t')'; - [error target_ indexes] = mbest(vqset, t, m); - % use first stage VQ as error driving adaptive EQ - eq_error = t - vq1(indexes(1),:); - eq = (1-gain)*eq + gain*eq_error; - errors = [errors; error]; targets_ = [targets_; target_]; - end -endfunction - - -% Given target and vq matrices, estimate eq via two metrics. First -% metric seems to work best. Both uses first stage VQ error for EQ - -function [eq1 eq2] = est_eq(vq, targets) - [ntargets K] = size(targets); - [nvq K] = size(vq); - - eq1 = zeros(1,K); eq2 = zeros(1,K); - for i=1:length(targets) - [mse_list index_list] = search_vq(vq, targets(i,:), 1); - - % eq metric 1: average of error for best VQ entry - eq1 += targets(i,:) - vq(index_list(1),:); - - % eq metric 2: average of error across all VQ entries - for j=1:nvq - eq2 += targets(i,:) - vq(j,:); - end - end - - eq1 /= ntargets; - eq2 /= (ntargets*nvq); -endfunction - -function [targets e] = load_targets(fn_target_f32) - nb_features = 41; - K = 20; - - % .f32 files are in scripts directory, first K values rate_K_no_mean vectors - [dir name ext] = fileparts(fn_target_f32); - fn = sprintf("../script/%s_feat.f32", name); - feat = load_f32(fn, nb_features); - e = feat(:,1); - targets = feat(:,2:K+1); -endfunction - -% rather simple EQ in front of VQ - -function [eqs ideal] = est_eq_front(targets) - [tmp K] = size(targets); - ideal = [ 8 10 12 14 14*ones(1,K-1-4) -20]; - eq = zeros(1,K); gain = 0.02; - eqs = []; - for i=1:length(targets) - update = targets(i,:) - ideal; - eq = (1-gain)*eq + gain*update; - eq(find(eq < 0)) = 0; - eqs = [eqs; eq]; - end -endfunction - -function table_across_samples - K = 20; - - % VQ is in .txt file in this directory, we have two to choose from. train_120 is the Codec 2 700C VQ, - % train_all_speech was trained up from a different, longer database, as a later exercise - vq_name = "train_120"; - #vq_name = "train_all_speech"; - vq1 = load(sprintf("%s_1.txt", vq_name)); - vq2 = load(sprintf("%s_2.txt", vq_name)); - - printf("----------------------------------------------------------------------------------\n"); - printf("Sample Initial vq1 vq1_eq2 vq1_eq2 vq2 vq2_eq1 vq2_eq2 \n"); - printf("----------------------------------------------------------------------------------\n"); - - fn_targets = { "cq_freedv_8k_lfboost" "cq_freedv_8k_hfcut" "cq_freedv_8k" "hts1a" "hts2a" "cq_ref" "ve9qrp_10s" "vk5qi" "c01_01_8k" "ma01_01"}; - #fn_targets = {"cq_freedv_8k_lfboost"}; - figs=1; - for i=1:length(fn_targets) - - % load target and estimate eq - [targets e] = load_targets(fn_targets{i}); - eq1 = est_eq(vq1, targets); - eq2s = est_eq_front(targets); - % for these simulation uses fixed EQ sample, rather than letting it vary frame by frame - eq2 = eq2s(end,:); - - % first stage VQ ----------------- - - errors1 = vq_targets(vq1, targets); - errors1_eq1 = vq_targets(vq1, targets-eq1); - errors1_eq2 = vq_targets(vq1, targets-eq2); - - % two stage mbest VQ -------------- - - [errors2 targets_] = vq_targets2(vq1, vq2, targets); - [errors2_eq1 targets_eq1_] = vq_targets2(vq1, vq2, targets-eq1); - [errors2_eq2 targets_eq2_] = vq_targets2(vq1, vq2, targets-eq2); - - % save to .f32 files for listening tests - if strcmp(vq_name,"train_120") - save_f32(sprintf("../script/%s_vq2.f32", fn_targets{i}), targets_); - save_f32(sprintf("../script/%s_vq2_eq1.f32", fn_targets{i}), targets_eq1_); - save_f32(sprintf("../script/%s_vq2_eq2.f32", fn_targets{i}), targets_eq2_); - else - save_f32(sprintf("../script/%s_vq2_as.f32", fn_targets{i}), targets_); - save_f32(sprintf("../script/%s_vq2_as_eq.f32", fn_targets{i}), targets_eq_); - end - printf("%-21s %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n", fn_targets{i}, - var(targets(:)), var(errors1(:)), var(errors1_eq1(:)), var(errors1_eq2(:)), - var(errors2(:)), var(errors2_eq1(:)), var(errors2_eq2(:))); - - figure(figs++); clf; - %plot(var(errors2'),'b;vq2;'); hold on; plot(var(errors2_eq1'),'g;vq2_eq1;'); plot(var(errors2_eq2'),'r;vq2_eq2;'); hold off; - plot(mean(targets),'g;mean(targets);'); hold on; plot(mean(vq1),'g;mean(vq1);'); plot(eq2,'r;eq2;'); hold off; - title(fn_targets{i}); axis([1 K -20 30]); - end -endfunction - - -% interactve, menu driven frame by frame plots - -function interactive(fn_vq_txt, fn_target_f32) - K = 20; - vq = load("train_120_1.txt"); - [targets e] = load_targets(fn_target_f32); - eq1 = est_eq(vq, targets); - - [errors1_eq2 eqs2] = vq_targets_adap_eq(vq, targets, zeros(1,K)); - [errors1_eq2 eqs2] = vq_targets_adap_eq(vq, targets, eqs2(end,:)); - eq2 = eqs2(end,:); - - figure(1); clf; - mesh(e+targets) - figure(2); clf; - plot(eq1,'b;eq1;') - hold on; - plot(mean(targets),'c;mean(targets);'); plot(eq2,'g;eq2;'); - hold off; - figure(3); clf; mesh(eqs2); title('eq2 evolving') - - % enter single step loop - f = 20; neq = 0; eq=zeros(1,K); - do - figure(4); clf; - t = targets(f,:) - eq; - [mse_list index_list] = search_vq(vq, t, 1); - error = t - vq(index_list(1),:); - plot(e(f)+t,'b;target;'); - hold on; - plot(e(f)+vq(index_list,:),'g;vq;'); - plot(error,'r;error;'); - plot(eq,'c;eq;'); - plot([1 K],[e(f) e(f)],'--') - hold off; - axis([1 K -20 80]) - % interactive menu - - printf("\r f: %2d eq: %d ind: %3d var: %3.1f menu: n-next b-back e-eq q-quit", f, neq, index_list(1), var(error)); - fflush(stdout); - k = kbhit(); - - if k == 'n' f+=1; end - if k == 'e' - neq++; - end - if neq == 3 neq = 0; end - if neq == 0 eq = zeros(1,K); end - if neq == 1 eq = eq1; end - if neq == 2 eq = eqs2(f,:); end - if k == 'b' f-=1; end - until (k == 'q') - printf("\n"); -endfunction - - -% Experiment to test iterative approach of block update and remove -% mean (ie frame energy), shows some promise at reducing HF energy -% over several iterations while not affecting already good samples - -function experiment_iterate_block(fn_vq_txt, fn_target_f32) - K = 20; - vq = load("train_120_1.txt"); - [targets e] = load_targets(fn_target_f32); - - figure(1); clf; - plot(mean(targets),'b;mean(targets);'); - hold on; - plot(mean(vq), 'g;mean(vq);'); - figure(2); clf; hold on; - eq = zeros(1,K); - for i=1:3 - [errors eq] = vq_targets_block_eq(vq, targets, eq, 10); - figure(1); plot(mean(targets-eq)); - figure(2); plot(eq); - printf("i: %d %6.2f\n", i, var(errors(:))) - end -endfunction - -% Experiment to test EQ of input (before) VQ. We set a threshold on -% when to equalise, so we don't upset already flat-ish samples. This -% is the algorithm used for C at the time of writing (newamp1.c, newamp_700c.m) - -function experiment_front_eq(fn_vq_txt, fn_target_f32) - K = 20; - vq = load("train_120_1.txt"); - [targets e] = load_targets(fn_target_f32); - - [eqs ideal] = est_eq_front(targets); - - figure(1); clf; - plot(mean(targets),'b;mean(targets);'); - hold on; - plot(ideal, 'g;ideal;'); - plot(eqs(end,:), 'r;eq;'); - plot(mean(targets)-eqs(end,:), 'c;equalised;'); - plot(mean(vq),'b--;mean(vq);'); - hold off; - figure(2); clf; mesh(eqs(1:100,:)); title('EQ weights over time'); - ylabel('Time (frames'); xlabel('Freq (mel)'); -endfunction - -more off - -% choose one of these to run first -% You'll need to run scripts/train_700C_quant.sh first to generate the .f32 files - -%interactive("train_120_1.txt", "cq_freedv_8k_lfboost.f32") -%table_across_samples; -%vq_700c_plots({"all_speech_8k.f32" "all_speech_8k_hp300.f32" "dev-clean-8k.f32" "train_8k.f32" } ) -%vq_700c_plots({"ve9qrp_10s.f32" "cq_freedv_8k_lfboost.f32" "cq_ref.f32" "hts1a.f32" "vk5qi.f32"}) -%experiment_iterate_block("train_120_1.txt", "ve9qrp_10s.f32") -%experiment_iterate_block("train_120_1.txt", "cq_freedv_8k_lfboost.f32") -%experiment_front_eq("train_120_1.txt", "cq_freedv_8k_lfboost.f32") diff --git a/octave/vq_compare.m b/octave/vq_compare.m deleted file mode 100644 index 0a156c1..0000000 --- a/octave/vq_compare.m +++ /dev/null @@ -1,349 +0,0 @@ -% vq_compare.m -% David Rowe Sep 2021 -% -% Compare the Eb/No performance of Vector Quantisers (robustness to bit errors) using -% Spectral Distortion (SD) measure. - -#{ - usage: - - 1. Generate the initial VQ (vq_stage1.f32) and input test vector file (all_speech_8k_lim.f32): - - cd codec2/build_linux - ../script/train_trellis.sh - - 2. Run the Psuedo-Gray binary switch tool to optimise the VQ against single bit errors: - - ./misc/vq_binary_switch -d 20 vq_stage1.f32 vq_stage1_bs001.f32 -m 5000 --st 2 --en 16 -f - - This can take a while, but if you ctrl-C at any time it will have saved the most recent optimised VQ. - - 3. Run this script to compare the two VQs: - - octave:34> vq_compare -#} - - -function vq_compare(action="run_curves", vq_fn, dec=1, EbNodB=3, in_fn, out_fn) - more off; - randn('state',1); - graphics_toolkit("gnuplot"); - - if strcmp(action, "run_curves") - run_curves(30*100); - end - if strcmp(action, "vq_file") - vq_file(vq_fn, dec, EbNodB, in_fn, out_fn) - end -endfunction - - -% ------------------------------------------------------------------- - -% converts a decimal value to a soft dec binary value -function c = dec2sd(dec, nbits) - - % convert to binary - - c = zeros(1,nbits); - for j=0:nbits-1 - mask = 2.^j; - if bitand(dec,mask) - c(nbits-j) = 1; - end - end - - % map to +/- 1 - - c = -1 + 2*c; -endfunction - -% fast version of vector quantiser -function [indexes target_] = vector_quantiser_fast(vq, target, verbose=1) - [vq_size K] = size(vq); - [ntarget tmp] = size(target); - target_ = zeros(ntarget,K); - indexes = zeros(1,ntarget); - - % pre-compute energy of each VQ vector - vqsq = zeros(vq_size,1); - for i=1:vq_size - vqsq(i) = vq(i,:)*vq(i,:)'; - end - - % use efficient matrix multiplies to search for best match to target - for i=1:ntarget - best_e = 1E32; - e = vqsq - 2*(vq * target(i,:)'); - [best_e best_ind] = min(e); - if verbose printf("best_e: %f best_ind: %d\n", best_e, best_ind), end; - target_(i,:) = vq(best_ind,:); indexes(i) = best_ind; - end -endfunction - - -% VQ a target sequence of frames then run a test using vanilla uncoded/trellis decoder -function results = run_test(target, vq, EbNo, verbose) - [frames tmp] = size(target); - [vq_length tmp] = size(vq); - nbits = log2(vq_length); - nerrors = 0; - tbits = 0; - nframes = 0; - nper = 0; - - % Vector Quantise target vectors sequence - [tx_indexes target_ ] = vector_quantiser_fast(vq, target, verbose); - % use convention of indexes starting from 0 - tx_indexes -= 1; - % mean SD of VQ with no errors - diff = target - target_; - mse_noerrors = mean(diff(:).^2); - - % construct tx symbol codewords from VQ indexes - tx_codewords = zeros(frames, nbits); - for f=1:frames - tx_codewords(f,:) = dec2sd(tx_indexes(f), nbits); - end - - rx_codewords = tx_codewords + randn(frames, nbits)*sqrt(1/(2*EbNo)); - rx_indexes = zeros(1,frames); - - for f=1:frames - tx_bits = tx_codewords(f,:) > 0; - rx_bits = rx_codewords(f,:) > 0; - rx_indexes(f) = sum(rx_bits .* 2.^(nbits-1:-1:0)); - errors = sum(xor(tx_bits, rx_bits)); - nerrors += errors; - if errors nper++;, end - tbits += nbits; - nframes++; - end - - EbNodB = 10*log10(EbNo); - target_ = vq(rx_indexes+1,:); - diff = target - target_; - mse = mean(diff(:).^2); - printf("Eb/No: %3.2f dB nframes: %3d nerrors: %4d BER: %4.3f PER: %3.2f mse: %3.2f %3.2f\n", - EbNodB, nframes, nerrors, nerrors/tbits, nper/nframes, mse_noerrors, mse); - results.ber = nerrors/tbits; - results.per = nper/nframes; - results.mse_noerrors = mse_noerrors; - results.mse = mse; - results.tx_indexes = tx_indexes; - results.rx_indexes = rx_indexes; -endfunction - -% VQ a target sequence of frames then run a test using a LDPC code -function results = run_test_ldpc(target, vq, EbNo, verbose) - [frames tmp] = size(target); - [vq_length tmp] = size(vq); - nbits = log2(vq_length); - nerrors = 0; - tbits = 0; - nframes = 0; - nper = 0; - - % init LDPC code - mod_order = 4; bps = 2; - modulation = 'QPSK'; - mapping = 'gray'; - max_iterations = 100; demod_type = 0; decoder_type = 0; - ldpc; init_cml(); - tempStruct = load("HRA_56_56.txt"); - b = fieldnames(tempStruct); - ldpcArrayName = b{1,1}; - % extract the array from the struct - HRA = tempStruct.(ldpcArrayName); - [code_param framesize rate] = ldpc_init_user(HRA, modulation, mod_order, mapping); - - % set up noise - EbNodB = 10*log10(EbNo); - EsNodB = EbNodB + 10*log10(rate) + 10*log10(bps); - EsNo = 10^(EsNodB/10); - variance = 1/EsNo; - - % Vector Quantise target vectors sequence - [tx_indexes target_ ] = vector_quantiser_fast(vq, target, verbose); - % use convention of indexes starting from 0 - tx_indexes -= 1; - % mean SD of VQ with no errors - diff = target - target_; - mse_noerrors = mean(diff(:).^2); - - % construct tx frames x nbit matrix using VQ indexes - tx_bits = zeros(frames, nbits); - for f=1:frames - tx_bits(f,:) = dec2sd(tx_indexes(f), nbits) > 0; - end - - % find a superframe size, that has an integer number of nbits and data_bits_per_frame frames - bits_per_superframe = nbits; - while mod(bits_per_superframe,nbits) || mod(bits_per_superframe,code_param.data_bits_per_frame) - bits_per_superframe += nbits; - end - - Nsuperframes = floor(frames*nbits/bits_per_superframe); - Nldpc_codewords = Nsuperframes*bits_per_superframe/code_param.data_bits_per_frame; - frames = Nsuperframes*bits_per_superframe/nbits; - %printf("bits_per_superframe: %d Nldpc_codewords: %d frames: %d\n", bits_per_superframe, Nldpc_codewords, frames); - - % reshape tx_bits matrix into Nldpc_codewords x data_bits_per_frame - tx_bits = tx_bits(1:frames,:); - tx_bits_ldpc = reshape(tx_bits',code_param.data_bits_per_frame, Nldpc_codewords)'; - - % modulate tx symbols - tx_symbols = []; - for nn=1:Nldpc_codewords - [tx_codeword atx_symbols] = ldpc_enc(tx_bits_ldpc(nn,:), code_param); - tx_symbols = [tx_symbols atx_symbols]; - end - - noise = sqrt(variance*0.5)*(randn(1,length(tx_symbols)) + j*randn(1,length(tx_symbols))); - rx_symbols = tx_symbols+noise; - - % LDPC decode - for nn = 1:Nldpc_codewords - st = (nn-1)*code_param.coded_syms_per_frame + 1; - en = (nn)*code_param.coded_syms_per_frame; - - arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo, ones(1,code_param.coded_syms_per_frame)); - rx_bits_ldpc(nn,:) = arx_codeword(1:code_param.data_bits_per_frame); - end - - % reshape rx_bits_ldpc matrix into frames x nbits - rx_bits = reshape(rx_bits_ldpc',nbits,frames)'; - - rx_indexes = tx_indexes; - for f=1:frames - rx_indexes(f) = sum(rx_bits(f,:) .* 2.^(nbits-1:-1:0)); - errors = sum(xor(tx_bits(f,:), rx_bits(f,:))); - nerrors += errors; - if errors nper++;, end - tbits += nbits; - nframes++; - end - - EbNodB = 10*log10(EbNo); - target_ = vq(rx_indexes+1,:); - diff = target - target_; - mse = mean(diff(:).^2); - printf("Eb/No: %3.2f dB nframes: %4d nerrors: %4d BER: %4.3f PER: %3.2f mse: %3.2f %3.2f\n", - EbNodB, nframes, nerrors, nerrors/tbits, nper/nframes, mse_noerrors, mse); - results.ber = nerrors/tbits; - results.per = nper/nframes; - results.mse = mse; - results.tx_indexes = tx_indexes; - results.rx_indexes = rx_indexes; -endfunction - -% Simulations --------------------------------------------------------------------- - -% top level function to set up and run a test with a specific vq -function [results target_] = run_test_vq(vq_fn, target_fn, nframes=100, dec=1, EbNodB=3, ldpc_en=0, verbose=0) - K = 20; K_st=2+1; K_en=16+1; - - % load VQ - vq = load_f32(vq_fn, K); - [vq_size tmp] = size(vq); - vqsub = vq(:,K_st:K_en); - - % load sequence of target vectors we wish to VQ - target = load_f32(target_fn, K); - - % limit test to the first nframes vectors - if nframes != -1 - last = nframes; - else - last = length(target); - end - target = target(1:dec:last, K_st:K_en); - - % run a test - EbNo=10^(EbNodB/10); - if ldpc_en - results = run_test_ldpc(target, vqsub, EbNo, verbose); - else - results = run_test(target, vqsub, EbNo, verbose); - end - if verbose - for f=2:nframes-1 - printf("f: %03d tx_index: %04d rx_index: %04d\n", f, results.tx_indexes(f), results.rx_indexes(f)); - end - end - - % return full band vq-ed vectors - target_ = zeros(last,K); - target_(1:dec:last,:) = vq(results.rx_indexes+1,:); - - % use linear interpolation to restore original frame rate - for f=1:dec:last-dec - prev = f; next = f + dec; - for g=prev+1:next-1 - cnext = (g-prev)/dec; cprev = 1 - cnext; - target_(g,:) = cprev*target_(prev,:) + cnext*target_(next,:); - %printf("f: %d g: %d cprev: %f cnext: %f\n", f, g, cprev, cnext); - end - end -endfunction - -% generate sets of curves -function run_curves(frames=100, dec=1) - target_fn = "../build_linux/all_speech_8k_lim.f32"; - EbNodB = 0:5; - - results1_ldpc_log = []; - for i=1:length(EbNodB) - results = run_test_vq("../build_linux/vq_stage1.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=1, verbose=0); - results1_ldpc_log = [results1_ldpc_log results]; - end - results4_ldpc_log = []; - for i=1:length(EbNodB) - results = run_test_vq("../build_linux/vq_stage1_bs004.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=1, verbose=0); - results4_ldpc_log = [results4_ldpc_log results]; - end - - results1_log = []; - for i=1:length(EbNodB) - results = run_test_vq("../build_linux/vq_stage1.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=0, verbose=0); - results1_log = [results1_log results]; - end - results4_log = []; - for i=1:length(EbNodB) - results = run_test_vq("../build_linux/vq_stage1_bs004.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=0, verbose=0); - results4_log = [results4_log results]; - end - for i=1:length(results1_log) - ber(i) = results1_log(i).ber; - per(i) = results1_log(i).per; - mse_noerrors(i) = sqrt(results1_log(i).mse_noerrors); - mse_vq1(i) = sqrt(results1_log(i).mse); - mse_vq4(i) = sqrt(results4_log(i).mse); - mse_vq1_ldpc(i) = sqrt(results1_ldpc_log(i).mse); - mse_vq4_ldpc(i) = sqrt(results4_ldpc_log(i).mse); - end - - figure(1); clf; - semilogy(EbNodB, ber, 'g+-;ber;','linewidth', 2); hold on; - semilogy(EbNodB, per, 'b+-;per;','linewidth', 2); - grid('minor'); xlabel('Eb/No(dB)'); - hold off; - - figure(2); clf; - plot(EbNodB, mse_noerrors, "b+-;no errors;"); hold on; - plot(EbNodB, mse_vq1, "g+-;vanilla AWGN;"); - plot(EbNodB, mse_vq4, "b+-;binary switch;"); - plot(EbNodB, mse_vq1_ldpc, "r+-;ldpc (112,56);"); - plot(EbNodB, mse_vq4_ldpc, "k+-;binary switch ldpc (112,56);"); - load trellis_dec3_nstage3.txt - plot(EbNodB, rms_sd, "c+-;binary switch trellis dec3;"); - hold off; grid; title("RMS SD (dB)"); xlabel('Eb/No(dB)'); -endfunction - - -function vq_file(vq_fn, dec, EbNodB, in_fn, out_fn) - [results target_] = run_test_vq(vq_fn, in_fn, nframes=-1, dec, EbNodB, verbose=0); - save_f32(out_fn, target_); -endfunction - - |