1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
|
% fsk_lib_ldpc_demo.m
%
% LDPC coded 4FSK modem demo, demonstrating soft dec using CML library functions
fsk_lib;
ldpc;
% set up waveform
function [states M] = modem_init(Rs,Fs,df)
M = 4;
states = fsk_init(Fs,Rs,M,P=8,nsym=100);
states.tx_real = 0;
states.tx_tone_separation = 250;
states.ftx = -2.5*states.tx_tone_separation + states.tx_tone_separation*(1:M);
states.fest_fmin = -Fs/2;
states.fest_fmax = +Fs/2;
states.fest_min_spacing = Rs/2;
states.df = df;
states.ber_valid_thresh = 0.1;
states.ber_invalid_thresh = 0.2;
end
% Run a complete modem (freq and timing estimators running) at a
% single Eb/No point. At low Eb/No the estimators occasionally fall
% over so we get complete junk, we consider that case a packet error
% and exclude it from the BER estimation.
function [states uber cber cper] = modem_run_test(HRA, EbNodB = 10, num_frames=10, Fs=8000, Rs=100, df=0, plots=0)
rand('seed',1); randn('seed',1);
[states M] = modem_init(Rs, Fs, df);
N = states.N;
if plots; states.verbose = 0x4; end
% set up LDPC code
Hsize=size(HRA);
Krate = (Hsize(2)-Hsize(1))/Hsize(2); states.rate = Krate;
code_param = ldpc_init_user(HRA, modulation='FSK', mod_order=states.M, mapping='gray');
states.coden = code_param.coded_bits_per_frame;
states.codek = code_param.data_bits_per_frame;
% set up AWGN noise
EcNodB = EbNodB + 10*log10(Krate);
EcNo = 10^(EcNodB/10);
variance = states.Fs/(states.Rs*EcNo*states.bitspersymbol);
data_bits = round(rand(1,code_param.data_bits_per_frame)); tx_bits = [];
for f=1:num_frames
codeword_bits = LdpcEncode(data_bits, code_param.H_rows, code_param.P_matrix);
tx_bits = [tx_bits codeword_bits];
end
% modulator and AWGN channel
tx = fsk_mod(states, tx_bits);
noise = sqrt(variance/2)*randn(length(tx),1) + j*sqrt(variance/2)*randn(length(tx),1);
rx = tx + noise;
% freq estimator and demod
run_frames = floor(length(rx)/N)-1;
st = 1; f_log = []; rx_bits = []; rx_filt = [];
for f=1:run_frames
% extract nin samples from input stream
nin = states.nin;
en = st + states.nin - 1;
% due to nin variations it's possible to overrun buffer
if en < length(rx)
sf = rx(st:en);
states = est_freq(states, sf, states.M); states.f = states.f2;
[arx_bits states] = fsk_demod(states, sf);
rx_bits = [rx_bits arx_bits];
rx_filt = [rx_filt abs(states.f_int_resample)];
f_log = [f_log; states.f];
st += nin;
end
end
% count bit errors in test frames
num_frames=floor(length(rx_bits)/code_param.coded_bits_per_frame);
log_nerrs = []; num_frames_rx = 0; Tbits = Terrs = Tperr = Tpackets = 0;
uber = cber = 0.5; cper = 1;
for f=1:num_frames-1
st = (f-1)*code_param.coded_bits_per_frame + 1; en = (f+1)*code_param.coded_bits_per_frame;
states = ber_counter(states, codeword_bits, rx_bits(st:en));
log_nerrs = [log_nerrs states.nerr];
if states.ber_state num_frames_rx++; end
% Using sync provided by ber_counter() state machine for LDPC frame alignment
if states.ber_state
st_bit = (f-1)*code_param.coded_bits_per_frame + states.coarse_offset;
st_symbol = (st_bit-1)/states.bitspersymbol + 1;
en_symbol = st_symbol + code_param.coded_bits_per_frame/states.bitspersymbol - 1;
%printf("coded_bits: %d bps: %d st_bit: %d st_symbol: %d en_symbol: %d\n",
%code_param.coded_bits_per_frame, states.bitspersymbol, st_bit, st_symbol, en_symbol);
% map FSK filter outputs to LLRs, then LDPC decode (see also fsk_cml_sam.m)
symL = DemodFSK(1/states.v_est*rx_filt(:,st_symbol:en_symbol), states.SNRest, 1);
llr = -Somap(symL);
[x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, max_iterations=100, decoder_type=0, 1, 1);
Niters = sum(PCcnt~=0);
detected_data = x_hat(Niters,:);
Nerrs = sum(xor(data_bits, detected_data(1:code_param.data_bits_per_frame)));
Terrs += Nerrs;
Tbits += code_param.data_bits_per_frame;
if Nerrs Tperr++; end
Tpackets++;
end
end
if states.Terrs
printf("Fs: %d Rs: %d df % 3.2f EbNo: %4.2f ftx: %3d frx: %3d\n",Fs, Rs, df, EbNodB, num_frames, num_frames_rx);
uber = states.Terrs/states.Tbits; cber = Terrs/Tbits; cper = Tperr/Tpackets;
printf(" Uncoded: nbits: %6d nerrs: %6d ber: %4.3f\n", states.Tbits, states.Terrs, uber);
printf(" Coded..: nbits: %6d nerrs: %6d ber: %4.3f\n", Tbits, Terrs, cber);
printf(" Coded..: npckt: %6d perrs: %6d per: %4.3f\n", Tpackets, Tperr, cper);
end
if plots
figure(1); clf;
ideal=ones(length(f_log),1)*states.ftx;
plot((1:length(f_log)),ideal(:,1),'bk;ideal;')
hold on; plot((1:length(f_log)),ideal(:,2:states.M),'bk'); hold off;
hold on;
plot(f_log(:,1), 'linewidth', 2, 'b;peak;');
plot(f_log(:,2:states.M), 'linewidth', 2, 'b');
hold off;
xlabel('Time (frames)'); ylabel('Frequency (Hz)');
figure(2); clf; plot(log_nerrs); title('Errors per frame');
end
end
function freq_run_curve_peak_mask(HRA, num_frames=100)
EbNodB = 4:10;
m4fsk_ber_theory = [0.23 0.18 0.14 0.09772 0.06156 0.03395 0.01579 0.00591 0.00168 3.39E-4];
uber_log = []; cber_log = []; cper_log = [];
for ne = 1:length(EbNodB)
[states uber cber cper] = modem_run_test(HRA, EbNodB(ne), num_frames);
uber_log = [uber_log uber]; cber_log = [cber_log cber]; cper_log = [cper_log cper];
end
figure(1); clf;
EbNodB_raw = EbNodB+10*log10(states.rate)
semilogy(EbNodB_raw, m4fsk_ber_theory(round(EbNodB_raw+1)), 'linewidth', 2, 'bk+-;uber theory;');
grid; hold on;
semilogy(EbNodB_raw, uber_log+1E-12, 'linewidth', 2, '+-;uber;');
semilogy(EbNodB, cber_log+1E-12, 'linewidth', 2, 'r+-;cber;');
semilogy(EbNodB, cper_log+1E-12, 'linewidth', 2, 'c+-;cper;'); hold off;
xlabel('Eb/No (info bits, dB)'); ylabel('BER/PER'); axis([min(EbNodB_raw) max(EbNodB) 1E-4 1]);
title(sprintf("%dFSK rate %3.1f (%d,%d) Ncodewords=%d NCodewordBits=%d Fs=%d Rs=%d",
states.M, states.rate, states.coden, states.codek, num_frames, states.Tbits, states.Fs, states.Rs));
print("fsk_lib_ldpc.png", "-dpng")
end
% Choose simulation here ---------------------------------------------------
init_cml();
load H_256_512_4.mat; HRA=H;
more off;
% single point
[states uber cber cper] = modem_run_test(HRA, EbNodB=8);
if cber == 0
printf("PASS\n");
end
% curve
%freq_run_curve_peak_mask(HRA, 200)
|
