1% mfsk.m
2% David Rowe Nov 2015
3
4% Simulation to test m=2 and m=4 FSK demod
5
6
71;
8
9function sim_out = fsk_ber_test(sim_in)
10  Fs        = 96000;
11  M         = sim_in.M;
12  Rs        = sim_in.Rs;
13  Ts        = Fs/Rs;
14  verbose   = sim_in.verbose;
15
16  nbits     = sim_in.nbits;
17  nsym      = sim_in.nbits*2/M;
18  nsam      = nsym*Ts;
19  EsNodB    = sim_in.EbNodB + 10*log10(M/2);
20
21  % printf("M: %d nbits: %d nsym: %d\n", M, nbits, nsym);
22
23  if M == 2
24    f(1) = -Rs/2;
25    f(2) =  Rs/2;
26  end
27  if M == 4
28    f(1) = -3*Rs/2;
29    f(2) = -Rs/2;
30    f(3) =  Rs/2;
31    f(4) =  3*Rs/2;
32  end
33
34  % simulate over a range of Eb/No values
35
36  for ne = 1:length(EsNodB)
37    Nerrs = Terrs = Tbits = 0;
38
39    aEsNodB = EsNodB(ne);
40    EsNo = 10^(aEsNodB/10);
41    variance = Fs/(Rs*EsNo);
42
43    % Modulator -------------------------------
44
45    tx_bits = round(rand(1, nbits));
46    tx = zeros(1,nsam);
47    tx_phase = 0;
48
49    for i=1:nsym
50      if M == 2
51        tone = tx_bits(i) + 1;
52      else
53        tone = (tx_bits(2*(i-1)+1:2*i) * [2 1]') + 1;
54      end
55
56      tx_phase_vec = tx_phase + (1:Ts)*2*pi*f(tone)/Fs;
57      tx((i-1)*Ts+1:i*Ts) = exp(j*tx_phase_vec);
58      tx_phase = tx_phase_vec(Ts) - floor(tx_phase_vec(Ts)/(2*pi))*2*pi;
59    end
60
61    % Channel ---------------------------------
62
63    % We use complex (single sided) channel simulation, as it's convenient
64    % for the FM simulation.
65
66    noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam));
67    rx    = tx + noise;
68    if verbose > 1
69      printf("EbNo: %f Eb: %f var No: %f EbNo (meas): %f\n",
70      EbNo, var(tx)*Ts/Fs, var(noise)/Fs, (var(tx)*Ts/Fs)/(var(noise)/Fs));
71    end
72
73    % Demodulator -----------------------------
74
75    % non-coherent FSK demod
76
77    rx_bb = rx;
78    dc = zeros(M,nsam);
79    for m=1:M
80      dc(m,:) = rx_bb .* exp(-j*(0:nsam-1)*2*pi*f(m)/Fs);
81    end
82
83    rx_bits = zeros(1, nsym);
84    for i=1:nsym
85      st = (i-1)*Ts+1;
86      en = st+Ts-1;
87      for m=1:M
88        int(m,i) = abs(sum(dc(m,st:en)));
89      end
90      if m == 2
91        rx_bits(i) = int(1,i) < int(2,i);
92      else
93        [max_amp tone] = max([int(1,i) int(2,i) int(3,i) int(4,i)]);
94        if tone == 1
95          rx_bits(2*(i-1)+1:2*i) = [0 0];
96        end
97        if tone == 2
98          rx_bits(2*(i-1)+1:2*i) = [0 1];
99        end
100        if tone == 3
101          rx_bits(2*(i-1)+1:2*i) = [1 0];
102        end
103        if tone == 4
104          rx_bits(2*(i-1)+1:2*i) = [1 1];
105        end
106     end
107    end
108
109    error_positions = xor(rx_bits, tx_bits);
110    Nerrs = sum(error_positions);
111    Terrs += Nerrs;
112    Tbits += length(error_positions);
113
114    TERvec(ne) = Terrs;
115    BERvec(ne) = Terrs/Tbits;
116
117    if verbose > 1
118      figure(2)
119      clf
120      Rx = 10*log10(abs(fft(rx)));
121      plot(Rx(1:Fs/2));
122      axis([1 Fs/2 0 50]);
123
124      figure(3)
125      clf;
126      subplot(211)
127      plot(real(rx_bb(1:Ts*20)))
128      subplot(212)
129      Rx_bb = 10*log10(abs(fft(rx_bb)));
130      plot(Rx_bb(1:3000));
131      axis([1 3000 0 50]);
132
133      figure(4);
134      subplot(211)
135      stem(abs(mark_int(1:100)));
136      subplot(212)
137      stem(abs(space_int(1:100)));
138    end
139
140    if verbose
141      printf("EbNo (db): %3.2f Terrs: %d BER: %4.3f \n", aEsNodB - 10*log10(M/2), Terrs, Terrs/Tbits);
142    end
143  end
144
145  sim_out.TERvec = TERvec;
146  sim_out.BERvec = BERvec;
147endfunction
148
149
150function run_fsk_curves
151  sim_in.M         = 2;
152  sim_in.Rs        = 1200;
153  sim_in.nbits     = 12000;
154  sim_in.EbNodB    = 0:2:20;
155  sim_in.verbose   = 1;
156
157  EbNo  = 10 .^ (sim_in.EbNodB/10);
158  fsk_theory.BERvec = 0.5*exp(-EbNo/2); % non-coherent BFSK demod
159  fsk2_sim = fsk_ber_test(sim_in);
160
161  sim_in.M = 4;
162  fsk4_sim = fsk_ber_test(sim_in);
163
164  % BER v Eb/No curves
165
166  figure(1);
167  clf;
168  semilogy(sim_in.EbNodB, fsk_theory.BERvec,'r;2FSK theory;')
169  hold on;
170  semilogy(sim_in.EbNodB, fsk2_sim.BERvec,'g;2FSK sim;')
171  semilogy(sim_in.EbNodB, fsk4_sim.BERvec,'b;4FSK sim;')
172  hold off;
173  grid("minor");
174  axis([min(sim_in.EbNodB) max(sim_in.EbNodB) 1E-4 1])
175  legend("boxoff");
176  xlabel("Eb/No (dB)");
177  ylabel("Bit Error Rate (BER)")
178
179end
180
181
182function run_fsk_single
183  sim_in.M         = 4;
184  sim_in.Rs        = 1200;
185  sim_in.nbits     = 5000;
186  sim_in.EbNodB    = 8;
187  sim_in.verbose   = 1;
188
189  fsk_sim = fsk_ber_test(sim_in);
190endfunction
191
192
193rand('state',1);
194randn('state',1);
195graphics_toolkit ("gnuplot");
196
197run_fsk_curves
198%run_fsk_single
199
200