--- /dev/null
+% test_dqpsk2.m\r
+% David Rowe April 2014\r
+%\r
+% DQPSK modem simulation inclduing filtering to test modulating modem\r
+% tx power based on speech energy. Unlike test_dpsk runs at sample\r
+% rate Fs.\r
+\r
+1;\r
+\r
+% main test function \r
+\r
+function sim_out = ber_test(sim_in)\r
+ Fs = 8000;\r
+\r
+ verbose = sim_in.verbose;\r
+ framesize = sim_in.framesize;\r
+ Ntrials = sim_in.Ntrials;\r
+ Esvec = sim_in.Esvec;\r
+ phase_offset = sim_in.phase_offset;\r
+ w_offset = sim_in.w_offset;\r
+ plot_scatter = sim_in.plot_scatter;\r
+ Rs = sim_in.Rs;\r
+ hf_sim = sim_in.hf_sim;\r
+ Nc = sim_in.Nc;\r
+ symbol_amp = sim_in.symbol_amp;\r
+\r
+ bps = 2;\r
+ Nsymb = framesize/bps;\r
+ for k=1:Nc\r
+ prev_sym_tx(k) = qpsk_mod([0 0]);\r
+ prev_sym_rx(k) = qpsk_mod([0 0]);\r
+ end\r
+\r
+ % design root nyquist (root raised cosine) filter and init tx and rx filter states\r
+\r
+ alpha = 0.5; T=1/Fs; Nfiltsym=7; M=Fs/Rs;\r
+ if floor(Fs/Rs) != Fs/Rs\r
+ printf("oversampling ratio must be an integer\n");\r
+ return;\r
+ end\r
+ hrn = gen_rn_coeffs(alpha, T, Rs, Nfiltsym, M);\r
+ Nfilter = length(hrn);\r
+\r
+ % Start Simulation ----------------------------------------------------------------\r
+\r
+ for ne = 1:length(Esvec)\r
+ EsNodB = Esvec(ne);\r
+ EsNo = 10^(EsNodB/10);\r
+ \r
+ variance = 1/EsNo;\r
+ if verbose > 1\r
+ printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance);\r
+ end\r
+ \r
+ Terrs = 0; Tbits = 0;\r
+\r
+ tx_symb_log = [];\r
+ rx_symb_log = [];\r
+ noise_log = [];\r
+ sim_out.errors_log = [];\r
+ sim_out.tx_baseband_log = [];\r
+ sim_out.rx_filt_log = [];\r
+\r
+ % init HF channel\r
+\r
+ hf_n = 1;\r
+ hf_angle_log = [];\r
+ hf_fading = ones(1,Nsymb); % default input for ldpc dec\r
+ hf_model = ones(Ntrials*Nsymb/Nc, Nc); % defaults for plotting surface\r
+ symbol_amp_index = 1;\r
+\r
+ % bunch of outputs we log for graphing\r
+\r
+ sim_out.errors_log = [];\r
+ sim_out.Nerrs = [];\r
+ sim_out.snr_log = [];\r
+ sim_out.hf_model_pwr = [];\r
+ sim_out.tx_fdm_log = [];\r
+\r
+ % init filter memories and LOs\r
+\r
+ tx_filter_memory = zeros(Nc, Nfilter);\r
+ rx_filter_memory = zeros(Nc, Nfilter);\r
+ s_delay_line_filt = zeros(Nc, Nfiltsym);\r
+ phase_tx = ones(1,Nc);\r
+ phase_rx = ones(1,Nc);\r
+ Fcentre = 1500; Fsep = (1+alpha)*Rs;\r
+ freq = (2*pi/Fs)*(Fcentre/2 + Fsep*(1:Nc));\r
+ freq = exp(j*freq);\r
+\r
+ for nn = 1: Ntrials\r
+ \r
+ tx_bits = round( rand( 1, framesize ) );\r
+ \r
+ % modulate --------------------------------------------\r
+\r
+ s = zeros(1, Nsymb);\r
+ for i=1:Nc:Nsymb\r
+ for k=1:Nc\r
+ tx_symb = qpsk_mod(tx_bits(2*(i-1+k-1)+1:2*(i+k-1)));\r
+ s_qpsk(i+k-1) = tx_symb;\r
+ tx_symb *= prev_sym_tx(k);\r
+ prev_sym_tx(k) = tx_symb;\r
+ s(i+k-1) = symbol_amp(symbol_amp_index)*tx_symb;\r
+ end\r
+ end\r
+ symbol_amp_index++;\r
+ s_ch = s;\r
+\r
+ for i=1:Nc:Nsymb\r
+\r
+ % Delay tx symbols to match delay due to filters qpsk\r
+ % (rather than dqpsk) symbols used for convenience as\r
+ % it's easy to shift symbols than pair of bits\r
+\r
+ s_delay_line_filt(:,1:Nfiltsym-1) = s_delay_line_filt(:,2:Nfiltsym);\r
+ s_delay_line_filt(:,Nfiltsym) = s_qpsk(i:i+Nc-1);\r
+ s_qpsk(i:i+Nc-1) = s_delay_line_filt(:,1); \r
+ for k=1:Nc\r
+ tx_bits(2*(i-1+k-1)+1:2*(i+k-1)) = qpsk_demod(s_qpsk(i+k-1));\r
+ end\r
+\r
+ % tx filter\r
+\r
+ tx_baseband = zeros(Nc,M);\r
+\r
+ % tx filter each symbol, generate M filtered output samples for each symbol.\r
+ % Efficient polyphase filter techniques used as tx_filter_memory is sparse\r
+\r
+ tx_filter_memory(:,Nfilter) = s(i:i+Nc-1);\r
+\r
+ for k=1:M\r
+ tx_baseband(:,k) = M*tx_filter_memory(:,M:M:Nfilter) * hrn(M-k+1:M:Nfilter)';\r
+ end\r
+ tx_filter_memory(:,1:Nfilter-M) = tx_filter_memory(:,M+1:Nfilter);\r
+ tx_filter_memory(:,Nfilter-M+1:Nfilter) = zeros(Nc,M);\r
+\r
+ sim_out.tx_baseband_log = [sim_out.tx_baseband_log tx_baseband];\r
+\r
+ % upcovert\r
+ \r
+ tx_fdm = zeros(1,M);\r
+\r
+ for c=1:Nc\r
+ for k=1:M\r
+ phase_tx(c) = phase_tx(c) * freq(c);\r
+ tx_fdm(k) = tx_fdm(k) + tx_baseband(c,k)*phase_tx(c);\r
+ end\r
+ end\r
+ \r
+ sim_out.tx_fdm_log = [sim_out.tx_fdm_log tx_fdm];\r
+ \r
+ % HF channel\r
+ \r
+ rx_fdm = tx_fdm;\r
+\r
+ % downconvert\r
+\r
+ rx_baseband = zeros(Nc,M);\r
+ for c=1:Nc\r
+ for k=1:M\r
+ phase_rx(c) = phase_rx(c) * freq(c);\r
+ rx_baseband(c,k) = rx_fdm(k)*phase_rx(c)';\r
+ end\r
+ end\r
+\r
+ % rx filter\r
+\r
+ rx_filter_memory(:,Nfilter-M+1:Nfilter) = rx_baseband;\r
+ rx_filt = rx_filter_memory * hrn';\r
+ rx_filter_memory(:,1:Nfilter-M) = rx_filter_memory(:,1+M:Nfilter);\r
+ sim_out.rx_filt_log = [sim_out.rx_filt_log rx_filt];\r
+\r
+ s_ch(i:i+Nc-1) = rx_filt;\r
+ end\r
+\r
+ % HF channel simulation ------------------------------------\r
+ \r
+ if hf_sim\r
+ end\r
+\r
+ % "genie" SNR estimate \r
+ \r
+ snr = (s_ch*s_ch')/(Nsymb*variance);\r
+ sim_out.snr_log = [sim_out.snr_log snr];\r
+ sim_out.hf_model_pwr = [sim_out.hf_model_pwr mean(hf_fading.^2)];\r
+\r
+ % AWGN noise and phase/freq offset channel simulation\r
+ % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im\r
+\r
+ noise = sqrt(variance*0.5)*(randn(1,Nsymb) + j*randn(1,Nsymb));\r
+ noise_log = [noise_log noise];\r
+ \r
+ % organise into carriers to apply frequency and phase offset\r
+\r
+ for i=1:Nc:Nsymb\r
+ for k=1:Nc\r
+ s_ch(i+k-1) = s_ch(i+k-1)*exp(j*phase_offset) + noise(i+k-1);\r
+ end \r
+ phase_offset += w_offset;\r
+ end\r
+ \r
+ % de-modulate\r
+\r
+ rx_bits = zeros(1, framesize);\r
+ for i=1:Nc:Nsymb\r
+ for k=1:Nc\r
+ rx_symb = s_ch(i+k-1);\r
+ tmp = rx_symb;\r
+ rx_symb *= conj(prev_sym_rx(k)/abs(prev_sym_rx(k)));\r
+ prev_sym_rx(k) = tmp;\r
+ rx_bits((2*(i-1+k-1)+1):(2*(i+k-1))) = qpsk_demod(rx_symb);\r
+ rx_symb_log = [rx_symb_log rx_symb];\r
+ end\r
+ end\r
+\r
+ % ignore data until we have enough frames to fill filter memory\r
+ % then count errors\r
+\r
+ if nn > ceil(Nfiltsym/(Nsymb/Nc))\r
+ error_positions = xor(rx_bits, tx_bits);\r
+ sim_out.errors_log = [sim_out.errors_log error_positions];\r
+ Nerrs = sum(error_positions);\r
+ sim_out.Nerrs = [sim_out.Nerrs Nerrs];\r
+ Terrs += Nerrs;\r
+ Tbits += length(tx_bits);\r
+ end\r
+\r
+ end\r
+\r
+ TERvec(ne) = Terrs;\r
+ BERvec(ne) = Terrs/Tbits;\r
+\r
+ if verbose \r
+ printf("EsNo (dB): %f Terrs: %d BER %f BER theory %f", EsNodB, Terrs,\r
+ Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)));\r
+ printf("\n");\r
+ end\r
+ if verbose > 1\r
+ printf("Terrs: %d BER %f BER theory %f C %f N %f Es %f No %f Es/No %f\n\n", Terrs,\r
+ Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), var(tx_symb_log), var(noise_log),\r
+ var(tx_symb_log), var(noise_log), var(tx_symb_log)/var(noise_log));\r
+ end\r
+ end\r
+ \r
+ Ebvec = Esvec - 10*log10(bps);\r
+\r
+ sim_out.BERvec = BERvec;\r
+ sim_out.Ebvec = Ebvec;\r
+ sim_out.TERvec = TERvec;\r
+\r
+ if plot_scatter\r
+ figure(2);\r
+ clf;\r
+ scat = rx_symb_log(Nfiltsym*Nc:length(rx_symb_log)) .* exp(j*pi/4);\r
+ plot(real(scat), imag(scat),'+');\r
+ title('Scatter plot');\r
+\r
+ figure(3);\r
+ clf; \r
+ y = 1:Rs*2;\r
+ x = 1:Nc;\r
+ EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance);\r
+ mesh(x,y,EsNodBSurface);\r
+ grid\r
+ title('HF Channel Es/No');\r
+ end\r
+\r
+endfunction\r
+\r
+% Gray coded QPSK modulation function\r
+\r
+function symbol = qpsk_mod(two_bits)\r
+ two_bits_decimal = sum(two_bits .* [2 1]); \r
+ switch(two_bits_decimal)\r
+ case (0) symbol = 1;\r
+ case (1) symbol = j;\r
+ case (2) symbol = -j;\r
+ case (3) symbol = -1;\r
+ endswitch\r
+endfunction\r
+\r
+% Gray coded QPSK demodulation function\r
+\r
+function two_bits = qpsk_demod(symbol)\r
+ if isscalar(symbol) == 0\r
+ printf("only works with scalars\n");\r
+ return;\r
+ end\r
+ bit0 = real(symbol*exp(j*pi/4)) < 0;\r
+ bit1 = imag(symbol*exp(j*pi/4)) < 0;\r
+ two_bits = [bit1 bit0];\r
+endfunction\r
+\r
+function sim_in = standard_init\r
+ sim_in.verbose = 1;\r
+ sim_in.plot_scatter = 0;\r
+\r
+ sim_in.Esvec = 5:15; \r
+ sim_in.Ntrials = 100;\r
+ sim_in.framesize = 64;\r
+ sim_in.Rs = 100;\r
+ sim_in.Nc = 8;\r
+\r
+ sim_in.phase_offset = 0;\r
+ sim_in.w_offset = 0;\r
+ sim_in.phase_noise_amp = 0;\r
+\r
+ sim_in.hf_delay_ms = 2;\r
+ sim_in.hf_sim = 0;\r
+ sim_in.hf_phase_only = 0;\r
+ sim_in.hf_mag_only = 0;\r
+endfunction\r
+\r
+function awgn_hf_ber_curves()\r
+ sim_in = standard_init();\r
+\r
+ Ebvec = sim_in.Esvec - 10*log10(2);\r
+ BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10)));\r
+\r
+ dpsk_awgn = ber_test(sim_in);\r
+ sim_in.hf_sim = 1;\r
+ dpsk_hf = ber_test(sim_in);\r
+\r
+ figure(1); \r
+ clf;\r
+ semilogy(Ebvec, BER_theory,'r;QPSK theory;')\r
+ hold on;\r
+ semilogy(dpsk_awgn.Ebvec, dpsk_awgn.BERvec,'g;DQPSK;')\r
+ semilogy(dpsk_hf.Ebvec, dpsk_hf.BERvec,'g;DQPSK HF;')\r
+ hold off;\r
+ xlabel('Eb/N0')\r
+ ylabel('BER')\r
+ grid("minor")\r
+ axis([min(Ebvec) max(Ebvec) 1E-3 1])\r
+end\r
+\r
+sim_in = standard_init();\r
+\r
+% energy file sampled every 10ms\r
+\r
+load ../src/ve9qrp.txt\r
+pdB=10*log10(ve9qrp);\r
+for i=1:length(pdB)\r
+ if pdB(i) < 0\r
+ pdB(i) = 0;\r
+ end\r
+end\r
+\r
+% Down sample to 40ms rate used for 1300 bit/s codec, every 4th sample is transmitted\r
+\r
+pdB = pdB(4:4:length(pdB));\r
+\r
+% Use linear mapping function in dB domain to map to symbol power\r
+\r
+power_map_x = [ 0 20 24 40 50 ];\r
+power_map_y = [-6 -6 0 6 6];\r
+mapped_pdB = interp1(power_map_x, power_map_y, pdB);\r
+\r
+%sim_in.symbol_amp = 10 .^ (mapped_pdB/20);\r
+sim_in.symbol_amp = ones(1,length(pdB));\r
+sim_in.plot_scatter = 1;\r
+sim_in.verbose = 2;\r
+sim_in.hf_sim = 1;\r
+sim_in.Esvec = 100;\r
+sim_in.Ntrials = 100;\r
+\r
+dqpsk_pwr_hf = ber_test(sim_in);\r
+\r
+% note: need way to test that power is aligned with speech\r
+\r
+figure(4)\r
+clf;\r
+plot((1:sim_in.Ntrials)*80*4, pdB(1:sim_in.Ntrials));\r
+hold on;\r
+plot((1:sim_in.Ntrials)*80*4, mapped_pdB(1:sim_in.Ntrials),'r');\r
+hold off;\r
+\r
+figure(5)\r
+clf;\r
+\r
+s = load_raw("../raw/ve9qrp.raw");\r
+M=320; M_on_2 = M/2; % processing delay between input speech and centre of analysis window\r
+plot(M_on_2:(M_on_2-1+sim_in.Ntrials*M),s(1:sim_in.Ntrials*M))\r
+hold on;\r
+plot((1:sim_in.Ntrials)*M, 5000*sim_in.symbol_amp(1:sim_in.Ntrials),'r');\r
+hold off;\r
+axis([1 sim_in.Ntrials*M -3E4 3E4]);\r
+\r
+figure(6)\r
+clf;\r
+plot((1:sim_in.Ntrials)*M, 20*log10(sim_in.symbol_amp(1:sim_in.Ntrials)),'b;Es (dB);');\r
+hold on;\r
+plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.hf_model_pwr),'g;Fading (dB);');\r
+plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.snr_log),'r;Es/No (dB);');\r
+\r
+ber = dqpsk_pwr_hf.Nerrs/sim_in.framesize;\r
+ber_clip = ber;\r
+ber_clip(find(ber > 0.2)) = 0.2;\r
+plot((1:length(ber_clip))*M, -20+100*ber_clip,'k;BER (0-20%);');\r
+hold off;\r
+axis([1 sim_in.Ntrials*M -20 20])\r
+\r
+fep=fopen("dqpsk_errors_pwr.bin","wb"); fwrite(fep, dqpsk_pwr_hf.errors_log, "short"); fclose(fep);\r
+fber=fopen("ber.bin","wb"); fwrite(fber, ber, "float"); fclose(fber);\r
projects / freetel-svn-tracking.git / commitdiff