1% mancyfsk.m 2% David Rowe October 2015 3% 4% Manchester encoded 2FSK & 4FSK simulation. 5% 6% Attempt to design a FSK waveform that can pass through legacy FM 7% radios but still be optimally demodulated by SDRs. It doesn't have 8% to be optimally demodulated by legacy radios. Trick is getting it 9% to pass through 300-3000Hz audio filters in legacy radios. 10% 11% [X] code up modulator 12% [X] manchester two bit symbols 13% [X] plot spectrum 14% [X] demodulate using analog FM and ideal demods 15% [X] measure BER compared to ideal coherent FSK 16 171; 18 19fm; % analog FM library 20 21 22function states = legacyfsk_init(M,Rs) 23 Fs = states.Fs = 96000; 24 states.Rs = Rs; % symbol rate over channel 25 Ts = states.Ts = Fs/Rs; % symbol period in samples 26 states.M = M; % mFSK, either 2 or 4 27 bpsym = state.Rb = log2(M); % bits per symbol over channel 28 rate = states.rate = 0.5; % Manchester code rate 29 nbits = 100; 30 nbits = states.nbits = 100; % number of payload data symbols/frame 31 nbits2 = states.nbits2 = nbits/rate; % number of symbols/frame over channel after manchester encoding 32 nsym = states.nsym = nbits2/log2(M); % number of symbols per frame 33 nsam = states.nsam = nsym*Ts; 34 35 %printf(" Rs: %d M: %d bpsym: %d nbits: %d nbits2: %d nsym: %d nsam: %d\n", Rs, M, bpsym, nbits, nbits2, nsym, nsam); 36 37 states.fc = states.Fs/4; 38 if states.M == 2 39 states.f(1) = states.fc - Rs/2; 40 states.f(2) = states.fc + Rs/2; 41 else 42 states.f(1) = states.fc - 3*Rs/2; 43 states.f(2) = states.fc - Rs/2; 44 states.f(3) = states.fc + Rs/2; 45 states.f(4) = states.fc + 3*Rs/2; 46 end 47endfunction 48 49 50% test modulator function 51 52function tx = legacyfsk_mod(states, tx_bits) 53 Fs = states.Fs; 54 Ts = states.Ts; 55 Rs = states.Rs; 56 f = states.f; 57 M = states.M; 58 nsym = states.nsym; 59 60 tx = zeros(Ts*length(tx_bits)/log2(M),1); 61 tx_phase = 0; 62 63 step = log2(M); 64 k = 1; 65 for i=1:step:length(tx_bits) 66 if M == 2 67 tone = tx_bits(i) + 1; 68 else 69 tone = (tx_bits(i:i+1) * [2 1]') + 1; 70 end 71 tx_phase_vec = tx_phase + (1:Ts)*2*pi*f(tone)/Fs; 72 tx((k-1)*Ts+1:k*Ts) = 2.0*cos(tx_phase_vec); k++; 73 tx_phase = tx_phase_vec(Ts) - floor(tx_phase_vec(Ts)/(2*pi))*2*pi; 74 end 75 76endfunction 77 78 79function run_sim(sim_in) 80 81 frames = sim_in.frames; 82 test_frame_mode = sim_in.test_frame_mode; 83 M = sim_in.M; 84 Rs = sim_in.Rs; 85 demod = sim_in.demod; 86 EbNodB = sim_in.EbNodB; 87 timing_offset = sim_in.timing_offset; 88 89 % rx timing has been adjusted experimentally 90 91 if Rs == 4800 92 if demod == 1 93 rx_timing = 4; 94 else 95 rx_timing = 0; 96 end 97 end 98 if Rs == 2400 99 if demod == 1 100 rx_timing = 40; 101 else 102 rx_timing = 0; 103 end 104 end 105 106 % init fsk modem 107 108 more off 109 rand('state',1); 110 randn('state',1); 111 states = legacyfsk_init(M,Rs); 112 Fs = states.Fs; 113 nbits = states.nbits; 114 nbits2 = states.nbits2; 115 Ts = states.Ts; 116 nsam = states.nsam; 117 rate = states.rate; 118 119 % init analog FM modem 120 121 fm_states.Fs = Fs; 122 fm_max = fm_states.fm_max = 3E3; 123 fd = fm_states.fd = 5E3; 124 fm_states.fc = states.fc; 125 126 fm_states.pre_emp = 0; 127 fm_states.de_emp = 1; 128 fm_states.Ts = 1; 129 fm_states.output_filter = 1; 130 fm_states = analog_fm_init(fm_states); 131 [b, a] = cheby1(4, 1, 300/Fs, 'high'); % 300Hz HPF to simulate FM radios 132 133 % init sim states 134 135 rx_bits_buf = zeros(1,2*nbits2); 136 Terrs = Tbits = 0; 137 state = 0; 138 nerr_log = []; 139 140 % set up the channel noise. We have log(M)*rate payload bits/symbol 141 % we have log2(M) bits/symbol, and rate bits per payload symbol 142 % TODO: explain this better as Im confused! 143 144 EbNo = 10^(EbNodB/10); 145 EsNo = EbNo*rate*log2(M); 146 variance = states.Fs/((states.Rs)*EsNo); 147 %printf("EbNodB: %3.1f EbNo: %3.2f EsNo: %3.2f\n", EbNodB, EbNo, EsNo); 148 149 % set up the input bits 150 151 if test_frame_mode == 1 152 % test frame of bits, which we repeat for convenience when BER testing 153 test_frame = round(rand(1, nbits)); 154 tx_bits = []; 155 for i=1:frames+1 156 tx_bits = [tx_bits test_frame]; 157 end 158 end 159 if test_frame_mode == 2 160 % random bits, just to make sure sync algs work on random data 161 tx_bits = round(rand(1, nbits*(frames+1))); 162 end 163 if test_frame_mode == 3 164 % ...10101... sequence 165 tx_bits = zeros(1, nbits*(frames+1)); 166 tx_bits(1:2:length(tx_bits)) = 1; 167 end 168 169 % Manchester Encoding ----------------------------------------------------------- 170 171 % Manchester encoding, which removes DC term in baseband signal, 172 % making the waveform friendly to old-school legacy FM radios with 173 % voiceband filtering. The "code rate" is 0.5, which means we have 174 % encode one input bit into 2 output bits. The 2FSK encoder takes 175 % one input bit, the 4FSK encoder two input bits. 176 177 tx_bits_encoded = zeros(1,length(tx_bits)*2); 178 fsk2_enc = [[1 0]; [0 1]]; 179 % -1.5 1.5 1.5 -1.5 -0.5 0.5 0.5 -0.5 180 % 0 3 3 0 1 2 2 1 181 fsk4_enc = [[0 0 1 1]; [1 1 0 0]; [0 1 1 0]; [1 0 0 1]]; 182 k=1; 183 if M == 2 184 for i=1:2:length(tx_bits_encoded) 185 input_bit = tx_bits(k); k++; 186 tx_bits_encoded(i:i+1) = fsk2_enc(input_bit+1,:); 187 end 188 else 189 for i=1:4:length(tx_bits_encoded) 190 input_bits = tx_bits(k:k+1) * [2 1]'; k+=2; 191 tx_bits_encoded(i:i+3) = fsk4_enc(input_bits+1,:); 192 end 193 end 194 195 % FSK Modulator -------------------------------------------------------------- 196 197 % use ideal FSK modulator (note: need to try using analog FM modulator) 198 199 tx = legacyfsk_mod(states, tx_bits_encoded); 200 noise = sqrt(variance)*randn(length(tx),1); 201 rx = tx + noise; 202 timing_offset_samples = round(timing_offset*Ts); 203 rx = [zeros(timing_offset_samples,1); rx]; 204 205 % Demodulator ---------------------------------------------------------------------------- 206 207 if demod == 1 208 % use analog FM demodulator, aka a $40 Baofeng 209 210 [rx_out rx_bb] = analog_fm_demod(fm_states, rx'); 211 if sim_in.hpf 212 rx_out_hp = filter(b,a,rx_out); 213 else 214 rx_out_hp = rx_out; 215 end 216 rx_filt = filter(ones(1,Ts),1,rx_out_hp); 217 rx_timing_sig = rx_filt; 218 219 % TODO: for 4FSK determine amplitude/decn boundaries, choose closest to demod each symbol 220 221 end 222 223 if demod == 2 224 225 % optimal non-coherent demod at Rs 226 227 rx_timing_sig = zeros(1,length(rx)); 228 for m=1:M 229 phi_vec = (1:length(rx))*2*pi*states.f(m)/Fs; 230 dc = rx' .* exp(-j*phi_vec); 231 rx_filt(m,:) = abs(filter(ones(1,Ts),1,dc)); 232 rx_timing_sig = rx_timing_sig + rx_filt(m,1:length(rx)); 233 end 234 end 235 236 % Fine timing estimation ------------------------------------------------------ 237 238 % Estimate fine timing using line at Rs/2 that Manchester encoding provides 239 % We need this to sync up to Manchester codewords. TODO plot signal and 240 % timing "line" we extract 241 242 Np = length(rx_timing_sig); 243 w = 2*pi*(Rs)/Fs; 244 x = (rx_timing_sig .^ 2) * exp(-j*w*(0:Np-1))'; 245 norm_rx_timing = angle(x)/(2*pi) - 0.42; 246 %rx_timing = round(norm_rx_timing*Ts); 247 %printf("norm_rx_timing: %4.4f rx_timing: %d\n", norm_rx_timing, rx_timing); 248 249 % Max likelihood decoding of Manchester encoded symbols. Search 250 % through all ML possibilities to extract bits. Use energy (filter 251 % output sq) 252 253 % Manchester Decoding -------------------------------------------------------- 254 255 if M == 2 256 if demod == 1 257 258 % sample at optimum instant 259 260 [tmp l] = size(rx_filt); 261 rx_filt_dec = rx_filt(:, Ts+rx_timing:Ts:l); 262 263 [tmp l] = size(rx_filt_dec); 264 rx_bits = zeros(1,l); 265 k = 1; 266 for i=1:2:l-1 267 ml = [rx_filt_dec(i)-rx_filt_dec(i+1) -rx_filt_dec(i)+rx_filt_dec(i+1)]; 268 [mx mx_ind] = max(ml); 269 rx_bits(k) = mx_ind-1; k++; 270 end 271 end 272 273 if demod == 2 274 275 % sample at optimum instant 276 277 [tmp l] = size(rx_filt); 278 rx_filt_dec = rx_filt(:, Ts+rx_timing:Ts:l); 279 280 [tmp l] = size(rx_filt_dec); 281 rx_bits = zeros(1,l); 282 k = 1; 283 for i=1:2:l-1 284 %ml = [rx_filt_dec(2,i)*rx_filt_dec(1,i+1) rx_filt_dec(1,i)*rx_filt_dec(2,i+1)]; 285 ml = [rx_filt_dec(2,i)+rx_filt_dec(1,i+1) rx_filt_dec(1,i)+rx_filt_dec(2,i+1)]; 286 [mx mx_ind] = max(ml); 287 rx_bits(k) = mx_ind-1; k++; 288 end 289 end 290 else % M == 4 291 if demod == 1 292 % TODO: 4FSK version of demod 293 rx_bits=tx_bits; 294 end 295 if demod == 2 296 % sample at optimal instant 297 298 [tmp l] = size(rx_filt); 299 rx_filt_dec = rx_filt(:, Ts+rx_timing:Ts:l); 300 [tmp l] = size(rx_filt_dec); 301 rx_bits = zeros(1,l); 302 303 k = 1; 304 fsk4_dec = [[0 0]; [0 1]; [1 0]; [1 1]]; 305 for i=1:2:l-1 306 %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)]; 307 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))]; 308 [mx mx_ind] = max(ml); 309 rx_bits(k:k+1) = fsk4_dec(mx_ind,:); k+=2; 310 end 311 end 312 end 313 314 % useful for getting decoding right 315 %tx_bits(1:20) 316 %rx_bits(1:20) 317 318 % Frame sync and BER logic ------------------------------------------------------------- 319 320 st = 1; 321 for f=1:frames 322 323 % extract nin bits 324 325 nin = nbits; 326 en = st + nin - 1; 327 328 rx_bits_buf(1:nbits) = rx_bits_buf(nbits+1:2*nbits); 329 rx_bits_buf(nbits+1:2*nbits) = rx_bits(st:en); 330 331 st += nin; 332 333 % frame sync based on min BER 334 335 if test_frame_mode == 1 336 nerrs_min = nbits; 337 next_state = state; 338 if state == 0 339 for i=1:nbits 340 error_positions = xor(rx_bits_buf(i:nbits+i-1), test_frame); 341 nerrs = sum(error_positions); 342 %printf("i: %d nerrs: %d nerrs_min: %d \n", i, nerrs, nerrs_min); 343 if nerrs < nerrs_min 344 nerrs_min = nerrs; 345 coarse_offset = i; 346 end 347 end 348 if nerrs_min < 3 349 next_state = 1; 350 %printf("%d %d\n", coarse_offset, nerrs_min); 351 end 352 end 353 354 if state == 1 355 error_positions = xor(rx_bits_buf(coarse_offset:coarse_offset+nbits-1), test_frame); 356 nerrs = sum(error_positions); 357 Terrs += nerrs; 358 Tbits += nbits; 359 nerr_log = [nerr_log nerrs]; 360 end 361 362 state = next_state; 363 364 end 365 end 366 367 if test_frame_mode == 1 368 if sim_in.verbose 369 printf(" demod: %d frames: %d EbNodB: %3.1f Tbits: %d Terrs: %d BER %4.3f\n", demod, frames, EbNodB, Tbits, Terrs, Terrs/Tbits); 370 else 371 printf(" EbNodB: %3.1f BER %4.3f\n", EbNodB, Terrs/Tbits); 372 end 373 end 374 375 % Bunch O'plots -------------------------------------------------------------- 376 377 close all; 378 379 st = 1; en=20; 380 381 Tx=fft(tx, Fs); 382 TxdB = 20*log10(abs(Tx(1:Fs/2))); 383 figure(1) 384 clf; 385 plot(TxdB) 386 axis([1 Fs/2 (max(TxdB)-100) max(TxdB)]) 387 title('Tx Spectrum'); 388 389 figure(2) 390 clf 391 if demod == 1 392 subplot(211) 393 plot(rx_filt(st*Ts:en*Ts)); 394 title('After integrator'); 395 subplot(212) 396 plot(rx_filt_dec(st:en),'+'); 397 title('Decimated output'); 398 end 399 if demod == 2 400 subplot(211); 401 plot(rx_filt(1,st*Ts:en*Ts)); 402 hold on; 403 plot(rx_filt(2,st*Ts:en*Ts),'g'); 404 if M == 4 405 plot(rx_filt(3,st*Ts:en*Ts),'c'); 406 plot(rx_filt(4,st*Ts:en*Ts),'r'); 407 end 408 hold off; 409 title('Output of each filter'); 410 subplot(212); 411 plot(rx_filt_dec(1,st:en),'+'); 412 hold on; 413 plot(rx_filt_dec(2,st:en),'g+'); 414 if M == 4 415 plot(rx_filt_dec(3,st:en),'c+'); 416 plot(rx_filt_dec(4,st:en),'r+'); 417 end 418 hold off; 419 title('Decimated output of each filter'); 420 end 421 422 figure(3) 423 clf; 424 subplot(211) 425 plot(rx_timing_sig(st*Ts:en*Ts).^2) 426 title('rx-timing-sig') 427 subplot(212) 428 F = abs(fft(rx_timing_sig(1:Fs))); 429 plot(F(100:8000)) 430 title('FFT of rx-timing-sig') 431 432 if demod == 1 433 figure(4); 434 clf; 435 h = fft(rx_out, Fs); 436 hdB = 20*log10(abs(h)); 437 plot(hdB(1:4000)) 438 title('Spectrum of baseband modem signal after analog FM demod'); 439 axis([1 4000 (max(hdB)-40) max(hdB)]) 440 end 441 442 if demod == 1 443 figure(5) 444 clf; 445 subplot(211) 446 plot(rx_out(st*Ts:en*Ts)); 447 title('baseband modem signal after analog FM demod'); 448 subplot(212) 449 plot(rx_out_hp(st*Ts:en*Ts)); 450 title('baseband modem signal after 300Hz filter'); 451 end 452end 453 454 455% Run various permutations of simulation here --------------------------------------- 456 457function run_single 458 459 sim_in.frames = 100; 460 sim_in.test_frame_mode = 1; 461 sim_in.M = 2; 462 sim_in.Rs = 2400; 463 sim_in.demod = 1; 464 sim_in.EbNodB = 15; 465 sim_in.timing_offset = 0.0; 466 sim_in.hpf = 1; 467 sim_in.verbose = 1; 468 469 run_sim(sim_in); 470endfunction 471 472 473function run_lots 474 475 % adjusted a few scenarios for about 2% BER so we can compare 476 477 sim_in.frames = 100; 478 sim_in.test_frame_mode = 1; 479 sim_in.M = 2; 480 sim_in.Rs = 4800; 481 sim_in.demod = 1; 482 sim_in.EbNodB = 12; 483 sim_in.timing_offset = 0.0; 484 sim_in.hpf = 1; 485 sim_in.verbose = 0; 486 487 printf("Rs=4800 2FSK ideal demod\n"); 488 sim_in.EbNodB = 8.5; sim_in.demod = 2; run_sim(sim_in); 489 printf("Rs=4800 2FSK analog FM demod, not too shabby and pushes 2400bit/s thru a $40 HT!\n"); 490 sim_in.EbNodB = 12; sim_in.demod = 1; run_sim(sim_in); 491 printf("Rs=2400 2FSK analog FM demod, needs more power for same BER! Che?\n"); 492 sim_in.Rs = 2400; sim_in.EbNodB = 15; run_sim(sim_in); 493 printf("Hmm, doesnt improve with no 300Hz HPF, maybe due to less deviation?\n"); 494 sim_in.hpf = 0; run_sim(sim_in); 495 printf("Rs=2400 4FSK ideal demod, nice low Eb/No!\n"); 496 sim_in.demod = 2; sim_in.M = 4; sim_in.Rs = 2400; sim_in.EbNodB = 6; run_sim(sim_in); 497endfunction 498 499%run_single; 500run_lots; 501