From: drowe67 Date: Sat, 11 Apr 2015 05:10:45 +0000 (+0000) Subject: octave version of modem using lin reg for phase est, curves look OK X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=976397e37b9e41acf05a2abc7f2b1692ea47dc2e;p=freetel-svn-tracking.git octave version of modem using lin reg for phase est, curves look OK git-svn-id: https://svn.code.sf.net/p/freetel/code@2112 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index 96fbec69..0a624694 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -65,9 +65,11 @@ function sim_in = symbol_rate_init(sim_in) sim_in.Npilotsframe = Npilotsframe = 2; sim_in.Nsymbrowpilot = Nsymbrowpilot = Nsymbrow + Npilotsframe; - printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc); - printf(" There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns); - printf(" Including pilots, the frame is %d symbols long by %d carriers.\n\n", Nsymbrowpilot, Nc); + if verbose == 2 + printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc); + printf(" There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns); + printf(" Including pilots, the frame is %d symbols long by %d carriers.\n\n", Nsymbrowpilot, Nc); + end sim_in.prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc*Nchip); sim_in.prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc*Nchip); @@ -226,6 +228,7 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_ Npilotsframe = cohpsk.Npilotsframe; pilot = cohpsk.pilot; verbose = cohpsk.verbose; + coh_en = cohpsk.coh_en; % average pilots to get phase and amplitude estimates % we assume there are two samples at the start of each frame and two at the end @@ -234,12 +237,20 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_ pilot2 = [cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);]; phi_ = zeros(Nsymbrow, Nc); amp_ = zeros(Nsymbrow, Nc); - + + % corr method was used initially, but was poor at tracking fast phase + % changes. Linear regression works better. + for c=1:Nc - corr = pilot2(:,c)' * ct_symb_buf(sampling_points,c); - mag = sum(abs(ct_symb_buf(sampling_points,c))); + %corr = pilot2(:,c)' * ct_symb_buf(sampling_points,c); + %phi_(:, c) = angle(corr); + + y = ct_symb_buf(sampling_points,c) .* pilot2(:,c); + [m b] = linreg(sampling_points, y, length(sampling_points)); + yfit = m*[3 4 5 6] + b; + phi_(:, c) = angle(yfit); - phi_(:, c) = angle(corr); + mag = sum(abs(ct_symb_buf(sampling_points,c))); amp_(:, c) = mag/length(sampling_points); end @@ -251,8 +262,11 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_ for c=1:Nc for r=1:Nsymbrow i = (c-1)*Nsymbrow + r; - rx_symb(r,c) = ct_symb_buf(2+r,c)*exp(-j*phi_(r,c)); - %rx_symb(r,c) = ct_symb_buf(2+r,c); + if coh_en == 1 + rx_symb(r,c) = ct_symb_buf(2+r,c)*exp(-j*phi_(r,c)); + else + rx_symb(r,c) = ct_symb_buf(2+r,c); + end rx_symb_linear(i) = rx_symb(i); rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb(r,c)); %printf("phi_ %d %d %f %f\n", r,c,real(exp(-j*phi_(r,c))), imag(exp(-j*phi_(r,c)))); @@ -682,6 +696,7 @@ function sim_out = ber_test(sim_in) hf_sim = sim_in.hf_sim; nhfdelay = sim_in.hf_delay_ms*Rs/1000; hf_mag_only = sim_in.hf_mag_only; + f_off = sim_in.f_off; [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, Nsymbrowpilot*Ntrials); @@ -719,8 +734,9 @@ function sim_out = ber_test(sim_in) hf_n = 1; - phase_offset = 0; - w_offset = pi/16; + phase_offset_rect = 1; + w_offset = 2*pi*f_off/Rs; + w_offset_rect = exp(j*w_offset); ct_symb_buf = zeros(2*Nsymbrowpilot, Nc); prev_tx_symb = prev_rx_symb = ones(1,Nc); @@ -750,7 +766,6 @@ function sim_out = ber_test(sim_in) tx_bits_buf(1:framesize) = tx_bits; end - s_ch = tx_symb; % HF channel simulation ------------------------------------ @@ -797,8 +812,12 @@ function sim_out = ber_test(sim_in) noise = sqrt(variance*0.5)*(randn(Nsymbrowpilot,Nc*Nchip) + j*randn(Nsymbrowpilot,Nc*Nchip)); noise_log = [noise_log noise]; - s_ch = s_ch + noise; - + for r=1:Nsymbrowpilot + s_ch(r,:) *= phase_offset_rect; + phase_offset_rect *= w_offset_rect; + end + s_ch += noise; + ct_symb_buf(1:Nsymbrowpilot,:) = ct_symb_buf(Nsymbrowpilot+1:2*Nsymbrowpilot,:); ct_symb_buf(Nsymbrowpilot+1:2*Nsymbrowpilot,:) = s_ch; @@ -887,7 +906,7 @@ function sim_out = ber_test(sim_in) EsNodBSurface(find(EsNodBSurface < -5)) = -5; mesh(x,y,EsNodBSurface); grid - axis([1 (Nc+1)*Nchip 1 Rs*5 -5 15]) + axis([1 (Nc+1)*Nchip 1 Rs*5 -5 25]) title('HF Channel Es/No'); if verbose @@ -907,7 +926,7 @@ function sim_out = ber_test(sim_in) p = Ns; for r=1:m1 if p == Ns - phi_x_counter++; + phi_x_counter += Npilotsframe; p = 0; end p++; @@ -919,6 +938,7 @@ function sim_out = ber_test(sim_in) figure(5); clf subplot(211) + [m n] = size(phi_log); plot(phi_x, phi_log(:,2),'r+;Estimated HF channel phase;') if hf_sim hold on; @@ -928,6 +948,7 @@ function sim_out = ber_test(sim_in) end ylabel('Phase (rads)'); legend('boxoff'); + axis([1 m -1.1*pi 1.1*pi]) subplot(212) plot(phi_x, amp_log(:,2),'r+;Estimated HF channel amp;') @@ -939,17 +960,13 @@ function sim_out = ber_test(sim_in) ylabel('Amplitude'); xlabel('Time (symbols)'); legend('boxoff'); + axis([1 m 0 3]) end figure(4) clf - subplot(211) stem(Nerrs_log) - subplot(212) - if ldpc_code - stem(ldpc_Nerrs_log) - end - + axis([1 length(Nerrs_log) 0 max(Nerrs_log)+1]) end endfunction @@ -958,6 +975,7 @@ endfunction function sim_in = standard_init sim_in.verbose = 1; + sim_in.do_write_pilot_file = 0; sim_in.plot_scatter = 0; sim_in.Esvec = 50; diff --git a/codec2-dev/octave/linreg.m b/codec2-dev/octave/linreg.m new file mode 100644 index 00000000..7a44b8b3 --- /dev/null +++ b/codec2-dev/octave/linreg.m @@ -0,0 +1,33 @@ +% linreg.m +% David Rowe April 2015 +% +% Based on: +% http://stackoverflow.com/questions/5083465/fast-efficient-least-squares-fit-algorithm-in-c + +function [m b] = linreg(x,y,n) + sumx = 0.0; % sum of x + sumx2 = 0.0; % sum of x^2 + sumxy = 0.0; % sum of x * y + sumy = 0.0; % sum of y + sumy2 = 0.0; % sum of y**2 + + for i=1:n + sumx += x(i); + sumx2 += x(i)^2; + sumxy += x(i) * y(i); + sumy += y(i); + sumy2 += y(i)^2; + end + + denom = (n * sumx2 - sumx*sumx); + + if denom == 0 + % singular matrix. can't solve the problem. + m = 0; + b = 0; + else + m = (n * sumxy - sumx * sumy) / denom; + b = (sumy * sumx2 - sumx * sumxy) / denom; + end + +endfunction diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 109f125d..d3cb2390 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -78,6 +78,7 @@ acohpsk.ldpc_code_rate = 1; acohpsk.Nc = Nc; acohpsk.Rs = Rs; acohpsk.Ns = 4; +acohpsk.coh_en = 1; acohpsk.Nchip = 1; acohpsk.modulation = 'qpsk'; acohpsk.do_write_pilot_file = 0; diff --git a/codec2-dev/octave/test_cohpsk.m b/codec2-dev/octave/test_cohpsk.m index 50832bce..cb83e688 100644 --- a/codec2-dev/octave/test_cohpsk.m +++ b/codec2-dev/octave/test_cohpsk.m @@ -36,9 +36,8 @@ function test_curves sim_in.verbose = 1; sim_in.plot_scatter = 1; - sim_in.Esvec = 10; + sim_in.Esvec = 20; sim_in.framesize = 32; - sim_in.hf_sim = 1; sim_in.Ntrials = 100; sim_in.Rs = 50; sim_in.Nc = 4; @@ -48,12 +47,17 @@ function test_curves sim_in.modulation = 'qpsk'; sim_in.ldpc_code_rate = 1; sim_in.ldpc_code = 0; + sim_in.coh_en = 1; + + sim_in.hf_sim = 1; + sim_in.hf_mag_only = 0; + sim_in.f_off = 0; sim_qpsk = ber_test(sim_in); % AWGN curves ---------------------------------------------------- - sim_in.Ntrials = 500; + sim_in.Ntrials = 400; sim_in.hf_sim = 0; sim_in.plot_scatter = 0; sim_in.Esvec = 5:10; @@ -65,26 +69,24 @@ function test_curves sim_dqpsk = ber_test(sim_in, 'dqpsk'); sim_in.modulation = 'qpsk'; - sim_in.Ns = 4; - sim_in.Np = 2; sim_qpsk_pilot = ber_test(sim_in, 'qpsk'); % HF curves ---------------------------------------------------- - sim_in.Ntrials = 200; + sim_in.Ntrials = 400; sim_in.hf_sim = 1; sim_in.plot_scatter = 0; sim_in.Esvec = 5:20; - - Ebvec = sim_in.Esvec - 10*log10(2); - BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); - sim_in.modulation = 'dqpsk'; sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); sim_in.modulation = 'qpsk'; - sim_in.Ns = 4; - sim_in.Np = 2; + sim_in.coh_en = 0; + sim_in.hf_mag_only = 1; + sim_qpsk_hf = ber_test(sim_in, 'qpsk'); + + sim_in.coh_en = 1; + sim_in.hf_mag_only = 0; sim_qpsk_pilot_hf = ber_test(sim_in, 'qpsk'); % plot results --------------------------------------------------- @@ -97,16 +99,16 @@ function test_curves semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') semilogy(sim_qpsk_pilot.Ebvec, sim_qpsk_pilot.BERvec,'b;QPSK pilot AWGN;') - %semilogy(sim_qpsk_hf_ideal.Ebvec, sim_qpsk_hf_ideal.BERvec,'b;QPSK HF ideal;') - semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'c;DQPSK HF;') + semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'r;QPSK HF ideal;') semilogy(sim_qpsk_pilot_hf.Ebvec, sim_qpsk_pilot_hf.BERvec,'b;QPSK pilot HF;') + semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'c;DQPSK HF;') hold off; xlabel('Eb/N0') ylabel('BER') grid("minor") - axis([min(Ebvec) max(Ebvec) 1E-3 1]) + axis([min(sim_qpsk_hf.Ebvec) max(sim_qpsk_hf.Ebvec) 1E-3 1]) legend("boxoff"); endfunction @@ -129,13 +131,13 @@ function test_single sim_in.Ntrials = 100; sim_in.Esvec = 10; - sim_in.hf_sim = 0; + sim_in.hf_sim = 1; sim_in.hf_mag_only = 0; sim_in.modulation = 'qpsk'; + sim_in.coh_en = 1; + sim_in.f_off = 0; - sim_in.modulation = 'dqpsk'; - - sim_in.do_write_pilot_file = 0; + %sim_in.modulation = 'dqpsk'; sim_qpsk_hf = ber_test(sim_in); @@ -593,6 +595,7 @@ endfunction % Start simulations --------------------------------------- more off; +%close all; test_curves(); %test_single(); %rate_Fs_tx("tx_zero.raw"); diff --git a/codec2-dev/octave/tlinreg.m b/codec2-dev/octave/tlinreg.m index 81656e68..1cde350d 100644 --- a/codec2-dev/octave/tlinreg.m +++ b/codec2-dev/octave/tlinreg.m @@ -1,33 +1,7 @@ - -1; - -function [m b] = linreg(x,y,n) - sumx = 0.0; % sum of x - sumx2 = 0.0; % sum of x^2 - sumxy = 0.0; % sum of x * y - sumy = 0.0; % sum of y - sumy2 = 0.0; % sum of y**2 - - for i=1:n - sumx += x(i); - sumx2 += x(i)^2; - sumxy += x(i) * y(i); - sumy += y(i); - sumy2 += y(i)^2; - end - - denom = (n * sumx2 - sumx*sumx); - - if denom == 0 - % singular matrix. can't solve the problem. - m = 0; - b = 0; - else - m = (n * sumxy - sumx * sumy) / denom; - b = (sumy * sumx2 - sumx * sumxy) / denom; - end - -endfunction +% tlinreg +% David Rowe April 2015 +% +% Unit test for linear regression x = [1 2 7 8]; y = [ -0.70702 + 0.70708i 0.77314 - 0.63442i -0.98083 + 0.19511i 0.99508 - 0.09799i] .* [1 -1 1 -1];