--- /dev/null
+% test_dqpsk.m\r
+% David Rowe March 2014\r
+%\r
+% Single sample/symbol DQPSK modem simulation to test modulating modem\r
+% tx power based on speech energy.\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
+ nhfdelay = sim_in.hf_delay_ms*Rs/1000;\r
+ hf_phase_only = sim_in.hf_phase_only;\r
+ hf_mag_only = sim_in.hf_mag_only;\r
+ Nc = sim_in.Nc;\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
+ rate = 1;\r
+\r
+ % Init HF channel model from stored sample files of spreading signal ----------------------------------\r
+\r
+ % convert "spreading" samples from 1kHz carrier at Fs to complex\r
+ % baseband, generated by passing a 1kHz sine wave through PathSim\r
+ % with the ccir-poor model, enabling one path at a time.\r
+ \r
+ Fc = 1000; M = Fs/Rs;\r
+ fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb");\r
+ spread1k = fread(fspread, "int16")/10000;\r
+ fclose(fspread);\r
+ fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb");\r
+ spread1k_2ms = fread(fspread, "int16")/10000;\r
+ fclose(fspread);\r
+\r
+ % down convert to complex baseband\r
+ spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))');\r
+ spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))');\r
+\r
+ % remove -2000 Hz image\r
+ b = fir1(50, 5/Fs);\r
+ spread = filter(b,1,spreadbb);\r
+ spread_2ms = filter(b,1,spreadbb_2ms);\r
+\r
+ % discard first 1000 samples as these were near 0, probably as\r
+ % PathSim states were ramping up\r
+\r
+ spread = spread(1000:length(spread));\r
+ spread_2ms = spread_2ms(1000:length(spread_2ms));\r
+\r
+ % decimate down to Rs\r
+\r
+ spread = spread(1:M:length(spread));\r
+ spread_2ms = spread_2ms(1:M:length(spread_2ms));\r
+\r
+ % Determine "gain" of HF channel model, so we can normalise\r
+ % carrier power during HF channel sim to calibrate SNR. I imagine\r
+ % different implementations of ccir-poor would do this in\r
+ % different ways, leading to different BER results. Oh Well!\r
+\r
+ hf_gain = 1.0/sqrt(var(spread)+var(spread_2ms));\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
+ \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
+\r
+ sim_out.errors_log = [];\r
+ sim_out.ldpc_errors_log = [];\r
+\r
+ for nn = 1: Ntrials\r
+ \r
+ tx_bits = round( rand( 1, framesize*rate ) );\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
+ tx_symb *= prev_sym_tx(k);\r
+ prev_sym_tx(k) = tx_symb;\r
+ s(i+k-1) = tx_symb;\r
+ end\r
+ end\r
+ s_ch = s;\r
+\r
+ % HF channel simulation ------------------------------------\r
+ \r
+ if hf_sim\r
+\r
+ % separation between carriers. Note this is\r
+ % effectively under samples at Rs, I dont think this\r
+ % matters. Equivalent to doing freq shift at Fs, then\r
+ % decimating to Rs.\r
+\r
+ wsep = 2*pi*(1+0.5); % e.g. 75Hz spacing at Rs=50Hz, alpha=0.5 filters\r
+\r
+ if Nsymb/Nc != floor(Nsymb/Nc)\r
+ printf("Error: Nsymb/Nc must be an integrer\n")\r
+ return;\r
+ end\r
+\r
+ % arrange symbols in Nsymb/Nc by Nc matrix\r
+\r
+ for i=1:Nc:Nsymb\r
+\r
+ % Determine HF channel at each carrier for this symbol\r
+\r
+ for k=1:Nc\r
+ hf_model(hf_n, k) = hf_gain*(spread(hf_n) + exp(-j*k*wsep*nhfdelay)*spread_2ms(hf_n));\r
+ hf_fading(i+k-1) = abs(hf_model(hf_n, k));\r
+ if hf_mag_only\r
+ s_ch(i+k-1) *= abs(hf_model(hf_n, k));\r
+ else\r
+ s_ch(i+k-1) *= hf_model(hf_n, k);\r
+ end\r
+ end\r
+ hf_n++;\r
+ end\r
+ end\r
+ \r
+ tx_symb_log = [tx_symb_log s_ch];\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
+ error_positions = xor(rx_bits, tx_bits);\r
+ Nerrs = sum(error_positions);\r
+ Terrs += Nerrs;\r
+ Tbits += length(tx_bits);\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 .* 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
+\r
+ figure(4);\r
+ clf;\r
+ subplot(211)\r
+ plot(y,abs(hf_model(y,1)))\r
+ title('HF Channel Carrier 1 Mag');\r
+ subplot(212)\r
+ plot(y,angle(hf_model(y,1)))\r
+ title('HF Channel Carrier 1 Phase');\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
+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
projects / freetel-svn-tracking.git / commitdiff