From d6445a80746a21ecc9382281263916c09c7f91ce Mon Sep 17 00:00:00 2001 From: baobrien Date: Thu, 17 Sep 2015 01:28:09 +0000 Subject: [PATCH] 4FSK demod now getting close to non-coherent theory. git-svn-id: https://svn.code.sf.net/p/freetel/code@2339 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fsk4.m | 135 +++++++++++++++++++++++++++++++-------- 1 file changed, 107 insertions(+), 28 deletions(-) diff --git a/codec2-dev/octave/fsk4.m b/codec2-dev/octave/fsk4.m index ef5d740d..a21ffd00 100644 --- a/codec2-dev/octave/fsk4.m +++ b/codec2-dev/octave/fsk4.m @@ -49,7 +49,7 @@ function d = idmp(data, M, offset) end endfunction -function [tx, tx_filt] = fsk4_mod(fsk4_states, tx_bits) +function [tx, tx_filt, tx_stream] = fsk4_mod(fsk4_states, tx_bits) hbits = tx_bits(1:2:length(tx_bits)); lbits = tx_bits(2:2:length(tx_bits)); %Pad odd bit lengths @@ -66,32 +66,104 @@ function [tx, tx_filt] = fsk4_mod(fsk4_states, tx_bits) tx_stream(1+(i-1)*M:i*M) = fsk4_states.symmap(tx_symbols(i)); end tx_filt = filter(fsk4_states.tx_filter, 1, tx_stream); - tx_filt = tx_filt / max(tx_filt); + %tx_filt = tx_filt / max(tx_filt); tx = analog_fm_mod(fsk4_states.fm_states, tx_filt); endfunction %non-coherent demod based on a paper I found on IEEE xplore. Paper claims it is ~1db under coherent. % I don't think it works % Paper is titled "ALL- DIGITAL PSEUDO- COHERENT (PC) FSK MODEMS" -function sym = fsk4_demod_thing(fsk4_states, rx) +function bits = fsk4_demod_thing(fsk4_states, rx) global fsk4_symbols; Fs = fsk4_states.Fs; + t = (1:length(rx)); - shiftup = exp(j*2*pi*(1/4)*t); - - rx_up = real(rx.*shiftup); - symup = fsk4_symbols + Fs/4; - sym1m = exp(j*2*pi*(symup(1)/Fs)*t).*rx_up; - sym2m = exp(j*2*pi*(symup(2)/Fs)*t).*rx_up; - sym3m = exp(j*2*pi*(symup(3)/Fs)*t).*rx_up; - sym4m = exp(j*2*pi*(symup(4)/Fs)*t).*rx_up; - sym1m = idmp(sym1m,10); sym1m = (real(sym1m).^2+imag(sym1m).^2); - sym2m = idmp(sym2m,10); sym2m = (real(sym2m).^2+imag(sym2m).^2); - sym3m = idmp(sym3m,10); sym3m = (real(sym3m).^2+imag(sym3m).^2); - sym4m = idmp(sym4m,10); sym4m = (real(sym4m).^2+imag(sym4m).^2); + %shiftup = exp(j*2*pi*(1/4)*t); + + ffilt = 5000 + rx_filter = fir1(300,ffilt/fsk4_states.Fs); + %rx_filter = [1]; + tx_filter = fsk4_states.tx_filter; + %rx_up = filter(rx_filter,1,rx); + + rx_up = rx;%real(rx.*shiftup); + symup = fsk4_symbols;% + Fs/4; + sym1m = filter(rx_filter,1,exp(j*2*pi*(symup(1)/Fs)*t).*rx_up); + sym2m = filter(rx_filter,1,exp(j*2*pi*(symup(2)/Fs)*t).*rx_up); + sym3m = filter(rx_filter,1,exp(j*2*pi*(symup(3)/Fs)*t).*rx_up); + sym4m = filter(rx_filter,1,exp(j*2*pi*(symup(4)/Fs)*t).*rx_up); + sym1m = idmp(sym1m,20); sym1m = abs(sym1m);%(real(sym1m).^2+imag(sym1m).^2); + sym2m = idmp(sym2m,20); sym2m = abs(sym2m);%(real(sym2m).^2+imag(sym2m).^2); + sym3m = idmp(sym3m,20); sym3m = abs(sym3m);%(real(sym3m).^2+imag(sym3m).^2); + sym4m = idmp(sym4m,20); sym4m = abs(sym4m);%(real(sym4m).^2+imag(sym4m).^2); sym = sym1m*-3 + sym2m*-1 + sym3m*1 + sym4m*3; - plot(abs(sym)(1:2000)); - hist(sym(1:2:length(sym)),30); + %figure(1); + %plot((1:2000),abs(sym1m)(1:2:4000),(1:2000),abs(sym2m)(1:2:4000),(1:2000),abs(sym3m)(1:2:4000),(1:2000),abs(sym4m)(1:2:4000)); + figure(2); + plot((1:2000),sym1m(1:2000),(1:2000),sym2m(1:2000),(1:2000),sym3m(1:2000),(1:2000),sym4m(1:2000)); + + [x iv] = max([sym1m; sym2m; sym3m; sym4m;]); + bits = zeros(1,length(iv*2)); + for i=1:length(iv) + bits(1+(i-1)*2:i*2) = [[1 1];[1 0];[0 1];[0 0]](iv(i),(1:2)); + end + figure(3); + hist(iv(1:2:length(iv)),30); +endfunction + +function bits = fsk4_demod_two(fsk4_states,rx) + global fsk4_symbols; + figure(4); + + Fs = fsk4_states.Fs; + rf = (0:(Fs/2)); + rx_filter = fir2(100 ,rf/(Fs/2),fsk4_rcf_resp(rf-1000)); + + plot(20*log10(abs(fft(rx)))); + Fs = fsk4_states.Fs; + t = (1:length(rx)); + fsk4_symbols + rx = filter(rx_filter, 1, rx); + sym1dc = exp(-j*2*pi*(fsk4_symbols(1)/Fs)*t) .* rx; + sym2dc = exp(-j*2*pi*(fsk4_symbols(2)/Fs)*t) .* rx; + sym3dc = exp(-j*2*pi*(fsk4_symbols(3)/Fs)*t) .* rx; + sym4dc = exp(-j*2*pi*(fsk4_symbols(4)/Fs)*t) .* rx; + + figure(1); + %plot(t(1:20:length(t)),abs(idmp(sym1dc,20)),t(1:20:length(t)),abs(idmp(sym2dc,20))); + + %figure(2); + %plot(t(1:20:length(t)),abs(idmp(sym3dc,20)),t(1:20:length(t)),abs(idmp(sym4dc,20))); + nsym = floor(length(rx)/fsk4_states.M) + bits = zeros(1,nsym*2); + syms = zeros(1,nsym); + + int1 = abs(idmp(sym1dc,10)); + int2 = abs(idmp(sym2dc,10)); + int3 = abs(idmp(sym3dc,10)); + int4 = abs(idmp(sym4dc,10)); + + plot((1:length(int1)),int1,(1:length(int1)),int2,(1:length(int1)),int3,(1:length(int1)),int4); + + for i=(1:nsym) + st = (i-1)*fsk4_states.M+1; + en = st+fsk4_states.M-1; + sym1i = sum(sym1dc(st:en)); + %sym1i = ; + sym2i = sum(sym2dc(st:en)); + %sym2i = real(sym2i)^2 + imag(sym2i)^2; + sym3i = sum(sym3dc(st:en)); + %sym3i = real(sym3i)^2 + imag(sym3i)^2; + sym4i = sum(sym4dc(st:en)); + %sym4i = real(sym4i)^2 + imag(sym4i)^2; + %[v iv] = max(abs([sym4i sym3i sym2i sym1i])); + [v iv] = max([int4(i*2) int3(i*2) int2(i*2) int1(i*2)]); + syms(i) = iv; + bits(1+(i-1)*2:i*2) = [[1 1];[1 0];[0 1];[0 0]](iv,(1:2)); + end + figure(3); + hist(syms); + endfunction %incoherent demod loosly based on another paper. Works, more or less. @@ -99,19 +171,20 @@ endfunction function [bits err] = fsk4_demod_fmrid(fsk4_states, rx) rxd = analog_fm_demod(fsk4_states.fm_states,rx); - %rx_filt = filter(fsk4_states.tx_filter, 1, rxd); + % rx_filt = filter(fsk4_states.tx_filter, 1, rxd); rx_filt=rxd; sym = afsym = idmp(rx_filt,fsk4_states.M/2); % Demod symbol map. I should probably figure a better way to do this. % After integrating, the high symbol tends to be about 7.5 - dmsyms = rot90(fsk4_states.symmap * 7.5); + dmsyms = rot90(fsk4_states.symmap * 10); oddsyms = afsym(1:2:length(afsym)); evensyms = afsym(2:2:length(afsym)); + hist(evensyms); [errseven,deceven] = min(abs(evensyms - dmsyms)); [errsodd ,decodd ] = min(abs(oddsyms - dmsyms)); - + terreven = mean(errseven); terrodd = mean(errsodd ); @@ -132,35 +205,41 @@ endfunction % Bit error rate test % for a noise-free channel % now supports noisy channels -function ber = nfbert(aEbNodB) +function ber = nfbert(aEsNodB) bitcnt = 48000; offset = 29; test_bits = [zeros(1,100) rand(1,bitcnt)>.5]; %Random bits. Pad with zeros to prime the filters fsk4_states.M = 1; fsk4_states = fsk4_init(fsk4_states,2400); Fs = fsk4_states.Fs; - Rs = fsk4_states.Rs * 2; %Multiply symbol rate by 2, since we have 2 bits per symbol + Rb = fsk4_states.Rs * 2; %Multiply symbol rate by 2, since we have 2 bits per symbol tx = fsk4_mod(fsk4_states,test_bits); %add noise here %shamelessly copied from gmsk.m - - EbNo = 10^(aEbNodB/10); - variance = Fs/(Rs*EbNo); + EsNo = 10^(aEsNodB/10); + EbNo = EsNo + variance = Fs/(Rb*EbNo); nsam = length(tx); noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); rx = tx*exp(j*pi/2) + noise; - [rx_bits,rx_err] = fsk4_demod_fmrid(fsk4_states,rx); + rx_bits = fsk4_demod_thing(fsk4_states,rx); ber = 1; %thing to account for offset from input data to output data %No preamble detection yet - for offset = (25:31) + ox = 1; + for offset = (1:100) bern = sum(xor(rx_bits(offset:length(rx_bits)),test_bits(1:length(rx_bits)+1-offset)))/(bitcnt-offset); + if(bern < ber) + ox = offset; + end ber = min([ber bern]); end - plot((1:1000),rx_bits(1:1000),(1:1000),rx_err(1:1000)); + offset = ox; + %plot(xor(rx_bits(offset:length(rx_bits)),test_bits(1:length(rx_bits)+1-offset))); + %plot((1:1000),rx_bits(1:1000),(1:1000),rx_err(1:1000)); endfunction -- 2.25.1