From: drowe67 Date: Thu, 19 Jan 2017 08:37:35 +0000 (+0000) Subject: working for a variety of channel impairments, can also BER draw curves X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=5afcfb83853382643612ba97b72a474d1dfd740b;p=freetel-svn-tracking.git working for a variety of channel impairments, can also BER draw curves git-svn-id: https://svn.code.sf.net/p/freetel/code@2991 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/oqpsk.m b/codec2-dev/octave/oqpsk.m index 897bd128..8ae571b5 100644 --- a/codec2-dev/octave/oqpsk.m +++ b/codec2-dev/octave/oqpsk.m @@ -8,8 +8,9 @@ rand('state',1); randn('state',1); graphics_toolkit ("gnuplot"); format +more off; -% init nodem states +% init modem states function oqpsk_states = oqpsk_init(oqpsk_states, Rs) @@ -146,7 +147,7 @@ function [rx_symb rx_int filt_log dco_log timing_adj Toff] = oqpsk_demod(oqpsk_s xr_log = []; xi_log = []; w_log = []; timing_clock_phase = 0; - timing_angle = pi; % XXX + timing_angle = 0; % XXX timing_angle_log = zeros(1,nsam); end @@ -167,8 +168,8 @@ function [rx_symb rx_int filt_log dco_log timing_adj Toff] = oqpsk_demod(oqpsk_s xi = abs(imag(rx_int(l:l+tw-1))); w = exp(j*(l:l+tw-1)*2*pi*Rs/Fs); X = xr * w'; - %timing_clock_phase = timing_angle = angle(X); - timing_clock_phase = timing_angle = pi; % XXX + timing_clock_phase = timing_angle = angle(X); + %timing_clock_phase = timing_angle = 0; % XXX k++; xr_log = [xr_log xr]; xi_log = [xi_log xi]; @@ -187,7 +188,7 @@ function [rx_symb rx_int filt_log dco_log timing_adj Toff] = oqpsk_demod(oqpsk_s ph_err = sign(real(rx_int(i))*imag(rx_int(i)))*cos(timing_clock_phase); lower = ph_err*G2 + lower; filt = ph_err*G1 + lower; - dco_log(i) = dco = pi; % XXX + dco_log(i) = dco; % XXX dco = dco + filt; filt_log(i) = filt; @@ -244,7 +245,8 @@ function sim_out = oqpsk_test(sim_in) M = oqpsk_states.M; Fs = oqpsk_states.Fs; Rs = oqpsk_states.Rs; - + sample_clock_offset_ppm = sim_in.sample_clock_offset_ppm; + tx_testframe = round(rand(1, bitspertestframe)); ntestframes = floor(nbits/bitspertestframe); tx_bits = []; @@ -258,25 +260,47 @@ function sim_out = oqpsk_test(sim_in) variance = Fs/(Rs*EbNo*oqpsk_states.bps); [tx tx_symb] = oqpsk_mod(oqpsk_states, tx_bits); + if sample_clock_offset_ppm + tx = resample(tx, 1E6, 1E6-sample_clock_offset_ppm); + end nsam = length(tx); + phi = sim_in.phase_offset + 2*pi*sim_in.freq_offset*(1:nsam)/M; + noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); st = 1+sim_in.timing_offset; en = length(tx); - rx = tx(st:en)*exp(j*sim_in.phase_offset) + noise(st:en); + rx = tx(st:en).*exp(j*phi(st:en)) + noise(st:en); [rx_symb rx_int filt_log dco_log timing_adj Toff] = oqpsk_demod(oqpsk_states, rx); - % Determine ambiguities: - % a) could be m*pi/2 rotations of phase by phase est + % OK so the phase and timing estimators get us close (e.g. a good + % scatter diagram), but no banana just yet. One problem is the + % PLL can lock up on mulitples of pi/2. Combinations of phase + % offsets can confuse the timing estimator. On tricky example is a + % phase offset of pi/2 which swaps I & Q, and with OQPSK (unlike + % MSK and friends) we can't easily tell which is I and which is Q + % after a phase rotation, e.g. could be IQIQIQI or QIQIQIQ + + % So we need to determine the ambiguities: + % a) could be m*pi/2 rotations of phase % b) could be I and Q swapped by timing est % c) time alignment of test frame - + nsymb = bitspertestframe/oqpsk_states.bps; nrx_symb = length(rx_symb); rx_bits = zeros(1, bitspertestframe); atx_symb = tx_symb(1:nsymb); - % correlate with I and Q tx sequences at various offsets + % Treat I and Q as separate sequences, each with their own unique + % word. In our case the UW is the whole test frame. Correlate rx + % sequence with tx sequence at each possible offset through the + % received symbols to find the test frames. Note we also + % correlate I of tx with Q of rx to trap any IQ swaps. + + % The sign of the I and Q correlation lets us sort out the pi/2 + % phase rotation issue. + + nerrs_tot = 0; nbits_tot = 0; max_corr = real(atx_symb) * real(atx_symb)'; for offset=2:nrx_symb-nsymb+1 @@ -288,121 +312,142 @@ function sim_out = oqpsk_test(sim_in) %offset, corr_ii(offset), corr_qq(offset), corr_iq(offset), corr_qi(offset)); if abs(corr_ii(offset)) > 0.8 + % no IQ swap, or time offset - arx_symb = real(rx_symb(offset:offset+nsymb-1)) + j*imag(rx_symb(offset:offset+nsymb-1)); + + i_sign = sign(corr_ii(offset)); + q_sign = sign(corr_qq(offset)); + arx_symb = i_sign*real(rx_symb(offset:offset+nsymb-1)) + j*q_sign*imag(rx_symb(offset:offset+nsymb-1)); + + for i=1:nsymb + rx_bits(2*i-1:2*i) = qpsk_demod(arx_symb(i)*exp(-j*pi/4)); + end + nerrs = sum(xor(tx_testframe, rx_bits)); + if verbose > 2 + printf("offset: %5d swap: %d i_sign: % 2.1f q_sign: % 2.1f nerr: %d\n", + offset, 0, i_sign, q_sign, nerrs); + end + nerrs_tot += nerrs; + nbits_tot += bitspertestframe; end + if abs(corr_qi(offset)) > 0.8 + % IQ swap, I part in Q part of symbol before + i_sign = sign(corr_iq(offset-1)); q_sign = sign(corr_qi(offset)); arx_symb = i_sign*imag(rx_symb(offset-1:offset+nsymb-2)) + j*q_sign*real(rx_symb(offset:offset+nsymb-1)); + for i=1:nsymb rx_bits(2*i-1:2*i) = qpsk_demod(arx_symb(i)*exp(-j*pi/4)); end - nerr = sum(xor(tx_testframe, rx_bits)); - printf("offset: %5d swap: %d i_sign: % 2.1f q_sign: % 2.1f nerr: %d\n", - offset, 1, i_sign, q_sign, nerr); - end - end - - -#{ - for i=1:nrx_symb-nsymb+1 - for k=0:3 - phase_amb = exp(j*k*pi/2); - arx_symb = rx_symb(i:i+nsymb-1) .* phase_amb; - for l=1:nsymb - rx_bits(2*l-1:2*l) = qpsk_demod(arx_symb(l).*exp(-j*pi/4)); - end - nerr = sum(xor(tx_testframe, rx_bits)); - if nerr == 0 - printf("i: %d k: %d nerr: %d\n", i, k, nerr); + nerrs = sum(xor(tx_testframe, rx_bits)); + if verbose > 1 + printf("offset: %5d swap: %d i_sign: % 2.1f q_sign: % 2.1f nerr: %d\n", + offset, 1, i_sign, q_sign, nerrs); end + nerrs_tot += nerrs; + nbits_tot += bitspertestframe; end end -#} - TERvec(ne) = 0; - BERvec(ne) = 0; -#{ - Nerrs_min = nbits; Nbits_min = nbits; l = length(rx_bits); - for i=1:1 - Nerrs = sum(xor(rx_bits(i:l), tx_bits(1:l-i+1))); - if Nerrs < Nerrs_min - Nerrs_min = Nerrs; - Nbits_min = l; - end - end - - TERvec(ne) = Nerrs_min; - BERvec(ne) = Nerrs_min/nbits; + TERvec(ne) = nerrs_tot; + BERvec(ne) = nerrs_tot/nbits_tot; if verbose > 0 printf("EbNo dB: %3.1f Nbits: %d Nerrs: %d BER: %4.3f BER Theory: %4.3f\n", - aEbNodB, nbits, Nerrs_min, BERvec(ne), 0.5*erfc(sqrt(EbNo))); + aEbNodB, nbits_tot, nerrs_tot, BERvec(ne), 0.5*erfc(sqrt(EbNo))); end -#} - figure(1); clf; - subplot(211) - stem(real(tx_symb)) - title('Tx symbols'); - subplot(212) - stem(imag(tx_symb)) - - figure(2); clf; - f = fftshift(fft(rx)); - Tx = 20*log10(abs(f)); - plot((1:length(f))*Fs/length(f) - Fs/2, Tx) - grid; - title('OQPSK Demodulator Input Spectrum'); - - figure(3); clf; - nplot = min(16, nbits/oqpsk_states.bps); - title('Rx Integrator'); - subplot(211) - stem(real(rx_int(1:nplot*M))) - axis([1 nplot*M -1 1]) - subplot(212) - stem(imag(rx_int(1:nplot*M))) - axis([1 nplot*M -1 1]) - - figure(4); clf; - subplot(211); - plot(filt_log); - title('PLL filter') - subplot(212); - plot(dco_log); - title('PLL DCO phase'); - - figure(5); clf; - subplot(211) - plot(timing_adj); - title('Timing est'); - subplot(212) - plot(Toff); - title('Timing est unwrap'); - - figure(6); clf; - st = floor(0.5*nrx_symb); - plot(rx_symb(st:nrx_symb), '+'); - title('Scatter Diagram'); - axis([-1.5 1.5 -1.5 1.5]) - - figure(7); clf; - subplot(211); - stem(real(arx_symb)); - title('Rx symbols') - subplot(212); - stem(imag(arx_symb)); - - figure(8); clf; - subplot(211) - stem(tx_testframe(1:min(20,length(rx_bits)))) - title('Tx Bits') - subplot(212) - stem(rx_bits(1:min(20,length(rx_bits)))) - title('Rx Bits') + if find(sim_in.plots == 1) + figure(1); clf; + subplot(211) + stem(real(tx_symb)) + title('Tx symbols'); + subplot(212) + stem(imag(tx_symb)) + end + + if find(sim_in.plots == 2) + figure(2); clf; + f = fftshift(fft(rx)); + Tx = 20*log10(abs(f)); + plot((1:length(f))*Fs/length(f) - Fs/2, Tx) + grid; + title('OQPSK Demodulator Input Spectrum'); + end + + if find(sim_in.plots == 3) + figure(3); clf; + nplot = min(16, nbits/oqpsk_states.bps); + title('Rx Integrator'); + subplot(211) + stem(real(rx_int(1:nplot*M))) + axis([1 nplot*M -1 1]) + subplot(212) + stem(imag(rx_int(1:nplot*M))) + axis([1 nplot*M -1 1]) + end + + if find(sim_in.plots == 4) + figure(4); clf; + subplot(211); + plot(filt_log); + title('PLL filter') + subplot(212); + plot(dco_log); + title('PLL DCO phase'); + end + + if find(sim_in.plots == 5) + figure(5); clf; + subplot(211) + plot(timing_adj); + title('Timing est'); + subplot(212) + plot(Toff); + title('Timing est unwrap'); + end + + if find(sim_in.plots == 6) + figure(6); clf; + st = floor(0.5*nrx_symb); + plot(rx_symb(st:nrx_symb), '+'); + title('Scatter Diagram'); + axis([-1.5 1.5 -1.5 1.5]) + end + + if find(sim_in.plots == 7) + figure(7); clf; + subplot(211) + plot(corr_ii); + axis([1 length(corr_ii) -1.2 1.2]); + title('corr ii'); + subplot(212) + plot(corr_qi); + axis([1 length(corr_ii) -1.2 1.2]); + title('corr qi'); + end + + if find(sim_in.plots == 8) + figure(8); clf; + subplot(211); + stem(real(arx_symb)); + title('Rx symbols') + subplot(212); + stem(imag(arx_symb)); + end + + if find(sim_in.plots == 9) + figure(9); clf; + subplot(211) + stem(tx_testframe(1:min(20,length(rx_bits)))) + title('Tx Bits') + subplot(212) + stem(rx_bits(1:min(20,length(rx_bits)))) + title('Rx Bits') + end end sim_out.TERvec = TERvec; @@ -416,11 +461,14 @@ function run_oqpsk_single sim_in.phase_est = 1; sim_in.timing_est = 1; sim_in.bitspertestframe = 100; - sim_in.nbits = 1000; - sim_in.EbNodB = 100; - sim_in.verbose = 2; - sim_in.phase_offset = pi/2; - sim_in.timing_offset = 0; + sim_in.nbits = 10000; + sim_in.EbNodB = 4; + sim_in.verbose = 1; + sim_in.phase_offset = 3*pi/4; % in radians + sim_in.timing_offset = 4; % in samples 0..M-1 + sim_in.freq_offset = 0.001; % fraction of Symbol Rate + sim_in.plots = [2 4 5 6 7]; + sim_in.sample_clock_offset_ppm = 100; sim_out = oqpsk_test(sim_in); endfunction @@ -430,211 +478,42 @@ endfunction function run_oqpsk_curves sim_in.coherent_demod = 1; - sim_in.nsym = 48000; - sim_in.EbNodB = 2:10; + sim_in.EbNodB = 2:8; sim_in.verbose = 1; + sim_in.phase_est = 1; + sim_in.timing_est = 1; + sim_in.bitspertestframe = 100; + sim_in.nbits = 50000; + sim_in.phase_offset = 3*pi/4; % in radians + sim_in.timing_offset = 4; % in samples 0..M-1 + sim_in.freq_offset = 0.001; % fraction of Symbol Rate + sim_in.plots = []; + sim_in.sample_clock_offset_ppm = 0; oqpsk_coh = oqpsk_test(sim_in); - sim_in.coherent_demod = 0; - oqpsk_noncoh = oqpsk_test(sim_in); - Rs = oqpsk_coh.Rs; EbNo = 10 .^ (sim_in.EbNodB/10); - alpha = 0.75; % guess for BT=0.5 OQPSK - oqpsk_theory.BERvec = 0.5*erfc(sqrt(alpha*EbNo)); + oqpsk_theory.BERvec = 0.5*erfc(sqrt(EbNo)); % BER v Eb/No curves figure; clf; - semilogy(sim_in.EbNodB, oqpsk_theory.BERvec,'r;OQPSK theory;') + semilogy(sim_in.EbNodB, oqpsk_theory.BERvec,'r+-;OQPSK theory;') hold on; - semilogy(sim_in.EbNodB, oqpsk_coh.BERvec,'g;OQPSK sim coherent;') - semilogy(sim_in.EbNodB, oqpsk_noncoh.BERvec,'b;OQPSK sim non-coherent;') + semilogy(sim_in.EbNodB, oqpsk_coh.BERvec,'g+-;OQPSK sim;') hold off; grid("minor"); axis([min(sim_in.EbNodB) max(sim_in.EbNodB) 1E-4 1]) legend("boxoff"); xlabel("Eb/No (dB)"); ylabel("Bit Error Rate (BER)") - - % BER v C/No (1 Hz noise BW and Eb=C/Rs=1/Rs) - % Eb/No = (C/Rs)/(1/(N/B)) - % C/N = (Eb/No)*(Rs/B) - - RsOnB_dB = 10*log10(Rs/1); - figure; - clf; - semilogy(sim_in.EbNodB+RsOnB_dB, oqpsk_theory.BERvec,'r;OQPSK theory;') - hold on; - semilogy(sim_in.EbNodB+RsOnB_dB, oqpsk_coh.BERvec,'g;OQPSK sim coherent;') - semilogy(sim_in.EbNodB+RsOnB_dB, oqpsk_noncoh.BERvec,'b;OQPSK sim non-coherent;') - hold off; - grid("minor"); - axis([min(sim_in.EbNodB+RsOnB_dB) max(sim_in.EbNodB+RsOnB_dB) 1E-4 1]) - legend("boxoff"); - xlabel("C/No for Rs=4800 bit/s and 1 Hz noise bandwidth (dB)"); - ylabel("Bit Error Rate (BER)") - endfunction - -% attempt to perform "coarse sync" sync with the received frames, we -% check each frame for the best coarse sync position. Brute force -% approach, that would be changed for a real demod which has some -% sort of unique word. Start looking for valid frames 1 frame -% after start of pre-amble to give PLL time to lock - -function [total_errors total_bits Nerrs_log Nerrs_all_log errors_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits) - - Nerrs_log = zeros(1, nframes_rx); - Nerrs_all_log = zeros(1, nframes_rx); - total_errors = 0; - total_bits = 0; - framesize = length(tx_frame); - errors_log = []; - - for f=2:nframes_rx-1 - Nerrs_min = framesize; - for i=1:framesize; - st = (f-1)*framesize+i; en = st+framesize-1; - errors = xor(rx_bits(st:en), tx_frame); - Nerrs = sum(errors); - if Nerrs < Nerrs_min - Nerrs_min = Nerrs; - errors_min = errors; - end - end - Nerrs_all_log(f) = Nerrs_min; - if Nerrs_min/framesize < 0.1 - errors_log = [errors_log errors_min]; - Nerrs_log(f) = Nerrs_min; - total_errors += Nerrs_min; - total_bits += framesize; - end - end -endfunction - - -function plot_spectrum(oqpsk_states, rx, preamble_location, title_str) - Fs = oqpsk_states.Fs; - st = preamble_location + oqpsk_states.npreamble*oqpsk_states.M; - sig = rx(st:st+Fs*0.5); - h = hanning(length(sig))'; - Rx=20*log10(abs(fftshift(fft(sig .* h, Fs)))); - figure; - plot(-Fs/2:Fs/2-1,Rx); - grid("minor"); - xlabel('Hz'); - ylabel('dB'); - topy = ceil(max(Rx)/10)*10; - axis([-4000 4000 topy-50 topy+10]) - title(title_str); -endfunction - - -% Give the demod a hard time: frequency, phase, time offsets, sample clock difference - -function run_test_channel_impairments - Rs = 4800; - verbose = 1; - aEbNodB = 10; - phase_offset = 0; - freq_offset = 0; - timing_offset = 1; - sample_clock_offset_ppm = 0; - nsym = Rs*10; - - oqpsk_states.verbose = verbose; - oqpsk_states.coherent_demod = 1; - oqpsk_states.phase_track = 1; - oqpsk_states = oqpsk_init(oqpsk_states, Rs); - Fs = oqpsk_states.Fs; - Rs = oqpsk_states.Rs; - M = oqpsk_states.M; - - % A frame consists of nsym random data bits. - - framesize = 480; - nframes = floor(nsym/framesize); - tx_frame = round(rand(1, framesize)); - tx_bits = []; - for i=1:nframes - tx_bits = [tx_bits tx_frame]; - end - - tx = oqpsk_mod(oqpsk_states, tx_bits); - - tx = resample(tx, 1E6, 1E6-sample_clock_offset_ppm); - tx = [zeros(1,timing_offset) tx]; - nsam = length(tx); - - EbNo = 10^(aEbNodB/10); - variance = Fs/(Rs*EbNo*oqpsk_states.bps); - noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam)); - w = (0:nsam-1)*2*pi*freq_offset/Fs + phase_offset; - - rx = sqrt(2)*tx.*exp(j*w) + noise; - - % optional dump to file - - if 0 - fc = 1500; gain = 10000; - wc = 2*pi*fc/Fs; - w1 = exp(j*wc*(1:nsam)); - rx1 = gain*real(rx .* w1); - fout = fopen("rx_6dB.raw","wb"); - fwrite(fout, rx1, "short"); - fclose(fout); - end - - % printf("ntx: %d nrx: %d ntx_bits: %d\n", length(tx), length(rx), length(tx_bits)); - - [rx_bits rx_out rx_symb] = oqpsk_demod(oqpsk_states, rx); - nframes_rx = length(rx_bits)/framesize; - - % printf("ntx: %d nrx: %d ntx_bits: %d nrx_bits: %d\n", length(tx), length(rx), length(tx_bits), length(rx_bits)); - - % try different possible phase ambiguities, a UW would sort this out in the real world - - for i=1:4 - phase_amb = i*pi/2; - rx_bits = []; - for k=1:length(rx_symb) - rx_bits = [rx_bits qpsk_demod(rx_symb(k)*exp(j*(-pi/4+phase_amb)))]; - end - - [total_errors total_bits Nerrs_log Nerrs_all_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits); - - if total_bits - ber = total_errors/total_bits; - printf("Eb/No: %3.1f f_off: %4.1f ph_off: %4.3f ph_amb: %4.3f Nframes: %d Nbits: %d Nerrs: %d BER: %f\n", - aEbNodB, freq_offset, phase_amb, phase_offset, nframes_rx, total_bits, total_errors, ber); - - figure(1); clf; - subplot(211) - plot(Nerrs_log,'r;errors/frame counted for BER;'); - hold on; - plot(Nerrs_all_log,'g;all errors/frame;'); - hold off; - legend("boxoff"); - title('Bit Errors') - subplot(212) - stem(real(cumsum(Nerrs_log))) - title('Cumulative Bit Errors') - - end - end - -endfunction - - % Choose one of these to run ------------------------------------------ -run_oqpsk_single -%run_test_channel_impairments -%run_oqpsk_curves -%run_oqpsk_init +%run_oqpsk_single +run_oqpsk_curves