FreeDV
/
codec2
mirror of https://github.com/drowe67/codec2.git
1
0
Fork 0
Code Issues Releases Wiki Activity

yet more Octave files rm-ed

ms-warning-cleanup-2
drowe67 2023-07-07 14:42:02 +09:30 committed by David Rowe
parent a291cfaf4c
commit fbbd89ddf8
15 changed files with 0 additions and 14440 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
octave/doppler_spread_ut.m
View File

@ -1,51 +0,0 @@
% doppler_spread_ut.m
% David Rowe Jan 2016
%
% Unit test script for doppler_spread
f = 1;
Fs = 8000;
N = Fs*10;
[spread states] = doppler_spread(f, Fs, N);
% use spreading samples to modulate 1000Hz sine wave
% You can listen to this with: sine1k_1Hz.raw
% $ play -t raw -r 8000 -s -2
s = cos(2*pi*(1:N)*1000/Fs);
s = s .* spread;
s = real(s)*5000;
fs = fopen("sine1k_1Hz.raw","wb"); fwrite(fs,s,"short"); fclose(fs);
% Some plots
x = states.x; y = states.y; b = states.b;
H = freqz(b,1,x);
figure(1)
clf
subplot(211)
plot(x,y,';target;')
title('Gaussian Filter Freq Resp Lin');
legend('boxoff');
subplot(212)
plot(x,20*log10(y),';target;')
hold on;
plot(x,20*log10(y),'g+;actual;')
hold off;
axis([0 f*10/2 -60 0])
title('Gaussian Filter Freq Resp dB');
xlabel('Freq (Hz)');
legend('boxoff');
figure(2);
subplot(211)
plot(abs(spread))
title('Spreading Function Magnitude');
subplot(212)
plot(s)
title('1000Hz Sine Wave');
xlabel('Time (samples)')

65
octave/estsnr.m
View File

@ -1,65 +0,0 @@
% estsnr.m
% David Rowe May 2017
%
% estimate SNR of a sinewave in noise
function snr_dB = estsnr(x, Fs=8000, Nbw = 3000)
[nr nc] = size(x);
if nr == 1
x = x';
end
% find peak in +ve side of spectrum, ignoring DC
L = length(x);
X = abs(fft(x));
st = floor(0.05*L); en = floor(0.45*L);
[A mx_ind]= max(X(st:en));
mx_ind += st;
% signal energy might be spread by doppler, so sum energy
% in frequencies +/- 1%
s_st = floor(mx_ind*0.99); s_en = floor(mx_ind*1.01);
S = sum(X(s_st:s_en).^2);
% real signal, so -ve power is the same
S = 2*S;
SdB = 10*log10(S);
printf("Signal Power S: %3.2f dB\n", SdB);
% locate a band of noise next to it and find power in band
st = floor(mx_ind+0.05*(L/2));
en = st + floor(0.1*(L/2));
N = sum(X(st:en).^2);
% scale this to obtain total noise power across total bandwidth
N *= L/(en-st);
NdB = 10*log10(N);
printf("Noise Power N: %3.2f dB\n", NdB);
% scale noise to designed noise bandwidth /2 fudge factor as its a
% real signal, wish I had a better way to explain that!
NodB = NdB - 10*log10(Fs/2);
NscaleddB = NodB + 10*log10(Nbw);
snr_dB = SdB - NscaleddB;
figure(1); clf;
plot(20*log10(X(1:L/2)),'b');
hold on;
plot([s_st s_en], [NdB NdB]- 10*log10(L), 'r');
plot([st en], [NdB NdB]- 10*log10(L), 'r');
hold off;
top = 10*ceil(SdB/10);
bot = NodB - 20;
axis([1 L/2 bot top]);
grid
grid("minor")
endfunction

272
octave/hf_modem_curves.m
View File

@ -1,272 +0,0 @@
% hf_modem_curves
% David Rowe Feb 2017
%
% Ideal implementations of a bunch of different HF modems, used to
% generate plots for a blog post.
#{
[X] ideal AWGN/HF curves
[X] exp AWGN QPSK curves
[X] exp AWGN DQPSK curves
[X] exp HF channel model
[ ] diversity
[ ] COHPSK frames
+ would require multiple carriers
+ filtering or OFDM
#}
1;
% 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)
bit0 = real(symbol*exp(j*pi/4)) < 0;
bit1 = imag(symbol*exp(j*pi/4)) < 0;
two_bits = [bit1 bit0];
endfunction
% Rate Rs modem simulation model -------------------------------------------------------
function sim_out = ber_test(sim_in)
bps = 2; % two bits/symbol for QPSK
Rs = 50; % symbol rate (needed for HF model)
verbose = sim_in.verbose;
EbNovec = sim_in.EbNovec;
hf_en = sim_in.hf_en;
% user can supply number of bits per point to get good results
% at high Eb/No
if length(sim_in.nbits) > 1
nbitsvec = sim_in.nbits;
nbitsvec += 100 - mod(nbitsvec,100); % round up to nearest 100
else
nbitsvec(1:length(EbNovec)) = sim_in.nbits;
end
% init HF model
if hf_en
% some typical values
dopplerSpreadHz = 1.0; path_delay = 1E-3*Rs;
nsymb = max(nbitsvec)/2;
spread1 = doppler_spread(dopplerSpreadHz, Rs, nsymb);
spread2 = doppler_spread(dopplerSpreadHz, Rs, nsymb);
hf_gain = 1.0/sqrt(var(spread1)+var(spread2));
% printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1));
end
for ne = 1:length(EbNovec)
% work out noise power -------------
EbNodB = EbNovec(ne);
EsNodB = EbNodB + 10*log10(bps);
EsNo = 10^(EsNodB/10);
variance = 1/EsNo;
nbits = nbitsvec(ne);
nsymb = nbits/bps;
% modulator ------------------------
tx_bits = rand(1,nbits) > 0.5;
tx_symb = [];
prev_tx_symb = 1;
for s=1:nsymb
atx_symb = qpsk_mod(tx_bits(2*s-1:2*s));
if sim_in.dqpsk
atx_symb *= prev_tx_symb;
prev_tx_symb = atx_symb;
end
tx_symb = [tx_symb atx_symb];
end
% channel ---------------------------
rx_symb = tx_symb;
if hf_en
% simplified rate Rs simulation model that doesn't include
% ISI, just freq filtering. We assume perfect phase estimation
% so it's just amplitude distortion.
hf_model1 = hf_model2 = zeros(1, nsymb);
for s=1:nsymb
hf_model1(s) = hf_gain*(spread1(s) + exp(-j*path_delay)*spread2(s));
hf_model = abs(hf_model1(s));
if sim_in.diversity
% include amplitude information from another frequency in channel model
w1 = 7*2*pi;
hf_model2(s) = hf_gain*(spread1(s) + exp(-j*w1*path_delay)*spread2(s));
hf_model = 0.5*abs(hf_model1(s)) + 0.5*abs(hf_model2(s));
end
rx_symb(s) = rx_symb(s).*hf_model;
end
end
% variance is noise power, which is divided equally between real and
% imag components of noise
noise = sqrt(variance*0.5)*(randn(1,nsymb) + j*randn(1,nsymb));
rx_symb += noise;
% demodulator ------------------------------------------
% demodulate rx symbols to bits
rx_bits = [];
prev_rx_symb = 1;
for s=1:nsymb
arx_symb = rx_symb(s);
if sim_in.dqpsk
tmp = arx_symb;
arx_symb *= prev_rx_symb';
prev_rx_symb = tmp;
end
two_bits = qpsk_demod(arx_symb);
rx_bits = [rx_bits two_bits];
end
% count errors -----------------------------------------
error_pattern = xor(tx_bits, rx_bits);
nerrors = sum(error_pattern);
bervec(ne) = nerrors/nbits;
if verbose
printf("EbNodB: % 3.1f nbits: %5d nerrors: %5d ber: %4.3f\n", EbNodB, nbits, nerrors, bervec(ne));
if verbose == 2
figure(2); clf;
plot(rx_symb*exp(j*pi/4),'+','markersize', 10);
mx = max(abs(rx_symb));
axis([-mx mx -mx mx]);
if sim_in.diversity && sim_in.hf_en
figure(3);
plot(1:nsymb, abs(hf_model1), 1:nsymb, abs(hf_model2), 'linewidth', 2);
end
end
end
end
sim_out.bervec = bervec;
endfunction
% -------------------------------------------------------------
function run_single
sim_in.verbose = 2;
sim_in.nbits = 1000;
sim_in.EbNovec = 4;
sim_in.dqpsk = 0;
sim_in.hf_en = 0;
sim_in.diversity = 0;
sim_qpsk = ber_test(sim_in);
endfunction
function run_curves
max_nbits = 1E5;
sim_in.verbose = 1;
sim_in.EbNovec = 0:10;
sim_in.dqpsk = 0;
sim_in.hf_en = 0;
sim_in.diversity = 0;
% AWGN -----------------------------
ber_awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNovec/10)));
sim_in.nbits = min(max_nbits, floor(500 ./ ber_awgn_theory));
sim_qpsk = ber_test(sim_in);
sim_in.dqpsk = 1;
sim_dqpsk = ber_test(sim_in);
% HF -----------------------------
hf_sim_in = sim_in; hf_sim_in.dqpsk = 0; hf_sim_in.hf_en = 1;
hf_sim_in.EbNovec = 0:16;
EbNoLin = 10.^(hf_sim_in.EbNovec/10);
ber_hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1)));
hf_sim_in.nbits = min(max_nbits, floor(500 ./ ber_hf_theory));
sim_qpsk_hf = ber_test(hf_sim_in);
hf_sim_in.dqpsk = 1;
sim_dqpsk_hf = ber_test(hf_sim_in);
hf_sim_in.dqpsk = 0;
hf_sim_in.diversity = 1;
sim_qpsk_hf_div = ber_test(hf_sim_in);
% Plot results --------------------
close all;
figure (1, 'position', [100, 10, 600, 400]); clf;
semilogy(sim_in.EbNovec, ber_awgn_theory,'r+-;QPSK AWGN theory;', 'linewidth', 2)
xlabel('Eb/No (dB)')
ylabel('BER')
grid("minor")
axis([min(sim_in.EbNovec) max(sim_in.EbNovec) 1E-3 1])
hold on;
semilogy([0 4 4], [ber_awgn_theory(5) ber_awgn_theory(5) 1E-3],'k--', 'linewidth', 2);
hold off;
figure (2, 'position', [300, 10, 600, 400]); clf;
semilogy(sim_in.EbNovec, ber_awgn_theory,'r+-;QPSK AWGN theory;','markersize', 10, 'linewidth', 2)
hold on;
semilogy(sim_in.EbNovec, sim_qpsk.bervec,'g+-;QPSK AWGN simulated;','markersize', 10, 'linewidth', 2)
semilogy(sim_in.EbNovec, sim_dqpsk.bervec,'b+-;DQPSK AWGN simulated;','markersize', 10, 'linewidth', 2)
xlabel('Eb/No (dB)')
ylabel('BER')
grid("minor")
axis([min(sim_in.EbNovec) max(sim_in.EbNovec) 1E-3 1])
figure (3, 'position', [400, 10, 600, 400]); clf;
semilogy(sim_in.EbNovec, ber_awgn_theory,'r+-;QPSK AWGN theory;','markersize', 10, 'linewidth', 2)
hold on;
semilogy(sim_in.EbNovec, sim_qpsk.bervec,'g+-;QPSK AWGN simulated;','markersize', 10, 'linewidth', 2)
semilogy(sim_in.EbNovec, sim_dqpsk.bervec,'b+-;DQPSK AWGN simulated;','markersize', 10, 'linewidth', 2)
semilogy(hf_sim_in.EbNovec, ber_hf_theory,'r+-;QPSK HF theory;','markersize', 10, 'linewidth', 2)
semilogy(hf_sim_in.EbNovec, sim_dqpsk_hf.bervec,'b+-;DQPSK HF simulated;','markersize', 10, 'linewidth', 2)
semilogy(hf_sim_in.EbNovec, sim_qpsk_hf.bervec,'g+-;QPSK HF simulated;','markersize', 10, 'linewidth', 2)
semilogy(hf_sim_in.EbNovec, sim_qpsk_hf_div.bervec,'c+-;QPSK Diversity HF simulated;','markersize', 10, 'linewidth', 2)
hold off;
xlabel('Eb/No (dB)')
ylabel('BER')
grid("minor")
axis([min(hf_sim_in.EbNovec) max(hf_sim_in.EbNovec) 1E-3 1])
endfunction
% -------------------------------------------------------------
more off;
rand('seed',1); randn('seed', 1);
run_curves
#run_single

142
octave/ldpc_gen_c_h_file.m
View File

@ -1,142 +0,0 @@
% ldpc_gen_c_h_file.m
% David Rowe Sep 2015, B. Van Slyke 2019
%
% Create .c and h files for use in LDPC decoders
%
% NOTE: You'll need to install the CML library as a number of functions involved
% in LDPC use it. See ldpc.m for instructions in installing the CML
% library.
%
% usage examples:
%
% 1/ Using codes defined in external files:
%
% octave:1> ldpc_gen_c_h_file("HRA_112_112.txt")
% octave:1> ldpc_gen_c_h_file(""H_4096_8192_3d.mat")
%
% 2/ Using built in CML codes:
%
% octave:1> ldpc_gen_c_h_file("dvbs2", 0.6, 16200)
%
% Output: Two files with the same filename as the LDPC input, but with .c and .h
% extensions.
function ldpc_gen_c_h_file(varargin)
ldpc % load ldpc functions
ldpc_fsk_lib % for ldpc_encode
% Assuming cml has been installed in the users' home folder, which is the
% default install location
init_cml();
if nargin == 0
printf("Error - you must specify a file containing the LDPC codes (e.g. HRA_112_112.txt).\n");
return;
end
loadStr = varargin{1};
max_iterations = 100;
decoder_type = 0;
% the tests are performed using BPSK modulation, but in practice codes can be used
% with other modulation, e.g. QPSK
mod_order = 2; modulation = 'BPSK'; mapping = 'gray';
if strcmp(loadStr, "dvbs2")
rate = varargin{2};
framesize = varargin{3};
code_param = ldpc_init_builtin(loadStr, rate, framesize, modulation, mod_order, mapping);
n = code_param.ldpc_coded_bits_per_frame;
k = code_param.ldpc_data_bits_per_frame;
ldpcArrayName = sprintf("H_%d_%d",n,k);
includeFileName = strcat(ldpcArrayName, '.h');
sourceFileName = strcat(ldpcArrayName, '.c');
else
% The ldpc variable name may not be what we want for a file/variable names, but
% the load filename will be, so use it.
[~,ldpcArrayName,ext] = fileparts(loadStr);
includeFileName = strcat(ldpcArrayName, '.h');
sourceFileName = strcat(ldpcArrayName, '.c');
% Get the ext of the file first. If it's a txt, then do what we
% are doing. If .mat, then just load, knowing the variable is HRA
if strcmp(ext, '.mat') == 1
load(loadStr);
if exist("H") & !exist("HRA")
printf("renaming H to HRA...\n");
HRA=H;
end
else
% When calling 'load' this way, it returns a struct. The code assumes the
% struct has one element, and the one/first element is the array
% to process
tempStruct = load(loadStr);
b = fieldnames(tempStruct);
ldpcArrayName = b{1,1};
% extract the array from the struct
HRA = tempStruct.(ldpcArrayName);
endif
code_param = ldpc_init_user(HRA, modulation, mod_order, mapping);
end
code_length = code_param.coded_syms_per_frame;
% First, create the H file
f = fopen(includeFileName, "wt");
printHeader(f, includeFileName, ldpcArrayName, mfilename());
fprintf(f,"#define %s_NUMBERPARITYBITS %d\n", ldpcArrayName, rows(code_param.H_rows));
fprintf(f,"#define %s_MAX_ROW_WEIGHT %d\n", ldpcArrayName, columns(code_param.H_rows));
fprintf(f,"#define %s_CODELENGTH %d\n", ldpcArrayName, code_param.coded_syms_per_frame);
fprintf(f,"#define %s_NUMBERROWSHCOLS %d\n", ldpcArrayName, rows(code_param.H_cols));
fprintf(f,"#define %s_MAX_COL_WEIGHT %d\n", ldpcArrayName, columns(code_param.H_cols));
fprintf(f,"#define %s_DEC_TYPE %d\n", ldpcArrayName, decoder_type);
fprintf(f,"#define %s_MAX_ITER %d\n", ldpcArrayName, max_iterations);
fprintf(f,"\n");
fprintf(f,"extern const uint16_t %s_H_rows[];\n", ldpcArrayName);
fprintf(f,"extern const uint16_t %s_H_cols[];\n", ldpcArrayName);
fclose(f);
% Then, the C file
f = fopen(sourceFileName, "wt");
printHeader(f, sourceFileName, ldpcArrayName, mfilename());
fprintf(f, "#include <stdint.h>\n");
fprintf(f, "#include \"%s\"\n", includeFileName);
% clock out 2D array to linear C array in row order ....
fprintf(f,"\nconst uint16_t %s_H_rows[] = {\n", ldpcArrayName);
[r c] = size(code_param.H_rows);
for j=1:c
for i=1:r
fprintf(f, "%d", code_param.H_rows(i,j));
if (i == r) && (j ==c) % weird, this does nothing
fprintf(f,"\n};\n");
else
fprintf(f,", ");
end
end
end
fprintf(f,"\nconst uint16_t %s_H_cols[] = {\n", ldpcArrayName);
[r c] = size(code_param.H_cols);
for j=1:c
for i=1:r
fprintf(f, "%d", code_param.H_cols(i,j));
if (i == r) && (j == c)
fprintf(f,"\n};\n");
else
fprintf(f,", ");
end
end
end
fclose(f);
endfunction
function printHeader(f, includeFileName, ldpcArrayName, mFilename)
fprintf(f, "/*\n FILE....: %s\n\n", includeFileName);
fprintf(f, " Static arrays for LDPC codec %s, generated by %s.m.\n*/\n\n", ldpcArrayName, mFilename);
endfunction

72
octave/ldpc_gen_h_file.m
View File

@ -1,72 +0,0 @@
% ldpc_gen_h_file.m
% David Rowe Sep 2015
%
% Create a C include file for use in mpdecode.c C cmd line LDPC decoder
function ldpc_gen_h_file(code_param, max_iterations, decoder_type, input_decoder_c, x_hat, detected_data)
f = fopen(code_param.c_include_file, "wt");
fprintf(f, "/*\n FILE....: %s\n\n Static arrays for LDPC codec, generated", code_param.c_include_file);
fprintf(f, "\n ldpc_gen_h_file.m.\n\n*/\n\n");
fprintf(f,"#define NUMBERPARITYBITS %d\n", rows(code_param.H_rows));
fprintf(f,"#define MAX_ROW_WEIGHT %d\n", columns(code_param.H_rows));
fprintf(f,"#define CODELENGTH %d\n", code_param.symbols_per_frame);
fprintf(f,"#define NUMBERROWSHCOLS %d\n", rows(code_param.H_cols));
fprintf(f,"#define MAX_COL_WEIGHT %d\n", columns(code_param.H_cols));
fprintf(f,"#define DEC_TYPE %d\n", decoder_type);
fprintf(f,"#define MAX_ITER %d\n", max_iterations);
fprintf(f,"\ndouble H_rows[] = {\n");
% clock out 2D array to linear C array in row order ....
[r c] = size(code_param.H_rows);
for j=1:c
for i=1:r
fprintf(f, "%d", code_param.H_rows(i,j));
if (i == r) && (j ==c)
fprintf(f,"\n};\n");
else
fprintf(f,", ");
end
end
end
fprintf(f,"\ndouble H_cols[] = {\n");
[r c] = size(code_param.H_cols);
for j=1:c
for i=1:r
fprintf(f, "%d", code_param.H_cols(i,j));
if (i == r) && (j == c)
fprintf(f,"\n};\n");
else
fprintf(f,", ");
end
end
end
fprintf(f,"\ndouble input[] = {\n");
for i=1:length(input_decoder_c)
fprintf(f, "%.17g", input_decoder_c(i));
if i == length(input_decoder_c)
fprintf(f,"\n};\n");
else
fprintf(f,", ");
end
end
fprintf(f,"\nchar detected_data[] = {\n");
for i=1:length(detected_data)
fprintf(f, "%d", detected_data(i));
if i == length(detected_data)
fprintf(f,"\n};\n");
else
fprintf(f,", ");
end
end
fclose(f);
end

14
octave/load_f32.m
View File

@ -1,14 +0,0 @@
% load_f32.m
% David Rowe Jan 2019
%
% load up .f32 binary files from dump_data
function features = load_f32(fn, ncols)
f=fopen(fn,"rb");
features_lin=fread(f, 'float32');
fclose(f);
nrows = length(features_lin)/ncols;
features = reshape(features_lin, ncols, nrows);
features = features';
endfunction

60
octave/make_hilb.m
View File

@ -1,60 +0,0 @@
% make_hilb.m
% David Rowe May 2015
%
% creates Hilber Transformer FIR coeffs
graphics_toolkit ("gnuplot");
% from https://www.dsprelated.com/freebooks/sasp/Hilbert_Transform_Design_Example.html
M = 257; % window length = FIR filter length (Window Method)
fs = 8000; % sampling rate assumed (Hz)
f1 = 100; % lower pass-band limit = transition bandwidth (Hz)
beta = 8; % beta for Kaiser window for decent side-lobe rejection
fn = fs/2; % Nyquist limit (Hz)
f2 = fn - f1; % upper pass-band limit
N = 2^(nextpow2(8*M)); % large FFT for interpolated display
k1 = round(N*f1/fs); % lower band edge in bins
if k1<2, k1=2; end; % cannot have dc or fn response
kn = N/2 + 1; % bin index at Nyquist limit (1-based)
k2 = kn-k1+1; % high-frequency band edge
f1 = k1*fs/N % quantized band-edge frequencies
f2 = k2*fs/N
w = kaiser(M,beta)'; % Kaiser window in "linear phase form"
H = [ ([0:k1-2]/(k1-1)).^8,ones(1,k2-k1+1),...
([k1-2:-1:0]/(k1-1)).^8, zeros(1,N/2-1)];
h = ifft(H); % desired impulse response
hodd = imag(h(1:2:N)); % This should be zero
% put window in zero-phase form:
wzp = [w((M+1)/2:M), zeros(1,N-M), w(1:(M-1)/2)];
hw = wzp .* h; % single-sideband FIR filter, zero-centered
Hw = fft(hw);
hh = [hw(N-(M-1)/2+1:N),hw(1:(M+1)/2)]; % causal FIR
hh *= 2;
figure(1);
HH = fft([hh,zeros(1,N-M)]);
plot(20*log10(abs(HH)));
figure(2);
subplot(211); plot(real(hh)); title('real imp resp');
subplot(212); plot(imag(hh)); title('imag imp resp');
% save coeffs to a C header file
f=fopen("../src/ht_coeff.h","wt");
fprintf(f,"/* Hilbert Transform FIR filter coeffs */\n");
fprintf(f,"/* Generated by make_hilb Octave script */\n");
fprintf(f,"\n#define HT_N %d\n\n", M);
fprintf(f,"COMP ht_coeff[]={\n");
for r=1:M
if r < M
fprintf(f, " {%f,%f},\n", real(hh(r)), imag(hh(r)));
else
fprintf(f, " {%f,%f}\n};", real(hh(r)), imag(hh(r)));
end
end
fclose(f);

39
octave/make_ssbfilt.m
View File

@ -1,39 +0,0 @@
% make_ssbfilt.m
% David Rowe May 2015
%
% Creates low pass filter coeff used to implement a SSB filter in ch
graphics_toolkit ("gnuplot");
ssbfilt_n = 100;
ssbfilt_bw = 2400;
ssbfilt_centre = 1500;
Fs = 8000;
ssbfilt_coeff = sbfilt_coeff = fir1(ssbfilt_n, ssbfilt_bw/Fs);
figure(1)
clf;
h = freqz(ssbfilt_coeff,1,Fs/2);
plot(20*log10(abs(h)))
grid minor
% save coeffs to a C header file
f=fopen("../src/ssbfilt_coeff.h","wt");
fprintf(f,"/* %d Hz LPF FIR filter coeffs */\n", ssbfilt_bw);
fprintf(f,"/* Generated by make_ssbfilt Octave script */\n");
fprintf(f,"\n#define SSBFILT_N %d\n\n", ssbfilt_n);
fprintf(f,"\n#define SSBFILT_CENTRE %d\n\n", ssbfilt_centre);
fprintf(f,"float ssbfilt_coeff[]={\n");
for r=1:ssbfilt_n
if r < ssbfilt_n
fprintf(f, " %f,\n", ssbfilt_coeff(r));
else
fprintf(f, " %f\n};", ssbfilt_coeff(r));
end
end
fclose(f);

13
octave/pl_scatter.m
View File

@ -1,13 +0,0 @@
% pl_scatter.m
% Render scatter plot from freedv_data_raw_rx --scatter
function pl_scatter(filename)
s=load(filename);
figure(1); clf;
for b=1:length(fieldnames(s))
field_name = fieldnames(s){b};
x = s.(field_name);
plot(x,'+');
end
print("scatter.png", "-dpng");
endfunction

141
octave/qam16_test.m
View File

@ -1,141 +0,0 @@
% qam16_test.m
% David Rowe May 2020
%
% Octave symbol rate QAM16/LDPC experiments
% Libraries we need
1;
qam16;
ldpc;
function test_qam16(fg=2)
printf("QAM16 ----------------------------------------\n");
mod_order = 16; bps = log2(mod_order);
modulation = 'QAM'; mapping = ""; demod_type = 0; decoder_type = 0;
max_iterations = 100; EsNo_dec = 10;
qam16_const = [
1 + j, 1 + j*3, 3 + j, 3 + j*3;
1 - j, 1 - j*3, 3 - j, 3 - j*3;
-1 + j, -1 + j*3, -3 + j, -3 + j*3;
-1 - j, -1 - j*3, -3 - j, -3 - j*3];
rms = sqrt(qam16_const(:)'*qam16_const(:)/16);
qam16_const = qam16_const/rms;
constellation_source = 'custom';
test_qam16_mod_demod(qam16_const);
load HRA_504_396.txt
if strcmp(constellation_source,'cml')
code_param = ldpc_init_user(HRA_504_396, modulation, mod_order, mapping);
else
code_param = ldpc_init_user(HRA_504_396, modulation, mod_order, mapping, reshape(qam16_const,1,16));
end
rate = code_param.ldpc_data_bits_per_frame/code_param.ldpc_coded_bits_per_frame;
printf("EbNodB Tbits Terrs BER Tcbits Tcerrs Perrs CBER CPER\n");
EbNodBvec = 3:10; Ntrials = 1000;
for i=1:length(EbNodBvec)
EbNodB = EbNodBvec(i);
EsNodB = EbNodB + 10*log10(rate) + 10*log10(bps); EsNodBvec(i) = EsNodB;
EsNo = 10^(EsNodB/10);
variance = 1/EsNo;
Terrs = Tbits = 0; Tcerrs = 0; Tcbits = 0; Perrs = 0; rx_symbols_log = [];
for nn = 1:Ntrials
tx_bits = round(rand(1, code_param.ldpc_data_bits_per_frame));
[tx_codeword, tx_symbols] = ldpc_enc(tx_bits, code_param);
noise = sqrt(variance*0.5)*(randn(1,length(tx_symbols)) + j*randn(1,length(tx_symbols)));
rx_symbols = tx_symbols + noise;
rx_symbols_log = [rx_symbols_log rx_symbols];
% uncoded decode/demod and count errors
rx_codeword = zeros(1,code_param.ldpc_coded_bits_per_frame);
for s=1:length(rx_symbols)
rx_codeword((s-1)*bps+1:s*bps) = qam16_demod(qam16_const,rx_symbols(s));
end
Nerr = sum(xor(tx_codeword,rx_codeword));
Terrs += Nerr;
Tbits += code_param.ldpc_coded_bits_per_frame;
% LDPC demod/decode and count errors
dec_rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols, EsNo_dec, ones(1,length(rx_symbols)));
errors_positions = xor(tx_bits, dec_rx_codeword(1:code_param.ldpc_data_bits_per_frame));
Ncerr = sum(errors_positions);
Tcbits += code_param.ldpc_data_bits_per_frame; Tcerrs += Ncerr;
if Ncerr Perrs++; end
end
figure(fg); clf; plot(rx_symbols_log,"."); axis([-1.5 1.5 -1.5 1.5]); drawnow;
printf("%5.1f %6d %6d %5.2f %6d %6d %6d %5.2f %5.2f\n",
EbNodB, Tbits, Terrs, Terrs/Tbits, Tcbits, Tcerrs, Perrs, Tcerrs/Tcbits, Perrs/Ntrials);
ber(i) = Terrs/Tbits; cber(i) = Tcerrs/Tcbits; cper(i) = Perrs/Ntrials;
end
print("qam64_scatter.png","-dpng");
figure(fg+1); clf; title('QAM16 Uncoded');
uncoded_EbNodBvec = EbNodBvec + 10*log10(rate);
ber_theory = ber_qam(uncoded_EbNodBvec);
semilogy(uncoded_EbNodBvec,ber_theory,'b+-;uncoded QAM16 BER theory;','markersize', 10, 'linewidth', 2); hold on;
semilogy(uncoded_EbNodBvec,ber+1E-10,'g+-;uncoded QAM16 BER;','markersize', 10, 'linewidth', 2); hold on;
grid; axis([min(uncoded_EbNodBvec) max(uncoded_EbNodBvec) 1E-5 1]); xlabel('Uncoded Eb/No (dB)');
print("qam16_uncoded_ber.png","-dpng");
figure(fg+2); clf; title('QAM16 with LDPC (504,396)');
semilogy(EbNodBvec,cber+1E-10,'b+-;QAM16 coded BER;','markersize', 10, 'linewidth', 2); hold on;
semilogy(EbNodBvec,cper+1E-10,'g+-;QAM16 coded PER;','markersize', 10, 'linewidth', 2); hold off;
grid; axis([min(EbNodBvec) max(EbNodBvec) 1E-5 1]); xlabel('Eb/No (dB)');
figure(fg+3); clf; title('QAM16 with LDPC (504,396)');
semilogy(EsNodBvec,cber+1E-10,'b+-;QAM16 coded BER;','markersize', 10, 'linewidth', 2); hold on;
semilogy(EsNodBvec,cper+1E-10,'g+-;QAM16 coded PER;','markersize', 10, 'linewidth', 2); hold off;
grid; axis([min(EsNodBvec) max(EsNodBvec) 1E-5 1]); xlabel('Es/No (dB)');
print("qam16_504_396.png","-dpng");
endfunction
% Thanks Bill VK5DSP, for the QAM BER functions
function p = ber_qam(ebn0)
% Calculate the bit error rate (BER) for square 16QAM in AWGN
% given the Eb/N0 in dB, ebn0 can be a scalar or vector
% (assuming coherent detection, uncoded)
% [section 5.3 Webb and Hanzo text]
e = 4*10.^(ebn0/10); % Es/N0 vector in linear
b2 = qfn(sqrt(e/5));
b1 = (qfn(sqrt(e/5)) + qfn(3*sqrt(e/5)))/2;
p = (b1+b2)/2;
endfunction
function tail=qfn(a)
% Usage: tail=qfn(a)
% where: tail=area under the tail of the normal dist. from a to inf.
% for zero mean, unit variance distribution
%
% If no argument is given, plot Q(x) for x = 0 to 5
% use erfc instead of 1-erf to avoid truncation errors! April 2010
fact = 1 / sqrt(2);
if exist('a')
% tail = 0.5 * ( 1 - erf(a * fact));
tail = 0.5 * erfc(a * fact);
else
x=(0: 0.1: 6); semilogy(x, 0.5*( erfc(x * fact)));
title('Q function plot');
xlabel('x'); ylabel('Q(x)');
end
endfunction
% --------------------------------------------------------------------------------
% START SIMULATIONS
% --------------------------------------------------------------------------------
more off;
format;
% Start CML library (see CML set up instructions in ldpc.m)
init_cml();
test_qam16(1)

851
octave/test_ldpc_fsk_lib.m
View File

@ -1,851 +0,0 @@
% test_ldpc_fsk_lib
% David Rowe 16 April 2016
%
% A series of tests for ldpc_fsk_lib, and C versions ldpc_enc and ldpc_dec.
% Gradually builds up complete C command line for SSTV balloon system,
% using Octave versions of LDPC and FSK modem as reference points.
1;
% encodes and decodes one frame, also writes codeword.bin for testing
% decode_from_file() below, and can optionally generate include file for
% C version of decoder.
function [data code_param] = simple_ut(c_include_file)
load('H2064_516_sparse.mat');
HRA = full(HRA);
max_iterations = 100;
decoder_type = 0;
EsNodB = 3;
mod_order = 2;
code_param = ldpc_init(HRA, mod_order);
data = round( rand( 1, code_param.data_bits_per_frame ) );
codeword = ldpc_encode(code_param, data);
f = fopen("codeword.bin","wt"); fwrite(f, codeword, "uint8"); fclose(f);
s = 1 - 2 * codeword;
code_param.symbols_per_frame = length( s );
EsNo = 10^(EsNodB/10);
variance = 1/(2*EsNo);
noise = sqrt(variance)* randn(1,code_param.symbols_per_frame);
rx = s + noise;
if nargin == 1
code_param.c_include_file = c_include_file;
end
[detected_data Niters] = ldpc_decode(rx, code_param, max_iterations, decoder_type);
error_positions = xor(detected_data(1:code_param.data_bits_per_frame), data);
Nerrs = sum(error_positions);
printf("Nerrs = %d\n", Nerrs);
end
% This version decodes from a file of bits
function detected_data = decode_from_file(filename)
max_iterations = 100;
decoder_type = 0;
load('H2064_516_sparse.mat');
HRA = full(HRA);
mod_order = 2;
f = fopen(filename,"rb"); codeword = fread(f, "uint8")'; fclose(f);
r = 1 - 2 * codeword;
code_param = ldpc_init(HRA, mod_order);
[detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type);
end
% plots a BER curve for the LDPC decoder. Takes a while to run, uses parallel cores
function plot_curve
num_cores = 4; % set this to the number of cores you have
load('H2064_516_sparse.mat');
HRA = full(HRA);
[Nr Nc] = size(HRA);
sim_in.rate = (Nc-Nr)/Nc;
sim_in.HRA = HRA;
sim_in.mod_order = 2;
sim_in.framesize = Nc;
sim_in.mod_order = 2;
sim_in.Lim_Ferrs = 100;
% note we increase number of trials as BER goes down
Esvec = [ 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 ];
Ntrials = [ 1E4 1E4 1E4 1E4 1E5 1E5 1E5 1E5 1E5 ];
num_runs = length(Esvec)
sim_in_vec(1:num_runs) = sim_in;
for i = 1:num_runs
sim_in_vec(i).Esvec = Esvec(i);
sim_in_vec(i).Ntrials = Ntrials(i);
end
%sim_out = ldpc5(sim_in_vec(1));
tstart = time();
sim_out = pararrayfun(num_cores, @ldpc5, sim_in_vec);
tend = time();
total_bits = sum(Ntrials)*sim_in.framesize;
total_secs = tend - tstart;
printf("%d bits in %4.1f secs, or %5f bits/s\n", total_bits, total_secs, total_bits/total_secs);
for i=1:num_runs
Ebvec(i) = sim_out(i).Ebvec;
BERvec(i) = sim_out(i).BERvec;
end
semilogy(Ebvec, BERvec, '+-')
xlabel('Eb/N0')
ylabel('BER')
title(['H2064 516 sparse.mat' ' ' date])
end
% Test C encoder
function test_c_encoder
load('H2064_516_sparse.mat');
HRA = full(HRA);
max_iterations = 100;
decoder_type = 0;
EsNodB = 3;
mod_order = 2;
frames = 100;
EsNo = 10^(EsNodB/10);
variance = 1/(2*EsNo);
code_param = ldpc_init(HRA, mod_order);
data = round(rand(1,frames*code_param.data_bits_per_frame));
f = fopen("data.bin","wt"); fwrite(f, data, "uint8"); fclose(f);
% Outboard C encoder
system("../src/ldpc_enc data.bin codewords.bin");
% Test with Octave decoder
f = fopen("codewords.bin","rb"); codewords = fread(f, "uint8")'; fclose(f);
Nerrs = 0;
for i=1:frames
st = (i-1)*code_param.symbols_per_frame+1; en = st+code_param.symbols_per_frame-1;
tx = 1 - 2 * codewords(st:en);
noise = sqrt(variance)*randn(1,code_param.symbols_per_frame);
rx = tx + noise;
[detected_data Niters] = ldpc_decode(rx, code_param, max_iterations, decoder_type);
st = (i-1)*code_param.data_bits_per_frame+1; en = st+code_param.data_bits_per_frame-1;
error_positions = xor(detected_data(1:code_param.data_bits_per_frame), data(st:en));
Nerrs += sum(error_positions);
end
printf("Nerrs = %d\n", Nerrs);
end
function test_c_decoder
load('H2064_516_sparse.mat');
HRA = full(HRA);
max_iterations = 100;
decoder_type = 0;
mod_order = 2;
frames = 10;
EsNodB = 2;
sdinput = 1;
EsNo = 10^(EsNodB/10);
variance = 1/(2*EsNo);
code_param = ldpc_init(HRA, mod_order);
data = round(rand(1,code_param.data_bits_per_frame*frames));
f = fopen("data.bin","wt"); fwrite(f, data, "uint8"); fclose(f);
system("../src/ldpc_enc data.bin codewords.bin");
f = fopen("codewords.bin","rb"); codewords = fread(f, "uint8")'; fclose(f);
s = 1 - 2 * codewords;
noise = sqrt(variance)*randn(1,code_param.symbols_per_frame*frames);
r = s + noise;
% calc LLRs frame by frame
for i=1:frames
st = (i-1)*code_param.symbols_per_frame+1;
en = st + code_param.symbols_per_frame-1;
llr(st:en) = sd_to_llr(r(st:en));
end
% Outboard C decoder
if sdinput
f = fopen("sd.bin","wb"); fwrite(f, r, "double"); fclose(f);
system("../src/ldpc_dec sd.bin data_out.bin --sd");
else
f = fopen("llr.bin","wb"); fwrite(f, llr, "double"); fclose(f);
system("../src/ldpc_dec llr.bin data_out.bin");
end
f = fopen("data_out.bin","rb"); data_out = fread(f, "uint8")'; fclose(f);
Nerrs = Nerrs2 = zeros(1,frames);
for i=1:frames
% Check C decoder
data_st = (i-1)*code_param.data_bits_per_frame+1;
data_en = data_st+code_param.data_bits_per_frame-1;
st = (i-1)*code_param.symbols_per_frame+1;
en = st+code_param.data_bits_per_frame-1;
data_out_c = data_out(st:en);
error_positions = xor(data_out_c, data(data_st:data_en));
Nerrs(i) = sum(error_positions);
% Octave decoder
st = (i-1)*code_param.symbols_per_frame+1; en = st+code_param.symbols_per_frame-1;
[detected_data Niters] = ldpc_decode(r(st:en), code_param, max_iterations, decoder_type);
st = (i-1)*code_param.data_bits_per_frame+1; en = st+code_param.data_bits_per_frame-1;
data_out_octave = detected_data(1:code_param.data_bits_per_frame);
error_positions = xor(data_out_octave, data(st:en));
Nerrs2(i) = sum(error_positions);
%printf("%4d ", Niters);
end
printf("Errors per frame:\nC.....:");
for i=1:frames
printf("%4d ", Nerrs(i));
end
printf("\nOctave:");
for i=1:frames
printf("%4d ", Nerrs2(i));
end
printf("\n");
end
% Saves a complex vector s to a file "filename" of IQ unsigned 8 bit
% chars, same as RTLSDR format.
function save_rtlsdr(filename, s)
mx = max(abs(s));
re = real(s); im = imag(s);
l = length(s);
iq = zeros(1,2*l);
%iq(1:2:2*l) = 127 + re*(127/mx);
%iq(2:2:2*l) = 127 + im*(127/mx);
iq(1:2:2*l) = 127 + 32*re;
iq(2:2:2*l) = 127 + 32*im;
figure(3); clf; plot(iq); title('simulated IQ signal from RTL SDR');
fs = fopen(filename,"wb");
fwrite(fs,iq,"uchar");
fclose(fs);
endfunction
% Oversamples by a factor of 2 using Octaves resample() function then
% uses linear interpolation to achieve fractional sample rate
function rx_resample_fract = fractional_resample(rx, resample_rate);
assert(resample_rate < 2, "keep resample_rate between 0 and 2");
rx_resample2 = resample(rx, 2, 1);
l = length(rx_resample2);
rx_resample_fract = zeros(1,l);
k = 1;
step = 2/resample_rate;
for i=1:step:l-1
i_low = floor(i);
i_high = ceil(i);
f = i - i_low;
rx_resample_fract(k) = (1-f)*rx_resample2(i_low) + f*rx_resample2(i_high);
%printf("i: %f i_low: %d i_high: %d f: %f\n", i, i_low, i_high, f);
k++;
end
rx_resample_fract = rx_resample_fract(1:k-1);
endfunction
% Using simulated SSTV packet, generate complex fsk mod signals, 8-bit
% unsigned IQ for feeding into C demod chain. Can also be used to
% generate BER curves. Found bugs in UW size and our use of csdr
% re-sampler using this function, and by gradually and carefully
% building up the C command line.
#{
todo: [X] uncoded BER
[X] octave fsk demod
[X] use C demod
[X] compare uncoded BER to unsigned 8 bit IQ to regular 16-bit
[X] generate complex rx signal with noise
[X] used cmd line utils to drive demod
[X] test with resampler
[X] measure effect on PER with coding
#}
function [n_uncoded_errs n_uncoded_bits] = run_sstv_sim(sim_in, EbNodB)
frames = sim_in.frames;
demod_type = sim_in.demod_type;
% init LDPC code
load('H2064_516_sparse.mat');
HRA = full(HRA);
max_iterations = 100;
decoder_type = 0;
mod_order = 2;
code_param = ldpc_init(HRA, mod_order);
% note fixed frame of bits used for BER testing
tx_codeword = gen_sstv_frame;
% init FSK modem
fsk_horus_as_a_lib = 1;
fsk_horus;
states = fsk_horus_init_hbr(9600, 8, 1200, 2, length(tx_codeword));
states.df(1:states.M) = 0;
states.dA(1:states.M) = 1;
states.tx_real = 0; % Octave fsk_mod generates complex valued output
% so we can simulate rtl_sdr complex output
% Set up simulated tx tones to sit in the middle of cdsr passband
filt_low_norm = 0.1; filt_high_norm = 0.4;
fc = states.Fs*(filt_low_norm + filt_high_norm)/2;
%fc = 1800;
f1 = fc - states.Rs/2;
f2 = fc + states.Rs/2;
states.ftx = [f1 f2];
% set up AWGN channel
EbNo = 10^(EbNodB/10);
variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol);
% start simulation ----------------------------------------
tx_bit_stream = [];
for i=1:frames
% uncomment for different data on each frame
%tx_codeword = gen_sstv_frame;
tx_bit_stream = [tx_bit_stream tx_codeword];
end
printf("%d bits at %d bit/s is a %3.1f second run\n", length(tx_bit_stream), 115200,length(tx_bit_stream)/115200);
% modulate and channel model
tx = fsk_horus_mod(states, tx_bit_stream);
noise_real = sqrt(variance)*randn(length(tx),1);
noise_complex = sqrt(variance/2)*(randn(length(tx),1) + j*randn(length(tx),1));
% demodulate -----------------------------------------------------
if demod_type == 1
% Octave demod
if states.tx_real
rx = tx + noise_real;
else
rx = tx + noise_complex;
end
SNRdB = 10*log10(var(tx)/var(noise_complex));
% demodulate frame by frame using Octave demod
st = 1;
run_frames = floor(length(rx)/states.N);
rx_bit_stream = [];
rx_sd_stream = [];
for f=1:run_frames
% extract nin samples from rx sample stream
nin = states.nin;
en = st + states.nin - 1;
if en <= length(rx) % due to nin variations its possible to overrun buffer
sf = rx(st:en);
st += nin;
% demodulate to stream of bits
states.f = [f1 f2]; % note that for Octave demod we cheat and use known tone frequencies
% allows us to determine if freq offset estimation in C demod is a problem
[rx_bits states] = fsk_horus_demod(states, sf);
rx_bit_stream = [rx_bit_stream rx_bits];
rx_sd_stream = [rx_sd_stream states.rx_bits_sd];
end
end
end
if demod_type == 2
% baseline C demod
if states.tx_real
rx = tx + noise_real;
else
rx = 2*real(tx) + noise_real;
end
SNRdB = 10*log10(var(tx)/var(noise_real));
rx_scaled = 1000*real(rx);
f = fopen("fsk_demod.raw","wb"); fwrite(f, rx_scaled, "short"); fclose(f);
system("../build_linux/src/fsk_demod 2X 8 9600 1200 fsk_demod.raw fsk_demod.bin");
f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f);
end
if demod_type == 3
% C demod driven by csdr command line kung fu
assert(states.tx_real == 0, "need complex signal for this test");
rx = tx + noise_complex;
SNRdB = 10*log10(var(tx)/var(noise_real));
save_rtlsdr("fsk_demod.iq", rx);
system("cat fsk_demod.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin");
f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f);
end
if demod_type == 4
% C demod with resampler ....... getting closer to Mark's real time cmd line
assert(states.tx_real == 0, "need complex signal for this test");
rx = tx + noise_complex;
SNRdB = 10*log10(var(tx)/var(noise_real));
printf("resampling ...\n");
rx_resample_fract = fractional_resample(rx, 1.08331);
%rx_resample_fract = fractional_resample(rx_resample_fract, 1/1.08331);
save_rtlsdr("fsk_demod_resample.iq", rx_resample_fract);
printf("run C cmd line chain ...\n");
% system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin");
system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../unittest/tsrc - - 0.9230968 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin");
% system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr fractional_decimator_ff 1.08331 | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin");
f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f);
end
if demod_type == 5
% C demod with resampler and use C code to measure PER, in this
% test we don't need to run state machine below as C code gives us
% the outputs we need
assert(states.tx_real == 0, "need complex signal for this test");
rx = tx + noise_complex;
SNRdB = 10*log10(var(tx)/var(noise_real));
printf("fract resampling ...\n");
rx_resample_fract = fractional_resample(rx, 1.08331);
save_rtlsdr("fsk_demod_resample.iq", rx_resample_fract);
% useful for HackRF
%printf("10X resampling ...\n");
%rx_resample_10M = resample(rx_resample_fract, 10, 1);
%save_rtlsdr("fsk_demod_10M.iq", rx_resample_10M);
printf("run C cmd line chain - uncoded PER\n");
system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../unittest/tsrc - - 0.9230968 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - - | ../src/drs232 - /dev/null -v");
printf("run C cmd line chain - LDPC coded PER\n");
system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../unittest/tsrc - - 0.9230968 | ../build_linux/src/fsk_demod 2XS 8 9600 1200 - - | ../src/drs232_ldpc - /dev/null -v");
end
if demod_type == 6
% C demod with complex input driven simplfied csdr command line, just measure BER of demod
assert(states.tx_real == 0, "need complex signal for this test");
rx = tx + noise_complex;
SNRdB = 10*log10(var(tx)/var(noise_real));
save_rtlsdr("fsk_demod.iq", rx);
system("cat fsk_demod.iq | csdr convert_u8_f | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin C");
f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f);
end
if demod_type == 7
% C demod with complex input, measure uncoded and uncoded PER
assert(states.tx_real == 0, "need complex signal for this test");
rx = tx + noise_complex;
SNRdB = 10*log10(var(tx)/var(noise_real));
save_rtlsdr("fsk_demod.iq", rx);
printf("run C cmd line chain - uncoded PER\n");
system("cat fsk_demod.iq | csdr convert_u8_f | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - - C | ../src/drs232 - /dev/null -v");
printf("run C cmd line chain - LDPC coded PER\n");
%system("cat fsk_demod.iq | csdr convert_u8_f | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2XS 8 9600 1200 - - C | ../src/drs232_ldpc - /dev/null -v");
system("cat fsk_demod.iq | ../build_linux/src/fsk_demod 2XS 8 9600 1200 - - CU8 | ../src/drs232_ldpc - /dev/null -v");
end
if (demod_type != 5) && (demod_type != 7)
% state machine. Look for SSTV UW. When found count bit errors over one frame of bits
state = "wait for uw";
start_uw_ind = 16*10+1; end_uw_ind = start_uw_ind + 5*10 - 1;
uw_rs232 = tx_codeword(start_uw_ind:end_uw_ind); luw = length(uw_rs232);
start_frame_ind = end_uw_ind + 1;
nbits = length(rx_bit_stream);
uw_thresh = 5;
n_uncoded_errs = 0;
n_uncoded_bits = 0;
n_packets_rx = 0;
last_i = 0;
% might as well include RS232 framing bits in uncoded error count
nbits_frame = code_param.data_bits_per_frame*10/8;
uw_errs = zeros(1, nbits);
for i=luw:nbits
uw_errs(i) = sum(xor(rx_bit_stream(i-luw+1:i), uw_rs232));
end
for i=luw:nbits
next_state = state;
if strcmp(state, 'wait for uw')
if uw_errs(i) <= uw_thresh
next_state = 'count errors';
tx_frame_ind = start_frame_ind;
rx_frame_ind = i + 1;
n_uncoded_errs_this_frame = 0;
%printf("%d %s %s\n", i, state, next_state);
if last_i
printf("i: %d i-last_i: %d ", i, i-last_i);
end
end
end
if strcmp(state, 'count errors')
n_uncoded_errs_this_frame += xor(rx_bit_stream(i), tx_codeword(tx_frame_ind));
n_uncoded_bits++;
tx_frame_ind++;
if tx_frame_ind == (start_frame_ind+nbits_frame)
n_uncoded_errs += n_uncoded_errs_this_frame;
printf("n_uncoded_errs_this_frame: %d\n", n_uncoded_errs_this_frame);
frame_rx232_rx = rx_bit_stream(rx_frame_ind:rx_frame_ind+nbits_frame-1);
%tx_codeword(start_frame_ind+1:start_frame_ind+10)
%frame_rx232_rx(1:10)
sstv_checksum(frame_rx232_rx);
last_i = i;
n_packets_rx++;
next_state = 'wait for uw';
end
end
state = next_state;
end
uncoded_ber = n_uncoded_errs/n_uncoded_bits;
printf("EbNodB: %4.1f SNRdB: %4.1f pkts: %d bits: %d errs: %d BER: %4.3f\n",
EbNodB, SNRdB, n_packets_rx, n_uncoded_bits, n_uncoded_errs, uncoded_ber);
figure(2);
plot(uw_errs);
title('Unique Word Hamming Distance')
end
endfunction
% Function to test flight mode software. Takes a rx stream of
% demodulated bits, and locates frames using UW detection. Extracts
% data and parity bits. Uses data bits to generate parity bits here
% and compare.
function compare_parity_bits(rx_bit_stream)
nframes = 500;
% init LDPC code
load('H2064_516_sparse.mat');
HRA = full(HRA);
max_iterations = 100;
decoder_type = 0;
mod_order = 2;
code_param = ldpc_init(HRA, mod_order);
% generate frame, this will have random bits not related to
% rx_stream, however we just use it for the UW
tx_codeword = gen_sstv_frame;
l = length(tx_codeword);
printf("expected rs232 frames codeword length: %d\n", l);
% state machine. Look for SSTV UW. When found count bit errors over one frame of bits
state = "wait for uw";
start_uw_ind = 16*10+1; end_uw_ind = start_uw_ind + 4*10 - 1;
uw_rs232 = tx_codeword(start_uw_ind:end_uw_ind); luw = length(uw_rs232);
start_frame_ind = end_uw_ind + 1;
nbits = nframes*l;
uw_thresh = 5;
n_uncoded_errs = 0;
n_uncoded_bits = 0;
n_packets_rx = 0;
last_i = 0;
% might as well include RS232 framing bits in uncoded error count
uw_errs = luw*ones(1, nbits);
for i=luw:nbits
uw_errs(i) = sum(xor(rx_bit_stream(i-luw+1:i), uw_rs232));
end
frame_start = find(uw_errs < 2)+1;
nframes = length(frame_start)
for i=1:nframes
% double check UW OK
st_uw = frame_start(i) - luw; en_uw = frame_start(i) - 1;
uw_err_check = sum(xor(rx_bit_stream(st_uw:en_uw), uw_rs232));
%printf("uw_err_check: %d\n", uw_err_check);
% strip off rs232 start/stop bits
nbits_rs232 = (256+2+65)*10;
nbits = (256+2+65)*8;
nbits_byte = 10;
rx_codeword = zeros(1,nbits);
pdb = 1;
for k=1:nbits_byte:nbits_rs232
for l=1:8
rx_codeword(pdb) = rx_bit_stream(frame_start(i)-1+k+l);
pdb++;
end
end
assert(pdb == (nbits+1));
data_bits = rx_codeword(1:256*8);
checksum_bits = rx_codeword(256*8+1:258*8);
parity_bits = rx_codeword(258*8+1:258*8+516);
padding_bits = rx_codeword(258*8+516+1:258*8+516+1);
% stopped here as we found bug lol!
end
figure(1); clf;
plot(uw_errs);
title('Unique Word Hamming Distance')
figure(2); clf;
lframe_start = length(frame_start);
plot(frame_start(2:lframe_start)-frame_start(1:lframe_start-1));
%title('Unique Word Hamming Distance')
endfunction
% Start simulation --------------------------------------------------------
more off;
currentdir = pwd;
thiscomp = computer;
if strcmpi(thiscomp, 'MACI64')==1
if exist('CMLSimulate')==0
cd '/Users/bill/Current/Projects/DLR_FSO/Visit2013_FSO_GEO/cml'
addpath '../' % assume the source files stored here
CmlStartup % note that this is not in the cml path!
disp('added MACI64 path and run CmlStartup')
end
end
if strfind(thiscomp, 'pc-linux-gnu')==8
if exist('LdpcEncode')==0,
cd '~/tmp/cml'
CmlStartup
disp('CmlStartup has been run')
% rmpath '/home/bill/cml/mexhelp' % why is this needed?
% maybe different path order in octave cf matlab ?
end
end
cd(currentdir)
ldpc_fsk_lib;
randn('state',1);
rand('state',1);
% ------------------ select which demo/test to run here ---------------
demo = 12;
if demo == 1
printf("simple_ut....\n");
data = simple_ut;
end
if demo == 2
printf("generate C header file....\n");
data = simple_ut("../src/H2064_516_sparse.h");
end
if demo == 3
printf("decode_from_file ......\n");
data = simple_ut;
detected_data = decode_from_file("codeword.bin");
error_positions = xor( detected_data(1:length(data)), data );
Nerrs = sum(error_positions);
printf(" Nerrs = %d\n", Nerrs);
end
if demo == 4
printf("plot a curve....\n");
plot_curve;
end
if demo == 5
% generate test data and save to disk
[data code_param] = simple_ut;
f = fopen("dat_in2064.bin","wb"); fwrite(f, data, "uint8"); fclose(f);
% Outboard C encoder
system("../src/ldpc_enc dat_in2064.bin dat_op2064.bin");
% Test with Octave decoder
detected_data = decode_from_file("dat_op2064.bin");
error_positions = xor(detected_data(1:length(data)), data);
Nerrs = sum(error_positions);
printf("Nerrs = %d\n", Nerrs);
end
if demo == 6
test_c_encoder;
end
if demo == 7
test_c_decoder;
end
% generates simulated demod soft decision symbols to drive C ldpc decoder with
if demo == 8
frames = 100;
EsNodB = 3;
EsNo = 10^(EsNodB/10);
variance = 1/(2*EsNo);
frame_rs232 = [];
for i=1:frames
frame_rs232 = [frame_rs232 gen_sstv_frame];
end
% write hard decn version to disk file, useful for fsk_mod input
f = fopen("sstv.bin","wb"); fwrite(f, frame_rs232, "char"); fclose(f);
% soft decision version (with noise)
s = 1 - 2*frame_rs232;
noise = sqrt(variance)*randn(1,length(frame_rs232));
r = s + noise;
f = fopen("sstv_sd.bin","wb"); fwrite(f, r, "float32"); fclose(f);
end
if demo == 9
frames = 100;
EbNodB = 11;
frame_rs232 = [];
for i=1:frames
frame_rs232 = [frame_rs232 gen_sstv_frame];
end
% Use C FSK modulator to generate modulated signal
f = fopen("sstv.bin","wb"); fwrite(f, frame_rs232, "char"); fclose(f);
system("../build_linux/src/fsk_mod 2 9600 1200 1200 2400 sstv.bin fsk_mod.raw");
% Add some channel noise here in Octave
f = fopen("fsk_mod.raw","rb"); tx = fread(f, "short")'; fclose(f); tx_pwr = var(tx);
Fs = 9600; Rs=1200; EbNolin = 10 ^ (EbNodB/10);
variance = (tx_pwr/2)*states.Fs/(states.Rs*EbNolin*states.bitspersymbol);
noise = sqrt(variance)*randn(1,length(tx));
SNRdB = 10*log10(var(tx)/var(noise));
rx = tx + noise;
f = fopen("fsk_demod.raw","wb"); tx = fwrite(f, rx, "short"); fclose(f);
% Demodulate using C modem and C de-framer/LDPC decoder
system("../build_linux/src/fsk_demod 2XS 8 9600 1200 fsk_demod.raw - | ../src/drs232_ldpc - dummy_out.bin");
end
% Plots uncoded BER curves for two different SSTV simulations. Used
% to compare results with different processing steps as we build up C
% command line. BER curves are powerful ways to confirm system is
% operating as expected, several bugs were found using this system.
if demo == 10
sim_in.frames = 10;
EbNodBvec = 7:10;
sim_in.demod_type = 3;
ber_test1 = [];
for i = 1:length(EbNodBvec)
[n_uncoded_errs n_uncoded_bits] = run_sstv_sim(sim_in, EbNodBvec(i));
ber_test1(i) = n_uncoded_errs/n_uncoded_bits;
end
sim_in.demod_type = 4;
ber_c = [];
for i = 1:length(EbNodBvec)
[n_uncoded_errs n_uncoded_bits] = run_sstv_sim(sim_in, EbNodBvec(i));
ber_test2(i) = n_uncoded_errs/n_uncoded_bits;
end
figure(1);
clf;
semilogy(EbNodBvec, ber_test1, '+-;first test;')
grid;
xlabel('Eb/No (dB)')
ylabel('BER')
hold on;
semilogy(EbNodBvec, ber_test2, 'g+-;second test;')
legend("boxoff");
hold off;
end
% Measure PER of complete coded and uncoded system
if demo == 11
sim_in.frames = 10;
EbNodB = 9;
sim_in.demod_type = 7;
run_sstv_sim(sim_in, EbNodB);
end
% Compare parity bits from an off-air stream of demodulated bits
% Use something like:
% cat ~/Desktop/923096fs_wenet.iq | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin CU8
% (note not soft dec mode)
if demo == 12
f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f);
compare_parity_bits(rx_bit_stream);
end

648
octave/tfsk_2400a.m
View File

@ -1,648 +0,0 @@
% tfsk.m
% Author: Brady O'Brien 8 January 2016
% Copyright 2016 David Rowe
%
% All rights reserved.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License version 2.1, as
% published by the Free Software Foundation. This program is
% distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or
% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program; if not, see <http://www.gnu.org/licenses/>.
% Octave script to check c port of fsk_horus against the fsk_horus.m
%
% [X] - Functions to wrap around fsk_mod and fsk_demod executables
% [X] - fsk_mod
% [X] - fsk_demod
% [X] - Functions to wrap around octave and c implementations, pass
% same dataset, compare outputs, and give clear go/no-go
% [X] - fsk_mod_test
% [X] - fsk_demod_test
% [X] - Port of run_sim and EbNodB curve test battery
% [X] - Extract and compare more parameters from demod
% [X] - Run some tests in parallel
#{
FSK Modem automated test instructions:
1. Use cmake to build in debug mode to ensure unittest/tfsk is built:
$ cd ~/codec2
$ rm -Rf build_linux && mkdir build_linux
$ cd build_linux
$ cmake -DCMAKE_BUILD_TYPE=Debug ..
$ make
2 - Change tfsk_location below if required
3 - Ensure Octave packages signal and parallel are installed
4 - Start Octave and run tfsk_2400a.m. It will perform all tests automatically
#}
%tfsk executable path/file
global tfsk_location = '../build_linux/unittest/tfsk';
%Set to 1 for verbose printouts
global print_verbose = 0;
fsk_horus_as_a_lib = 1; % make sure calls to test functions at bottom are disabled
%fsk_horus_2fsk;
fsk_horus
pkg load signal;
pkg load parallel;
graphics_toolkit('gnuplot');
global mod_pass_fail_maxdiff = 1e-3/5000;
function mod = fsk_mod_c(Fs,Rs,f1,fsp,bits,M)
global tfsk_location;
%command to be run by system to launch the modulator
command = sprintf('%s MX %d %d %d %d %d fsk_mod_ut_bitvec fsk_mod_ut_modvec fsk_mod_ut_log.txt',tfsk_location,M,f1,fsp,Fs,Rs);
%save input bits into a file
bitvecfile = fopen('fsk_mod_ut_bitvec','wb+');
fwrite(bitvecfile,bits,'uint8');
fclose(bitvecfile);
%run the modulator
system(command);
modvecfile = fopen('fsk_mod_ut_modvec','rb');
mod = fread(modvecfile,'single');
fclose(modvecfile);
endfunction
%Compare 2 vectors, fail if they are not close enough
function pass = vcompare(vc,voct,vname,tname,tol,pnum)
global print_verbose;
%Get delta of vectors
dvec = abs(abs(vc)-abs(voct));
%Normalize difference
dvec = dvec ./ abs(max(abs(voct))+1e-8);
maxdvec = abs(max(dvec));
pass = maxdvec<tol;
if print_verbose == 1
printf(' Comparing vectors %s in test %s. Diff is %f\n',vname,tname,maxdvec);
end
if pass == 0
printf('\n*** vcompare failed %s in test %s. Diff: %f Tol: %f\n\n',vname,tname,maxdvec,tol);
titlestr = sprintf('Diff between C and Octave of %s for %s',vname,tname)
figure(10+pnum*2)
plot(abs(dvec))
title(titlestr)
figure(11+pnum*2)
plot((1:length(vc)),abs(vc),(1:length(voct)),abs(voct))
end
endfunction
function test_stats = fsk_demod_xt(Fs,Rs,f1,fsp,mod,tname,M=2)
global tfsk_location;
global print_verbose;
%Name of executable containing the modulator
fsk_demod_ex_file = '../build/unittest/tfsk';
modvecfilename = sprintf('fsk_demod_ut_modvec_%d',getpid());
bitvecfilename = sprintf('fsk_demod_ut_bitvec_%d',getpid());
tvecfilename = sprintf('fsk_demod_ut_tracevec_%d.txt',getpid());
%command to be run by system to launch the demod
command = sprintf('%s DX %d %d %d %d %d %s %s %s',tfsk_location,M,f1,fsp,Fs,Rs,modvecfilename,bitvecfilename,tvecfilename);
%save modulated input into a file
modvecfile = fopen(modvecfilename,'wb+');
fwrite(modvecfile,mod,'single');
fclose(modvecfile);
%run the modulator
system(command);
bitvecfile = fopen(bitvecfilename,'rb');
bits = fread(bitvecfile,'uint8');
fclose(bitvecfile);
bits = bits!=0;
%Load test vec dump
load(tvecfilename);
%Clean up files
delete(bitvecfilename);
delete(modvecfilename);
delete(tvecfilename);
o_f1_dc = [];
o_f2_dc = [];
o_f3_dc = [];
o_f4_dc = [];
o_f1_int = [];
o_f2_int = [];
o_f3_int = [];
o_f4_int = [];
o_f1 = [];
o_f2 = [];
o_f3 = [];
o_f4 = [];
o_EbNodB = [];
o_ppm = [];
o_Sf = [];
o_fest = [];
o_rx_timing = [];
o_norm_rx_timing = [];
o_nin = [];
%Run octave demod, dump some test vectors
states = fsk_horus_init_hbr(Fs,10,Rs,M);
Ts = states.Ts;
P = states.P;
states.ftx(1) = f1;
states.ftx(2) = f1+fsp;
states.ftx(3) = f1+fsp*2;
states.ftx(4) = f1+fsp*3;
states.dA = 1;
states.dF = 0;
modin = mod;
obits = [];
while length(modin)>=states.nin
ninold = states.nin;
states = est_freq(states, modin(1:states.nin), states.M);
[bitbuf,states] = fsk_horus_demod(states, modin(1:states.nin));
modin=modin(ninold+1:length(modin));
obits = [obits bitbuf];
%Save other parameters
o_f1_dc = [o_f1_dc states.f_dc(1,1:states.Nmem-Ts/P)];
o_f2_dc = [o_f2_dc states.f_dc(2,1:states.Nmem-Ts/P)];
o_f1_int = [o_f1_int states.f_int(1,:)];
o_f2_int = [o_f2_int states.f_int(2,:)];
o_EbNodB = [o_EbNodB states.EbNodB];
o_ppm = [o_ppm states.ppm];
o_rx_timing = [o_rx_timing states.rx_timing];
o_norm_rx_timing = [o_norm_rx_timing states.norm_rx_timing];
o_Sf = [o_Sf states.Sf'];
o_f1 = [o_f1 states.f(1)];
o_f2 = [o_f1 states.f(2)];
o_fest = [o_fest states.f];
o_nin = [o_nin states.nin];
if M==4
o_f3_dc = [o_f3_dc states.f_dc(3,1:states.Nmem-Ts/P)];
o_f4_dc = [o_f4_dc states.f_dc(4,1:states.Nmem-Ts/P)];
o_f3_int = [o_f3_int states.f_int(3,:)];
o_f4_int = [o_f4_int states.f_int(4,:)];
o_f3 = [o_f1 states.f(3)];
o_f4 = [o_f1 states.f(4)];
end
end
%close all
pass = 1;
pass = vcompare(o_Sf, t_fft_est(1:length(o_Sf)),'fft est',tname,1,1) && pass;
pass = vcompare(o_fest, t_f_est,'f est',tname,1,2) && pass;
pass = vcompare(o_rx_timing, t_rx_timing,'rx timing',tname,1,3) && pass;
if M==4
pass = vcompare(o_f3_dc, t_f3_dc, 'f3 dc', tname,1,4) && pass;
pass = vcompare(o_f4_dc, t_f4_dc, 'f4 dc', tname,1,5) && pass;
pass = vcompare(o_f3_int, t_f3_int, 'f3 int', tname,1,6) && pass;
pass = vcompare(o_f4_int, t_f4_int, 'f4 int', tname,1,7) && pass;
end
pass = vcompare(o_f1_dc, t_f1_dc, 'f1 dc', tname,1,8) && pass;
pass = vcompare(o_f2_dc, t_f2_dc, 'f2 dc', tname,1,9) && pass;
pass = vcompare(o_f2_int, t_f2_int, 'f2 int', tname,1,10) && pass;
pass = vcompare(o_f1_int, t_f1_int, 'f1 int', tname,1,11) && pass;
pass = vcompare(o_ppm , t_ppm, 'ppm', tname,1,12) && pass;
pass = vcompare(o_EbNodB, t_EbNodB,'EbNodB', tname,1,13) && pass;
pass = vcompare(o_nin, t_nin, 'nin', tname,1,14) && pass;
pass = vcompare(o_norm_rx_timing, t_norm_rx_timing,'norm rx timing',tname,1,15) && pass;
diffpass = sum(xor(obits,bits'))<5;
diffbits = sum(xor(obits,bits'));
if print_verbose == 1
printf('%d bit diff in test %s\n',diffbits,tname);
end
if diffpass==0
printf('\n***bitcompare test failed test %s diff %d\n\n',tname,sum(xor(obits,bits')))
figure(15)
plot(xor(obits,bits'))
title(sprintf('Bitcompare failure test %s',tname))
end
pass = pass && diffpass;
assert(pass);
test_stats.pass = pass;
test_stats.diff = sum(xor(obits,bits'));
test_stats.cbits = bits';
test_stats.obits = obits;
endfunction
function [dmod,cmod,omod,pass] = fsk_mod_test(Fs,Rs,f1,fsp,bits,tname,M=2)
global mod_pass_fail_maxdiff;
%Run the C modulator
cmod = fsk_mod_c(Fs,Rs,f1,fsp,bits,M);
%Set up and run the octave modulator
states.M = M;
states = fsk_horus_init_hbr(Fs,10,Rs,M);
states.ftx(1) = f1;
states.ftx(2) = f1+fsp;
if states.M == 4
states.ftx(3) = f1+fsp*2;
states.ftx(4) = f1+fsp*3;
end
states.dA = [1 1 1 1];
states.dF = 0;
omod = fsk_horus_mod(states,bits);
dmod = cmod-omod;
pass = max(dmod)<(mod_pass_fail_maxdiff*length(dmod));
if !pass
printf('Mod failed test %s!\n',tname);
end
endfunction
% Random bit modulator test
% Pass random bits through the modulators and compare
function pass = test_mod_2400a_randbits
rand('state',1);
randn('state',1);
bits = rand(1,96000)>.5;
[dmod,cmod,omod,pass] = fsk_mod_test(48000,1200,1200,1200,bits,"mod 2400a randbits",4);
if(!pass)
figure(1)
plot(dmod)
title("Difference between octave and C mod impl");
end
endfunction
% A big ol' channel impairment tester
% Shamlessly taken from fsk_horus
% This throws some channel imparment or another at the C and octave modem so they
% may be compared.
function stats = tfsk_run_sim(test_frame_mode,EbNodB,timing_offset,fading,df,dA,M=2)
global print_verbose;
frames = 190;
%EbNodB = 10;
%timing_offset = 2.0; % see resample() for clock offset below
%fading = 0; % modulates tx power at 2Hz with 20dB fade depth,
% to simulate balloon rotating at end of mission
%df = 0; % tx tone freq drift in Hz/s
%dA = 1; % amplitude imbalance of tones (note this affects Eb so not a gd idea)
more off
rand('state',10);
randn('state',10);
% ----------------------------------------------------------------------
% sm2000 config ------------------------
%states = fsk_horus_init(96000, 1200);
%states.f1_tx = 4000;
%states.f2_tx = 5200;
if test_frame_mode == 2
% 2400A config
states = fsk_horus_init_hbr(48000,10, 1200, M);
states.f1_tx = 1200;
states.f2_tx = 2400;
states.f3_tx = 3600;
states.f4_tx = 4800;
states.ftx(1) = 1200;
states.ftx(2) = 2400;
states.ftx(3) = 3600;
states.ftx(4) = 4800;
end
if test_frame_mode == 4
% horus rtty config ---------------------
states = fsk_horus_init_hbr(48000,10, 1200, M);
states.f1_tx = 1200;
states.f2_tx = 2400;
states.f3_tx = 3600;
states.f4_tx = 4800;
states.ftx(1) = 1200;
states.ftx(2) = 2400;
states.ftx(3) = 3600;
states.ftx(4) = 4800;
states.tx_bits_file = "horus_tx_bits_rtty.txt"; % Octave file of bits we FSK modulate
end
if test_frame_mode == 5
% horus binary config ---------------------
states = fsk_horus_init_hbr(48000,10, 1200, M);
states.f1_tx = 1200;
states.f2_tx = 2400;
states.f3_tx = 3600;
states.f4_tx = 4800;
states.ftx(1) = 1200;
states.ftx(2) = 2400;
states.ftx(3) = 3600;
states.ftx(4) = 4800;
%%%states.tx_bits_file = "horus_tx_bits_binary.txt"; % Octave file of bits we FSK modulate
states.tx_bits_file = "horus_payload_rtty.txt";
end
% ----------------------------------------------------------------------
states.verbose = 0;%x1;
N = states.N;
P = states.P;
Rs = states.Rs;
nsym = states.nsym;
Fs = states.Fs;
states.df = df;
states.dA = [dA dA dA dA];
states.M = M;
EbNo = 10^(EbNodB/10);
variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol);
% set up tx signal with payload bits based on test mode
if test_frame_mode == 1
% test frame of bits, which we repeat for convenience when BER testing
test_frame = round(rand(1, states.nsym));
tx_bits = [];
for i=1:frames+1
tx_bits = [tx_bits test_frame];
end
end
if test_frame_mode == 2
% random bits, just to make sure sync algs work on random data
tx_bits = round(rand(1, states.nbit*(frames+1)));
end
if test_frame_mode == 3
% ...10101... sequence
tx_bits = zeros(1, states.nsym*(frames+1));
tx_bits(1:2:length(tx_bits)) = 1;
end
if (test_frame_mode == 4) || (test_frame_mode == 5)
% load up a horus msg from disk and modulate that
test_frame = load(states.tx_bits_file);
ltf = length(test_frame);
ntest_frames = ceil((frames+1)*nsym/ltf);
tx_bits = [];
for i=1:ntest_frames
tx_bits = [tx_bits test_frame];
end
end
f1 = states.f1_tx;
fsp = states.f2_tx-f1;
states.dA = [dA dA dA dA];
states.ftx(1) = f1;
states.ftx(2) = f1+fsp;
if states.M == 4
states.ftx(3) = f1+fsp*2;
states.ftx(4) = f1+fsp*3;
end
tx = fsk_horus_mod(states, tx_bits);
if timing_offset
tx = resample(tx, 1000, 1001); % simulated 1000ppm sample clock offset
end
if fading
ltx = length(tx);
tx = tx .* (1.1 + cos(2*pi*2*(0:ltx-1)/Fs))'; % min amplitude 0.1, -20dB fade, max 3dB
end
noise = sqrt(variance)*randn(length(tx),1);
rx = tx + noise;
test_name = sprintf("tfsk EbNodB:%d frames:%d timing_offset:%d fading:%d df:%d",EbNodB,frames,timing_offset,fading,df);
tstats = fsk_demod_xt(Fs,Rs,states.f1_tx,fsp,rx,test_name,M);
pass = tstats.pass;
obits = tstats.obits;
cbits = tstats.cbits;
stats.name = test_name;
if tstats.pass
printf("Test %s passed\n",test_name);
else
printf("Test %s failed\n",test_name);
end
% Figure out BER of octave and C modems
bitcnt = length(tx_bits);
rx_bits = obits;
ber = 1;
ox = 1;
for offset = (1:100)
nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset)));
bern = nerr/(bitcnt-offset);
if(bern < ber)
ox = offset;
best_nerr = nerr;
end
ber = min([ber bern]);
end
offset = ox;
bero = ber;
ber = 1;
rx_bits = cbits;
ox = 1;
for offset = (1:100)
nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset)));
bern = nerr/(bitcnt-offset);
if(bern < ber)
ox = offset;
best_nerr = nerr;
end
ber = min([ber bern]);
end
offset = ox;
berc = ber;
stats.berc = berc;
stats.bero = bero;
% coherent BER theory calculation
if print_verbose == 1
printf("C BER: %f Oct BER: %f Test %s\n",berc,bero,test_name);
end
stats.thrcoh = .5*(M-1)*erfc(sqrt( (log2(M)/2) * EbNo ));
% non-coherent BER theory calculation
% It was complicated, so I broke it up
ms = M;
ns = (1:ms-1);
as = (-1).^(ns+1);
bs = (as./(ns+1));
cs = ((ms-1)./ns);
ds = ns.*log2(ms);
es = ns+1;
fs = exp( -(ds./es)*EbNo );
thrncoh = ((ms/2)/(ms-1)) * sum(bs.*((ms-1)./ns).*exp( -(ds./es)*EbNo ));
stats.thrncoh = thrncoh;
stats.pass = pass;
endfunction
function pass = ebno_battery_test(timing_offset,fading,df,dA,M)
%Range of EbNodB over which to test
ebnodbrange = fliplr(5:2:13);
ebnodbs = length(ebnodbrange);
mode = 2;
%Replication of other parameters for parcellfun
modev = repmat(mode,1,ebnodbs);
timingv = repmat(timing_offset,1,ebnodbs);
fadingv = repmat(fading,1,ebnodbs);
dfv = repmat(df,1,ebnodbs);
dav = repmat(dA,1,ebnodbs);
mv = repmat(M,1,ebnodbs);
statv = pararrayfun(floor(1.25*nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,mv);
%statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,mv);
passv = zeros(1,length(statv));
for ii=(1:length(statv))
passv(ii)=statv(ii).pass;
if statv(ii).pass
printf("Test %s passed\n",statv(ii).name);
else
printf("Test %s failed\n",statv(ii).name);
end
end
%All pass flags are '1'
pass = sum(passv)>=length(passv);
%and no tests died
pass = pass && length(passv)==ebnodbs;
passv;
assert(pass)
endfunction
%Test with and without sample clock offset
function pass = test_timing_var(df,dA,M)
pass = ebno_battery_test(1,0,df,dA,M)
assert(pass)
pass = pass && ebno_battery_test(0,0,df,dA,M)
assert(pass)
endfunction
%Test with and without 1 Hz/S freq drift
function pass = test_drift_var(M)
pass = test_timing_var(1,1,M)
assert(pass)
pass = pass && test_timing_var(0,1,M)
assert(pass)
endfunction
function pass = test_fsk_battery()
pass = 1;
pass = pass && test_mod_2400a_randbits;
assert(pass)
pass = pass && test_drift_var(4);
assert(pass)
if pass
printf("***** All tests passed! *****\n");
end
endfunction
function plot_fsk_bers(M=2)
%Range of EbNodB over which to plot
ebnodbrange = (4:13);
berc = ones(1,length(ebnodbrange));
bero = ones(1,length(ebnodbrange));
berinc = ones(1,length(ebnodbrange));
beric = ones(1,length(ebnodbrange));
ebnodbs = length(ebnodbrange)
mode = 2;
%Replication of other parameters for parcellfun
modev = repmat(mode,1,ebnodbs);
timingv = repmat(1,1,ebnodbs);
fadingv = repmat(0,1,ebnodbs);
dfv = repmat(1,1,ebnodbs);
dav = repmat(1,1,ebnodbs);
Mv = repmat(M,1,ebnodbs);
statv = pararrayfun(floor(nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv);
%statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv);
for ii = (1:length(statv))
stat = statv(ii);
berc(ii)=stat.berc;
bero(ii)=stat.bero;
berinc(ii)=stat.thrncoh;
beric(ii) = stat.thrcoh;
end
clf;
figure(M)
semilogy(ebnodbrange, berinc,sprintf('r;%dFSK non-coherent theory;',M))
hold on;
semilogy(ebnodbrange, beric ,sprintf('g;%dFSK coherent theory;',M))
semilogy(ebnodbrange, bero ,sprintf('b;Octave fsk horus %dFSK Demod;',M))
semilogy(ebnodbrange, berc,sprintf('+;C fsk horus %dFSK Demod;',M))
hold off;
grid("minor");
axis([min(ebnodbrange) max(ebnodbrange) 1E-5 1])
legend("boxoff");
xlabel("Eb/No (dB)");
ylabel("Bit Error Rate (BER)")
endfunction
xpass = test_fsk_battery
%plot_fsk_bers(2)
plot_fsk_bers(4)
if xpass
printf("***** All tests passed! *****\n");
else
printf("***** Some test failed! Look back through output to find failed test *****\n");
end

10247
octave/vq
View File

File diff suppressed because it is too large Load Diff

210
octave/vq_binary_switch.m
View File

@ -1,210 +0,0 @@
% vq_binary_switch.m
% David Rowe Sep 2021
%
% Experiments in making VQs robust to bit errors, this is an Octave
% implementation of [1].
%
% [1] Pseudo Gray Coding, Zeger & Gersho 1990
1;
% returns indexes of hamming distance 1 neighbours
function index_neighbours = distance_one_neighbours(N,k)
log2N = log2(N);
index_neighbours = [];
for b=0:log2N-1
index_neighbour = bitxor(k-1,2.^b) + 1;
index_neighbours = [index_neighbours index_neighbour];
end
end
% equation (33) of [1], for hamming distance 1
function c = cost_of_distance_one(vq, prob, k, verbose=0)
[N K] = size(vq);
log2N = log2(N);
c = 0;
for b=0:log2N-1
index_neighbour = bitxor(k-1,2.^b) + 1;
diff = vq(k,:) - vq(index_neighbour, :);
dist = sum(diff*diff');
c += prob(k)*dist;
if verbose
printf("k: %d b: %d index_neighbour: %d dist: %f prob: %f c: %f \n", k, b, index_neighbour, dist, prob(k), c);
end
end
endfunction
% equation (39) of [1]
function d = distortion_of_current_mapping(vq, prob, verbose=0)
[N K] = size(vq);
d = 0;
for k=1:N
c = cost_of_distance_one(vq, prob, k);
d += c;
if verbose
printf("k: %2d c: %f d: %f\n", k, c, d);
end
end
endfunction
function [vq distortion] = binary_switching(vq, prob, max_iteration, fast_en=1)
[N K] = size(vq);
iteration = 0;
i = 1;
finished = 0;
switches = 0;
distortion0 = distortion_of_current_mapping(vq, prob)
while !finished
% generate a list A(i) of which vectors have the largest cost of bit errors
c = zeros(1,N);
for k=1:N
c(k) = cost_of_distance_one(vq, prob, k);
end
[tmp A] = sort(c,"descend");
% Try switching each vector with A(i)
best_delta = 0;
for j=2:N
% we can't switch with ourself
if j != A(i)
if fast_en
delta = -cost_of_distance_one(vq, prob, A(i)) - cost_of_distance_one(vq, prob, j);
n1 = [distance_one_neighbours(N,A(i)) distance_one_neighbours(N,j)];
n1(n1 == A(i)) = [];
n1(n1 == j) = [];
for l=1:length(n1)
delta -= cost_of_distance_one(vq, prob, n1(l));
end
else
distortion1 = distortion_of_current_mapping(vq, prob);
end
% switch vq entries A(i) and j
tmp = vq(A(i),:);
vq(A(i),:) = vq(j,:);
vq(j,:) = tmp;
if fast_en
delta += cost_of_distance_one(vq, prob, A(i)) + cost_of_distance_one(vq, prob, j);
for l=1:length(n1)
delta += cost_of_distance_one(vq, prob, n1(l));
end
else
distortion2 = distortion_of_current_mapping(vq, prob);
delta = distortion2 - distortion1;
end
if delta < 0
if abs(delta) > best_delta
best_delta = abs(delta);
best_j = j;
end
end
% unswitch
tmp = vq(A(i),:);
vq(A(i),:) = vq(j,:);
vq(j,:) = tmp;
end
end % next j
% printf("best_delta: %f best_j: %d\n", best_delta, best_j);
if best_delta == 0
% Hmm, no improvement, lets try the next vector in the sorted cost list
if i == N
finished = 1;
else
i++;
end
else
% OK keep the switch that minimised the distortion
tmp = vq(A(i),:);
vq(A(i),:) = vq(best_j,:);
vq(best_j,:) = tmp;
switches++;
% set up for next iteration
iteration++;
distortion = distortion_of_current_mapping(vq, prob);
printf("it: %3d dist: %f %3.2f i: %3d sw: %3d\n", iteration, distortion,
distortion/distortion0, i, switches);
if iteration >= max_iteration, finished = 1, end
i = 1;
end
end
endfunction
% return indexes of hamming distance one vectors
function ind = neighbour_indexes(vq, k)
[N K] = size(vq);
log2N = log2(N);
ind = [];
for b=0:log2N-1
index_neighbour = bitxor(k-1,2.^b) + 1;
ind = [ind index_neighbour];
end
endfunction
function test_binary_switch
vq1 = [1 1; -1 1; -1 -1; 1 -1];
%f=fopen("vq1.f32","wb"); fwrite(f, vq1, 'float32'); fclose(f);
[vq2 distortion] = binary_switching(vq1, ones(1,4), 10);
% algorithm should put hamming distance 1 neighbours in adjacent quadrants
distance_to_closest_neighbours = 2;
% there are two hamming distance 1 neighbours
target_distortion = 2^2*distance_to_closest_neighbours*length(vq1);
assert(target_distortion == distortion);
printf("test_binary_switch OK!\n");
endfunction
function test_fast
N=16; % Number of VQ codebook vectors
K=2; % Vector length
Ntrain=10000;
training_data = randn(Ntrain,K);
[idx vq1] = kmeans(training_data, N);
f=fopen("vq1.f32","wb");
for r=1:rows(vq1)
fwrite(f,vq1(r,:),"float32");
end
fclose(f);
[vq2 distortion] = binary_switching(vq1, [1 ones(1,N-1)], 1000, fast_en = 0);
[vq3 distortion] = binary_switching(vq1, [1 ones(1,N-1)], 1000, fast_en = 1);
assert(vq2 == vq3);
printf("test_fast OK!\n");
endfunction
function demo
N=16; % Number of VQ codebook vectors
K=2; % Vector length
Ntrain=10000;
training_data = randn(Ntrain,K);
[idx vq1] = kmeans(training_data, N);
[vq2 distortion] = binary_switching(vq1, [1 ones(1,N-1)], 1000, 1);
figure(1); clf; plot(training_data(:,1), training_data(:,2),'+');
hold on;
plot(vq1(:,1), vq1(:,2),'og','linewidth', 2);
plot(vq2(:,1), vq2(:,2),'or','linewidth', 2);
% plot hamming distance 1 neighbours
k = 1;
ind = neighbour_indexes(vq2, k);
for i=1:length(ind)
plot([vq2(k,1) vq2(ind(i),1)],[vq2(k,2) vq2(ind(i),2)],'r-','linewidth', 2);
end
hold off;
endfunction
pkg load statistics
%test_binary_switch;
test_fast;
%demo

1615
octave/vq_nm_xx_11_40_hpf150
View File

File diff suppressed because it is too large Load Diff