--- /dev/null
+% test_qpsk.m\r
+% David Rowe Feb 2014\r
+%\r
+% Single sample per symbol QPSK modem simulation, based on code by Bill Cowley\r
+% Generates curves BER versus E/No curves for different modems. Design to\r
+% test perform initial tests on coherent demodulation for HF channels without\r
+% building a full blown modem. Lacks filtering, timing estimation, frame sync.\r
+% \r
+\r
+1;\r
+\r
+% main test function \r
+\r
+function sim_out = ber_test(sim_in, modulation)\r
+\r
+ framesize = sim_in.framesize;\r
+ Ntrials = sim_in.Ntrials;\r
+ Esvec = sim_in.Esvec;\r
+ phase_offset = sim_in.phase_offset;\r
+ phase_est = sim_in.phase_est;\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
+ hf_spread_hz = sim_in.hf_spread_hz;\r
+ hf_delay_samples = sim_in.hf_delay_samples;\r
+\r
+ bps = 2;\r
+ Nsymb = framesize/bps;\r
+ prev_sym_tx = qpsk_mod([0 0]);\r
+ prev_sym_rx = qpsk_mod([0 0]);\r
+ rx_symb_log = [];\r
+\r
+ Np = 5;\r
+ r_delay_line = zeros(1,Np);\r
+ s_delay_line = zeros(1,Np);\r
+ hf_delay_line = zeros(1,Nsymb+hf_delay_samples);\r
+ spread_main_phi = 0;\r
+ spread_delay_phi = 0;\r
+ spread_main_phi_log = [];\r
+\r
+ % convert "spreading" samples from 1kHz carrier at Fs to complex baseband at Rs\r
+\r
+ Fs = 8000; Fc = 1000;\r
+ fspread = fopen("../unittest/sine1k_2Hz_spread.raw","rb");\r
+ spread1k = fread(fspread, "int16")/10000;\r
+\r
+ % down convert to complex baseband\r
+ spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))');\r
+\r
+ % remove -2000 Hz image\r
+ b = fir1(50, 5/Fs);\r
+ spreadlpf = filter(b,1,spreadbb);\r
+\r
+ % decimate to symbol rate\r
+ spread = spreadlpf(1:floor(Fs/Rs):length(spreadlpf));\r
+\r
+ sc = 1;\r
+\r
+ % design root nyquist (root raised cosine) filter and init tx and rx filter states\r
+\r
+ alpha = 0.5; T=1/Fs; Nsym=7; M=Fs/Rs;\r
+ if floor(Fs/Rs) != Fs/Rs\r
+ printf("oversampling ratio must be an integer\n");\r
+ exit;\r
+ end\r
+ hrn = gen_rn_coeffs(alpha, T, Rs, Nsym, M);\r
+ Nfilter = length(hrn);\r
+ tx_filter_memory = zeros(1, Nfilter);\r
+ tx_baseband_log = [];\r
+ rx_filter_memory = zeros(1, Nfilter);\r
+ rx_baseband_log = [];\r
+ s_delay_line_filt = zeros(1,Nsym);\r
+\r
+ for ne = 1:length(Esvec)\r
+ Es = Esvec(ne);\r
+ EsNo = 10^(Es/10);\r
+ \r
+ variance = Fs/(2*Rs*EsNo);\r
+ Terrs = 0; Tbits = 0; Ferrs = 0;\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:Nsymb\r
+ tx_symb = qpsk_mod(tx_bits(2*(i-1)+1:2*i));\r
+ %printf("shift: %f prev_sym: %f ", tx_symb, prev_sym_tx);\r
+ if strcmp(modulation,'dqpsk')\r
+ tx_symb *= prev_sym_tx;\r
+ %printf("tx_symb: %f\n", tx_symb);\r
+ prev_sym_tx = tx_symb;\r
+ end \r
+ s(i) = tx_symb;\r
+ end\r
+\r
+ s_ch = s;\r
+\r
+ % root nyquist filter symbols\r
+\r
+ for k=1:Nsym\r
+\r
+ % tx filter symbols\r
+\r
+ tx_baseband = zeros(1,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_ch(k);\r
+\r
+ for i=1:M\r
+ tx_baseband(i) = M*tx_filter_memory(M:M:Nfilter) * hrn(M-i+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(1,M);\r
+ tx_baseband_log = [tx_baseband_log tx_baseband];\r
+\r
+ noise = sqrt(variance)*( randn(1,M) + j*randn(1,M) );\r
+ rx_baseband = tx_baseband.*exp(j*phase_offset) + noise;\r
+ phase_offset += w_offset;\r
+\r
+ % rx filter symbol\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
+ rx_baseband_log = [rx_baseband_log rx_filt];\r
+\r
+ % delay in tx data to compensate for filtering\r
+\r
+ s_delay_line_filt(1:Nsym-1) = s_delay_line_filt(2:Nsym);\r
+ s_delay_line_filt(Nsym) = s(k);\r
+ tx_bits(2*(k-1)+1:2*k) = qpsk_demod(s_delay_line_filt(1));\r
+ s(k) = s_delay_line_filt(1); % input to phase est later\r
+\r
+ s_ch(k) = rx_filt;\r
+ end\r
+\r
+ % Channel simulation\r
+\r
+ if hf_sim\r
+ s_ch = s_ch.*spread(sc:sc+Nsymb-1)';\r
+ sc += Nsymb;\r
+ end\r
+\r
+ % coherent demod phase estimation and correction\r
+ % need sliding window\r
+ % use known (pilot) symbols\r
+ % start with all symbols pilots, then gradually decimate, e.g. 1 in 5 pilots\r
+\r
+ if phase_est\r
+ for i=1:Nsymb\r
+\r
+ % delay line for phase est window\r
+\r
+ r_delay_line(1:Np-1) = r_delay_line(2:Np);\r
+ r_delay_line(Np) = s_ch(i);\r
+\r
+ % delay in tx data to compensate data for phase est window\r
+\r
+ s_delay_line(1:Np-1) = s_delay_line(2:Np);\r
+ s_delay_line(Np) = s(i);\r
+ tx_bits(2*(i-1)+1:2*i) = qpsk_demod(s_delay_line(floor(Np/2)+1));\r
+\r
+ % estimate phase from known symbols and correct\r
+\r
+ corr = s_delay_line * r_delay_line';\r
+ s_ch(i) = r_delay_line(floor(Np/2)+1).*exp(j*angle(corr));\r
+ end \r
+ %printf("corr: %f angle: %f\n", corr, angle(corr));\r
+ end\r
+\r
+ % de-modulate\r
+\r
+ rx_bits = zeros(1, framesize);\r
+ for i=1:Nsymb\r
+ rx_symb = s_ch(i);\r
+ if strcmp(modulation,'dqpsk')\r
+ tmp = rx_symb;\r
+ rx_symb *= conj(prev_sym_rx/abs(prev_sym_rx));\r
+ prev_sym_rx = tmp;\r
+ end\r
+ rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb);\r
+ rx_symb_log = [rx_symb_log rx_symb];\r
+ end\r
+\r
+ % Measure BER\r
+\r
+ % discard bits from first 2*Nsym symbols as tx and rx filter memories not full\r
+\r
+ if nn == 1\r
+ tx_bits = tx_bits(2*bps*Nsym+1:length(tx_bits));\r
+ rx_bits = rx_bits(2*bps*Nsym+1:length(rx_bits));\r
+ end\r
+\r
+ error_positions = xor( rx_bits, tx_bits );\r
+ Nerrs = sum(error_positions);\r
+ Terrs += Nerrs;\r
+ nerr(nn) = Nerrs;\r
+ \r
+ if Nerrs>0, Ferrs = Ferrs +1; end\r
+ Terrs = Terrs + Nerrs;\r
+ Tbits = Tbits + framesize;\r
+\r
+ end\r
+ \r
+ TERvec(ne) = Terrs;\r
+ FERvec(ne) = Ferrs;\r
+ BERvec(ne) = Terrs/Tbits;\r
+ end\r
+ \r
+ Ebvec = Esvec - 10*log10(bps);\r
+ sim_out.BERvec = BERvec;\r
+ sim_out.Ebvec = Ebvec;\r
+ sim_out.FERvec = FERvec;\r
+ sim_out.TERvec = TERvec;\r
+\r
+ if plot_scatter\r
+ figure(2);\r
+ scat = rx_symb_log(2*Nsym:length(rx_symb_log)) .* exp(j*pi/4);\r
+\r
+ plot(real(scat), imag(scat),'+');\r
+ figure(3);\r
+ clf\r
+ plot(real(spread(1:100)));\r
+ hold on\r
+ plot(imag(spread(1:100)),'r')\r
+ hold off;\r
+ figure(4)\r
+ subplot(211)\r
+ plot(imag(tx_baseband_log(1:30*M)));\r
+ subplot(212)\r
+ plot(imag(rx_symb_log(2*Nsym:length(rx_symb_log))));\r
+ end\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
+ 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
+% Start simulation ---------------------------------------\r
+\r
+sim_in.Esvec = 1:10; \r
+sim_in.Ntrials = 10;\r
+sim_in.framesize = 100;\r
+sim_in.phase_offset = 0;\r
+sim_in.phase_est = 0;\r
+sim_in.w_offset = 0;\r
+sim_in.plot_scatter = 0;\r
+sim_in.Rs = 100;\r
+sim_in.hf_sim = 0;\r
+sim_in.hf_spread_hz = 2;\r
+sim_in.hf_delay_samples = 5;\r
+\r
+sim_qpsk = ber_test(sim_in, 'qpsk');\r
+\r
+sim_in.phase_offset = 0;\r
+sim_in.phase_est = 0;\r
+sim_in.w_offset = 0; \r
+%sim_qpsk_coh = ber_test(sim_in, 'qpsk');\r
+\r
+sim_in.phase_offset = 0;\r
+sim_in.phase_est = 1;\r
+sim_in.w_offset = 0; \r
+sim_in.plot_scatter = 1;\r
+sim_in.Esvec = 10;\r
+sim_in.hf_sim = 0;\r
+sim_qpsk_scatter = ber_test(sim_in, 'qpsk');\r
+\r
+figure(1); \r
+clf;\r
+semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec)\r
+hold on;\r
+%semilogy(sim_qpsk_coh.Ebvec, sim_qpsk_coh.BERvec,'r;coherent;')\r
+hold off;\r
+xlabel('Eb/N0')\r
+ylabel('BER')\r
+grid\r
+\r
projects / freetel-svn-tracking.git / commitdiff