From 3d6f3fb753e60d21a7b9e6dba2bc33e7c724b466 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Mon, 10 Apr 2017 22:51:01 +0000 Subject: [PATCH] converted to QPSK git-svn-id: https://svn.code.sf.net/p/freetel/code@3095 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/bpsk_hf_fs.m | 193 +++++++++++++++++++++++++-------- 1 file changed, 145 insertions(+), 48 deletions(-) diff --git a/codec2-dev/octave/bpsk_hf_fs.m b/codec2-dev/octave/bpsk_hf_fs.m index 942c016b..36715881 100644 --- a/codec2-dev/octave/bpsk_hf_fs.m +++ b/codec2-dev/octave/bpsk_hf_fs.m @@ -20,11 +20,37 @@ 1; +% Gray coded QPSK modulation function + +function symbol = qpsk_mod(two_bits) + two_bits_decimal = sum(two_bits .* [2 1]); + switch(two_bits_decimal) + case (0) symbol = 1; + case (1) symbol = j; + case (2) symbol = -j; + case (3) symbol = -1; + endswitch +endfunction + + +% Gray coded QPSK demodulation function + +function two_bits = qpsk_demod(symbol) + bit0 = real(symbol*exp(j*pi/4)) < 0; + bit1 = imag(symbol*exp(j*pi/4)) < 0; + two_bits = [bit1 bit0]; +endfunction + + function sim_out = run_sim(sim_in) - Rs = 100; + Rs = 62.5; Fs = 8000; M = Fs/Rs; + Tcp = 0.004; Ncp = Tcp*M; + foffset = 0; + woffset = 2*pi*foffset/Fs; + bps = sim_in.bps; EbNodB = sim_in.EbNodB; verbose = sim_in.verbose; hf_en = sim_in.hf_en; @@ -32,10 +58,15 @@ function sim_out = run_sim(sim_in) Ns = sim_in.Ns; % step size for pilots Nc = sim_in.Nc; % Number of cols, aka number of carriers - Nbitsperframe = (Ns-1)*Nc; - printf("Nbitsperframe: %d\n", Nbitsperframe); - Nrowsperframe = Nbitsperframe/Nc; - printf("Nrowsperframe: %d\n", Nrowsperframe); + Nbitsperframe = (Ns-1)*Nc*bps; + Nrowsperframe = Nbitsperframe/(Nc*bps); + if verbose + printf("Rs:.........: %4.2f\n", Rs); + printf("M:..........: %d\n", M); + printf("bps:.........: %d\n", bps); + printf("Nbitsperframe: %d\n", Nbitsperframe); + printf("Nrowsperframe: %d\n", Nrowsperframe); + end % Important to define run time in seconds so HF model will evolve the same way % for different pilot insertion rates. So lets work backwards from approx @@ -52,7 +83,7 @@ function sim_out = run_sim(sim_in) Nframes = floor((Nrows-1)/Ns); Nbits = Nframes * Nbitsperframe; % number of payload data bits - Nr = Nbits/Nc; % Number of data rows to get Nbits total + Nr = Nbits/(Nc*bps); % Number of data rows to get Nbits total if verbose printf("Nc.....: %d\n", Nc); @@ -63,7 +94,7 @@ function sim_out = run_sim(sim_in) % double check if Nbits fit neatly into carriers - assert(Nbits/Nc == floor(Nbits/Nc), "Nbits/Nc must be an integer"); + assert(Nbits/(Nc*bps) == floor(Nbits/(Nc*bps)), "Nbits/(Nc*bps) must be an integer"); printf("Nframes: %d\n", Nframes); @@ -71,77 +102,119 @@ function sim_out = run_sim(sim_in) % extra row of pilots at end printf("Nrp....: %d (number of rows including pilots)\n", Nrp); + % set up HF model --------------------------------------------------------------- + + if hf_en + + % some typical values, or replace with user supplied + + dopplerSpreadHz = 1.0; path_delay_ms = 1; + + if isfield(sim_in, "dopplerSpreadHz") + dopplerSpreadHz = sim_in.dopplerSpreadHz; + end + if isfield(sim_in, "path_delay_ms") + path_delay_ms = sim_in.path_delay_ms; + end + printf("Doppler Spread: %3.2f Hz Path Delay: %3.2f ms\n", dopplerSpreadHz, path_delay_ms); + + path_delay_samples = path_delay_ms*Fs/1000; + + randn('seed',1); + spread1 = doppler_spread(dopplerSpreadHz, Fs, sim_in.Nsec*Fs*1.1); + spread2 = doppler_spread(dopplerSpreadHz, Fs, sim_in.Nsec*Fs*1.1); + + % sometimes doppler_spread() doesn't return exactly the number of samples we need + + assert(length(spread1) >= Nrp*M, "not enough doppler spreading samples"); + assert(length(spread2) >= Nrp*M, "not enough doppler spreading samples"); + + hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); + % printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1)); + end + % simulate for each Eb/No point ------------------------------------ for nn=1:length(EbNodB) rand('seed',1); randn('seed',1); - EbNo = 10 .^ (EbNodB(nn)/10); + EbNo = bps * (10 .^ (EbNodB(nn)/10)); variance = 1/(M*EbNo/2); - % generate tx bits and insert pilot rows and border cols - - tx_bits = []; tx_bits_np = []; - for f=1:Nframes - tx_bits = [tx_bits; ones(1,Nc+2)]; - for r=1:Nrowsperframe - arowofbits = rand(1,Nc) > 0.5; - tx_bits = [tx_bits; [1 arowofbits 1]]; - tx_bits_np = [tx_bits_np; arowofbits]; - end - end - tx_bits = [tx_bits; [1 ones(1,Nc) 1]]; % final row of pilots - [nr nc] = size(tx_bits); - assert(nr == Nrp); - % map to BPSK symbols + % generate tx bits - tx_sym = 2*tx_bits - 1; + tx_bits = rand(1,Nbits) > 0.5; -#{ - % upsample to rate Fs using Zero Order Hold (ZOH) + % map to symbols in linear array - tx_sym_os = zeros(Nrp*M, Nc+2); - for r=1:Nrp - for rr=(r-1)*M+1:r*M - tx_sym_os(rr,:) = tx_sym(r,:); + if bps == 1 + tx_sym_lin = 2*tx_bits - 1; + end + if bps == 2 + for s=1:Nbits/bps + tx_sym_lin(s) = qpsk_mod(tx_bits(2*(s-1)+1:2*s)); end end -#} + + % place symbols in multi-carrier frame with pilots and boundary carriers + + tx_sym = []; s = 1; + for f=1:Nframes + aframe = zeros(Nrowsperframe,Nc+2); + aframe(1,:) = 1; + for r=1:Nrowsperframe + arowofsymbols = tx_sym_lin(s:s+Nc-1); + s += Nc; + aframe(r+1,2:Nc+1) = arowofsymbols; + end + tx_sym = [tx_sym; aframe]; + end + tx_sym = [tx_sym; ones(1,Nc+2)]; % final row of pilots + [nr nc] = size(tx_sym); + assert(nr == Nrp); % OFDM up conversion and upsampling to rate Fs w = (0:Nc+1)*2*pi*Rs/Fs; - W = zeros(M,Nc+2); + W = zeros(Nc+2,M); for c=1:Nc+2 - W(:,c) = exp(j*w(c)*(0:M-1)); + W(c,:) = exp(j*w(c)*(0:M-1)); end Nsam = Nrp*M; - tx = zeros(Nrp*M,1); + tx = zeros(1,Nrp*M); for r=1:Nrp for c=1:Nc+2 - acarrier = tx_sym(r,c) * W(:,c); + acarrier = tx_sym(r,c) * W(c,:); tx((r-1)*M+1:r*M) += acarrier/M; end end rx = tx; + + if hf_en + %rx = hf_gain * tx(1:Nsam) .* spread1(1:Nsam); + %rx = hf_gain * [tx(path_delay_samples+1:Nsam) zeros(1,path_delay_samples)] .* spread2(1:Nsam); + rx = hf_gain * [tx(path_delay_samples+1:Nsam) zeros(1,path_delay_samples)]; + end - noise = sqrt(variance)*(0.5*randn(Nrp*M,1) + j*0.5*randn(Nrp*M,1)); - rx += noise; + rx = rx .* exp(j*woffset*(1:Nsam)); + noise = sqrt(variance)*(0.5*randn(1,Nsam) + j*0.5*randn(1,Nsam)); + rx += noise; + % downconvert, downsample and integrate using OFDM rx_sym = zeros(Nrp, Nc+2); for r=1:Nrp for c=1:Nc+2 - acarrier = rx((r-1)*M+1:r*M) .* conj(W(:,c)); + acarrier = rx((r-1)*M+1:r*M) .* conj(W(c,:)); rx_sym(r,c) = sum(acarrier); end end - + % pilot based phase est, we use known tx symbols as pilots ---------- phase_est_pilot_log = 10*ones(Nrp,Nc+2); @@ -180,20 +253,33 @@ function sim_out = run_sim(sim_in) end % r=1:Ns:Nrp-Ns end % c=2:Nc+1 - % remove pilots to give us just data symbols - - rx_np = []; + + % remove pilots to give us just data symbols and demodulate + + rx_bits = []; rx_np = []; for r=1:Nrp if mod(r-1,Ns) != 0 - rx_np = [rx_np; rx_corr(r,2:Nc+1)]; + arowofsymbols = rx_corr(r,2:Nc+1); + rx_np = [rx_np arowofsymbols]; + if bps == 1 + arowofbits = real(arowofsymbols) > 0; + end + if bps == 2 + arowofbits = zeros(1,Nc); + for c=1:Nc + arowofbits((c-1)*2+1:c*2) = qpsk_demod(arowofsymbols(c)); + end + end + rx_bits = [rx_bits arowofbits]; end end + assert(length(rx_bits) == Nbits); + % calculate BER stats as a block, after pilots extracted - rx_bits_np = real(rx_np) > 0; - errors = xor(tx_bits_np, rx_bits_np); - Nerrs = sum(sum(errors)); + errors = xor(tx_bits, rx_bits); + Nerrs = sum(errors); printf("EbNodB: %3.2f BER: %4.3f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits, Nbits, Nerrs); @@ -201,7 +287,7 @@ function sim_out = run_sim(sim_in) figure(1) plot(real(tx)) figure(2) - Tx = abs(fft(tx.*hanning(Nsam))); + Tx = abs(fft(tx.*hanning(Nsam)')); Tx_dB = 20*log10(Tx); dF = Fs/Nsam; plot((1:Nsam)*dF, Tx_dB); @@ -212,18 +298,28 @@ function sim_out = run_sim(sim_in) plot(rx_np,'+'); axis([-2 2 -2 2]); + if hf_en figure(4); clf; - plot(abs(hf_model(:,2:Nc+1))); + subplot(211) + plot(abs(spread1(1:Nsam))); + %hold on; plot(abs(spread2(1:Nsam)),'g'); hold off; + subplot(212) + plot(angle(spread1(1:Nsam))); end + figure(5); clf; plot(phase_est_log(:,2:Nc+1),'+', 'markersize', 10); hold on; plot(phase_est_pilot_log(:,2:Nc+1),'g+', 'markersize', 5); + +#{ if sim_in.hf_en plot(angle(hf_model(:,2:Nc+1))); end +#} + axis([1 Nrp -pi pi]); end @@ -270,6 +366,7 @@ end function run_single + sim_in.bps = 2; sim_in.Nc = 8; sim_in.Ns = 8; sim_in.Nsec = 10; -- 2.25.1