getting correct BER courves with filter on AWGN

git-svn-id: https://svn.code.sf.net/p/freetel/code@1397 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/octave/test_qpsk.m b/codec2-dev/octave/test_qpsk.m
index d6baf3aa193b217f2ae89cca9f4eb38e3113e53e..3ddd5b407bff52517b8b121fd23214e2a88d8c5f 100644 (file)
--- a/codec2-dev/octave/test_qpsk.m
+++ b/codec2-dev/octave/test_qpsk.m
@@ -59,18 +59,18 @@ function sim_out = ber_test(sim_in, modulation)
 \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
+    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
         exit;\r
     end\r
-    hrn = gen_rn_coeffs(alpha, T, Rs, Nsym, M);\r
+    hrn = gen_rn_coeffs(alpha, T, Rs, Nfiltsym, 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
+    s_delay_line_filt = zeros(1,Nfiltsym);\r
+ \r
+    wc = 2*pi*1500/Fs;\r
 \r
     for ne = 1:length(Esvec)\r
         Es = Esvec(ne);\r
@@ -78,6 +78,15 @@ function sim_out = ber_test(sim_in, modulation)
     \r
         variance = Fs/(2*Rs*EsNo);\r
         Terrs = 0;  Tbits = 0;  Ferrs = 0;\r
+        printf("EsNo (dB): %f EsNo: %f variance: %f\n", Es, EsNo, variance);\r
+\r
+        tx_filt_log = [];\r
+        rx_filt_log = [];\r
+        rx_baseband_log = [];\r
+        tx_baseband_log = [];\r
+        noise_log = [];\r
+\r
+        tx_phase = rx_phase = 0;\r
 \r
         for nn = 1: Ntrials\r
                   \r
@@ -101,11 +110,11 @@ function sim_out = ber_test(sim_in, modulation)
 \r
             % root nyquist filter symbols\r
 \r
-            for k=1:Nsym\r
+            for k=1:Nsymb\r
 \r
                % tx filter symbols\r
 \r
-               tx_baseband = zeros(1,M);\r
+               tx_filt = 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
@@ -113,33 +122,40 @@ function sim_out = ber_test(sim_in, modulation)
                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
+                   tx_filt(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
+               tx_filt_log = [tx_filt_log tx_filt];\r
+               \r
+               % AWGN noise and phase/freq offset channel simulation\r
 \r
                noise = sqrt(variance)*( randn(1,M) + j*randn(1,M) );\r
-               rx_baseband = tx_baseband.*exp(j*phase_offset) + noise;\r
+               noise_log = [noise_log noise];\r
+               rx_baseband = tx_filt.*exp(j*phase_offset) + noise;\r
                phase_offset += w_offset;\r
-\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
+               rx_filt_log = [rx_filt_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
+               s_delay_line_filt(1:Nfiltsym-1) = s_delay_line_filt(2:Nfiltsym);\r
+               s_delay_line_filt(Nfiltsym) = 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
+               s_ch(k) = rx_filt;               \r
             end\r
 \r
+            %noise = sqrt(variance)*( randn(1,Nsymb) + j*randn(1,Nsymb) );\r
+            %s_ch = s_ch.*exp(j*phase_offset) + noise;\r
+            %phase_offset += w_offset;\r
+\r
             % Channel simulation\r
 \r
             if hf_sim\r
@@ -190,51 +206,40 @@ function sim_out = ber_test(sim_in, modulation)
 \r
             % Measure BER\r
 \r
-            % discard bits from first 2*Nsym symbols as tx and rx filter memories not full\r
+            % discard bits from first 2*Nfiltsym 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
+                tx_bits = tx_bits(2*bps*Nfiltsym+1:length(tx_bits));\r
+                rx_bits = rx_bits(2*bps*Nfiltsym+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
+            Tbits = Tbits + length(tx_bits);\r
 \r
         end\r
     \r
         TERvec(ne) = Terrs;\r
         FERvec(ne) = Ferrs;\r
         BERvec(ne) = Terrs/Tbits;\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_filt_log), var(noise_log),\r
+               var(tx_filt_log)/Rs, var(noise_log)/Fs, (var(tx_filt_log)/Rs)/(var(noise_log)/Fs));\r
     end\r
     \r
     Ebvec = Esvec - 10*log10(bps);\r
     sim_out.BERvec = BERvec;\r
+    sim_out.BER_theoryvec = 0.5*erfc(sqrt(10.^(Ebvec/10)));\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
+        clf;\r
+        scat = rx_symb_log(2*Nfiltsym:length(rx_symb_log)) .* exp(j*pi/4);\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
@@ -261,8 +266,8 @@ endfunction
 % Start simulation ---------------------------------------\r
 \r
 sim_in.Esvec            = 1:10; \r
-sim_in.Ntrials          = 10;\r
-sim_in.framesize        = 100;\r
+sim_in.Ntrials          = 100;\r
+sim_in.framesize        = 30;\r
 sim_in.phase_offset     = 0;\r
 sim_in.phase_est        = 0;\r
 sim_in.w_offset         = 0;\r
@@ -280,10 +285,10 @@ sim_in.w_offset         = 0;
 %sim_qpsk_coh            = ber_test(sim_in, 'qpsk');\r
 \r
 sim_in.phase_offset     = 0;\r
-sim_in.phase_est        = 1;\r
+sim_in.phase_est        = 0;\r
 sim_in.w_offset         = 0;  \r
 sim_in.plot_scatter     = 1;\r
-sim_in.Esvec            = 10;\r
+sim_in.Esvec            = 7;\r
 sim_in.hf_sim           = 0;\r
 sim_qpsk_scatter        = ber_test(sim_in, 'qpsk');\r
 \r
@@ -291,6 +296,7 @@ figure(1);
 clf;\r
 semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec)\r
 hold on;\r
+semilogy(sim_qpsk.Ebvec, sim_qpsk.BER_theoryvec,'r;coherent;')\r
 %semilogy(sim_qpsk_coh.Ebvec, sim_qpsk_coh.BERvec,'r;coherent;')\r
 hold off;\r
 xlabel('Eb/N0')\r