From 66b11650456434fe3e19bbeab7930a427cc34a47 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Mon, 2 Apr 2012 01:13:45 +0000 Subject: [PATCH] fixed some differences between fdmdv_ut and fdmdv_mod, fdmdv_demod git-svn-id: https://svn.code.sf.net/p/freetel/code@362 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/README_fdmdv.txt | 11 ++ codec2-dev/octave/fdmdv.m | 64 ++++++---- codec2-dev/octave/fdmdv_demod.m | 190 ++++++++++++++++++----------- codec2-dev/octave/fdmdv_mod.m | 33 ++--- codec2-dev/octave/fdmdv_ut.m | 55 +++++---- 5 files changed, 215 insertions(+), 138 deletions(-) diff --git a/codec2-dev/octave/README_fdmdv.txt b/codec2-dev/octave/README_fdmdv.txt index 5d5c0e02..ecf22445 100644 --- a/codec2-dev/octave/README_fdmdv.txt +++ b/codec2-dev/octave/README_fdmdv.txt @@ -3,3 +3,14 @@ Modeling sample clock errors using sox: sox -r 8000 -s -2 mod_dqpsk.raw -s -2 mod_dqpsk_8008hz.raw rate -h 8008 + +TODO + +[ ] Get file based mod and demod working again +[ ] write file based ch simulator +[ ] demod outputs ber (maybe after settling time) +[ ] try to run from shell script +[ ] run a few tests +[ ] start coding in C and repeat tests +[ ] arb bit stream input to Octave mod and demod + diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m index fbecfcfe..af8b8ee1 100644 --- a/codec2-dev/octave/fdmdv.m +++ b/codec2-dev/octave/fdmdv.m @@ -28,7 +28,6 @@ global Fsep = 75; % Separation between carriers (Hz) global Fcentre = 1200; % Centre frequency, Nc/2 below this, N/c above (Hz) global Nt = 5; % number of symbols we estimate timing over global P = 4; % oversample factor used for rx symbol filtering -global pilot_bit = 0; % current phase of symbol used to make pilot carrier % Generate root raised cosine (Root Nyquist) filter --------------- @@ -165,9 +164,10 @@ endfunction % Construct FDM signal by frequency shifting each filtered symbol -% stream +% stream. Returns complex signal so we can apply frequency offsets +% easily. -function [tx_fdm pilot] = fdm_upconvert(tx_filt) +function tx_fdm = fdm_upconvert(tx_filt) global Fs; global M; global Nc; @@ -200,7 +200,7 @@ function [tx_fdm pilot] = fdm_upconvert(tx_filt) c = Nc+1; for i=1:M phase_tx(c) = phase_tx(c) * freq(c); - pilot(i) = sqrt(2)*2*tx_filt(c,i)*phase_tx(c); + pilot(i) = 2*tx_filt(c,i)*phase_tx(c); tx_fdm(i) = tx_fdm(i) + pilot(i); end @@ -209,7 +209,7 @@ function [tx_fdm pilot] = fdm_upconvert(tx_filt) % We return the complex (single sided) signal to make frequency % shifting for the purpose of testing easier - tx_fdm = sqrt(2)*tx_fdm; + tx_fdm = 2*tx_fdm; endfunction @@ -292,7 +292,6 @@ function foff = rx_est_freq_offset(rx_fdm, pilot) global Fcentre; global freq; global freq_rx_pilot; - global phase_rx_pilot; global Npilotbaseband; global Npilotlpf; global Npilotcoeff; @@ -476,14 +475,14 @@ function bits = get_test_bits(nbits) global Ntest_bits; % length of test sequence global current_test_bit; global test_bits; - + for i=1:nbits bits(i) = test_bits(current_test_bit++); if (current_test_bit > Ntest_bits) current_test_bit = 1; endif end - + endfunction @@ -518,12 +517,18 @@ endfunction % Initialise ---------------------------------------------------- -global tx_filter_memory = zeros(Nc+1, Nfilter); -global rx_filter_memory = zeros(Nc+1, Nfilter); +global pilot_bit; +pilot_bit = 0; % current value of pilot bit + +global tx_filter_memory; +tx_filter_memory = zeros(Nc+1, Nfilter); +global rx_filter_memory; +rx_filter_memory = zeros(Nc+1, Nfilter); % phasors used for up and down converters -global freq = zeros(Nc+1,1);; +global freq; +freq = zeros(Nc+1,1); for c=1:Nc/2 carrier_freq = (-Nc/2 - 1 + c)*Fsep + Fcentre; freq(c) = exp(j*2*pi*carrier_freq/Fs); @@ -539,31 +544,44 @@ freq(Nc+1) = exp(j*2*pi*Fcentre/Fs); % This really helped PAPR. We don't need to adjust rx % phase a DPSK takes care of that -%global phase_tx = ones(Nc+1,1); -global phase_tx = exp(j*2*pi*(0:Nc)/(Nc+1)); -global phase_rx = ones(Nc+1,1); +global phase_tx; +%phase_tx = ones(Nc+1,1); +phase_tx = exp(j*2*pi*(0:Nc)/(Nc+1)); +global phase_rx; +phase_rx = ones(Nc+1,1); -% Freq offset estimator states +% Freq offset estimator constants global Npilotcoeff = 30; % number of pilot LPF coeffs global pilot_coeff = fir1(Npilotcoeff, 200/(Fs/2))'; % 200Hz LPF global Npilotbaseband = Npilotcoeff + 4*M; % number of pilot baseband samples reqd for pilot LPF global Npilotlpf = 4*M; % number of samples we DFT pilot over, pilot est window -global pilot_baseband = zeros(1, Npilotbaseband); % pilot baseband samples -global pilot_lpf = zeros(1, Npilotlpf); % LPC pilot samples -global phase_rx_pilot = [1 1]; + +% Freq offset estimator states constants + +global pilot_baseband; +pilot_baseband = zeros(1, Npilotbaseband); % pilot baseband samples +global pilot_lpf +pilot_lpf = zeros(1, Npilotlpf); % LPF pilot samples % Timing estimator states -global rx_filter_mem_timing = zeros(Nc, Nt*P); -global rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming); +global rx_filter_mem_timing; +rx_filter_mem_timing = zeros(Nc, Nt*P); +global rx_baseband_mem_timing; +rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming); -% Test bit stream state variables +% Test bit stream constants global Ntest_bits = Nc*Nb*4; % length of test sequence -global current_test_bit = 1; global test_bits = rand(1,Ntest_bits) > 0.5; -global rx_test_bits_mem = zeros(1,Ntest_bits); + +% Test bit stream state variables + +global current_test_bit = 1; +current_test_bit = 1; +global rx_test_bits_mem; +rx_test_bits_mem = zeros(1,Ntest_bits); % Generate M samples of DBPSK pilot signal for Freq offset estimation diff --git a/codec2-dev/octave/fdmdv_demod.m b/codec2-dev/octave/fdmdv_demod.m index 7aff775b..29cd6dc6 100644 --- a/codec2-dev/octave/fdmdv_demod.m +++ b/codec2-dev/octave/fdmdv_demod.m @@ -9,90 +9,136 @@ function fdmdv_demod(rawfilename, nbits) -fdmdv; % include modem code - -modulation = 'dqpsk'; -Foff_hz = 0; -phase_offset = 1; -freq_offset = exp(j*2*pi*Foff_hz/Fs); - -fin = fopen(rawfilename, "rb"); -rx_fdm = fread(fin, Inf, "short"); -gain = 1000; -rx_fdm /= gain; -if (nargin == 1) - frames = floor(length(rx_fdm)/M); -else + fdmdv; % include modem code + + modulation = 'dqpsk'; + + fin = fopen(rawfilename, "rb"); + rx_fdm_buf = fread(fin, Inf, "short"); + gain = 1000; + rx_fdm_buf /= gain; + if (nargin == 1) + frames = floor(length(rx_fdm_buf)/M); + else frames = nbits/(Nc*Nb); -endif + endif -prev_rx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4); + prev_rx_symbols = ones(Nc+1,1); + foff_phase = 1; -total_bit_errors = 0; -total_bits = 0; + % pilot states, used for copy of pilot at rx -rx_timing_log = []; -rx_symbols_log = []; -rx_phase_log = []; -sync_log = []; + pilot_rx_bit = 0; + pilot_symbol = sqrt(2); + pilot_freq = freq(Nc+1); + pilot_phase = 1; + pilot_filter_mem = zeros(1, Nfilter); -% Main loop ---------------------------------------------------- + % BER stats -for i=1:frames - for n=1:M - phase_offset *= freq_offset; - rx_fdm((i-1)*M+1+n) *= phase_offset; - end - rx_baseband = fdm_downconvert(rx_fdm((i-1)*M+1:i*M)); - rx_filt = rx_filter(rx_baseband); + total_bit_errors = 0; + total_bits = 0; + bit_errors_log = []; + sync_log = []; - [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband); - rx_timing_log = [rx_timing_log rx_timing]; + rx_symbols_log = []; + rx_timing_log = []; + foff_log = []; - %rx_phase = rx_est_phase(rx_symbols); - %rx_phase_log = [rx_phase_log rx_phase]; - %rx_symbols = rx_symbols*exp(j*rx_phase); + % Main loop ---------------------------------------------------- + for f=1:frames + rx_fdm = rx_fdm_buf((f-1)*M+1:f*M); - if strcmp(modulation,'dqpsk') - rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)]; - else - rx_symbols_log = [rx_symbols_log rx_symbols]; - endif - rx_bits = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation); - prev_rx_symbols = rx_symbols; + % frequency offset estimation and correction - [sync bit_errors] = put_test_bits(rx_bits); - sync_log = [sync_log sync]; + [pilot pilot_rx_bit pilot_symbol pilot_filter_mem pilot_phase] = generate_pilot_fdm(pilot_rx_bit, pilot_symbol, pilot_filter_mem, pilot_phase, pilot_freq); + foff = rx_est_freq_offset(rx_fdm, pilot); + foff_log = [ foff_log foff ]; + foff_rect = exp(j*2*pi*foff/Fs); - if sync == 1 - total_bit_errors = total_bit_errors + bit_errors; - total_bits = total_bits + Ntest_bits; - endif + for i=1:M + foff_phase *= foff_rect'; + rx_fdm(i) = rx_fdm(i)*foff_phase; + end + + % baseband processing + + rx_baseband = fdm_downconvert(rx_fdm); + rx_filt = rx_filter(rx_baseband); + + [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband); + rx_timing_log = [rx_timing_log rx_timing]; + + if strcmp(modulation,'dqpsk') + rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)]; + else + rx_symbols_log = [rx_symbols_log rx_symbols]; + endif + [rx_bits sync] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation); + prev_rx_symbols = rx_symbols; + sync_log = [sync_log sync]; + + % count bit errors if we find a test frame + + [test_frame_sync bit_errors] = put_test_bits(rx_bits); + if (test_frame_sync == 1 && f > 15) + total_bit_errors = total_bit_errors + bit_errors; + total_bits = total_bits + Ntest_bits; + bit_errors_log = [bit_errors_log bit_errors]; + else + bit_errors_log = [bit_errors_log -1]; + end end -end - -ber = total_bit_errors/total_bits; -printf("%d bits %d errors Meas BER: %1.4f\n", total_bits, total_bit_errors,ber); - -figure(1) -clf; -[n m] = size(rx_symbols_log); -plot(real(rx_symbols_log(:,20:m)),imag(rx_symbols_log(:,20:m)),'+') -%plot(rx_fdm); -figure(2) -clf; -%subplot(211) -plot(rx_timing_log) -%subplot(212) -%plot(sync_log) -figure(3) -clf; -for i=1:Nc - subplot(4,4,i); - mx = max(abs(rx_symbols_log(i,20:m))); - plot(real(rx_symbols_log(i,20:m)),imag(rx_symbols_log(i,20:m)),'+') - axis([-mx mx -mx mx]); -end + % --------------------------------------------------------------------- + % Print Stats + % --------------------------------------------------------------------- + + ber = total_bit_errors / total_bits; + + printf("%d bits %d errors BER: %1.4f\n",total_bits, total_bit_errors, ber); + + % --------------------------------------------------------------------- + % Plots + % --------------------------------------------------------------------- + + figure(1) + clf; + [n m] = size(rx_symbols_log); + plot(real(rx_symbols_log(1:Nc+1,10:m)),imag(rx_symbols_log(1:Nc+1,10:m)),'+') + axis([-2 2 -2 2]); + title('Scatter Diagram'); + + figure(2) + clf; + subplot(211) + plot(rx_timing_log) + title('timing offset (samples)'); + subplot(212) + plot(foff_log) + title('Freq offset (Hz)'); + + figure(3) + clf; + subplot(211) + plot(rx_fdm_buf); + title('FDM Rx Signal'); + subplot(212) + Nfft=Fs; + S=fft(rx_fdm_buf,Nfft); + SdB=20*log10(abs(S)); + plot(SdB(1:Fs/4)) + title('FDM Tx Spectrum'); + + figure(4) + clf; + subplot(211) + stem(sync_log) + title('BPSK Sync') + subplot(212) + stem(bit_errors_log); + title('Bit Errors for test data') + +endfunction diff --git a/codec2-dev/octave/fdmdv_mod.m b/codec2-dev/octave/fdmdv_mod.m index 589047f6..99b17143 100644 --- a/codec2-dev/octave/fdmdv_mod.m +++ b/codec2-dev/octave/fdmdv_mod.m @@ -9,22 +9,23 @@ function tx_fdm = fdmdv_mod(rawfilename, nbits) -fdmdv; % include modem code + fdmdv; % include modem code -frames = floor(nbits/(Nc*Nb)); -tx_fdm = []; -gain = 1000; % Scale up to 16 bit shorts -prev_tx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4); + frames = floor(nbits/(Nc*Nb)) + tx_fdm = []; + gain = 1000; % Scale up to 16 bit shorts + prev_tx_symbols = ones(Nc+1,1); -for i=1:frames - tx_bits = get_test_bits(Nc*Nb); - tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits,'qpsk'); - prev_tx_symbols = tx_symbols; - tx_baseband = tx_filter(tx_symbols); -tx_fdm = [tx_fdm real(fdm_upconvert(tx_baseband))]; -end + for i=1:frames + tx_bits = get_test_bits(Nc*Nb); + tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits,'dqpsk'); + prev_tx_symbols = tx_symbols; + tx_baseband = tx_filter(tx_symbols); + tx_fdm = [tx_fdm real(fdm_upconvert(tx_baseband))]; + end -tx_fdm *= gain; -fout = fopen(rawfilename,"wb"); -fwrite(fout, tx_fdm, "short"); -fclose(fout); + tx_fdm *= gain; + fout = fopen(rawfilename,"wb"); + fwrite(fout, tx_fdm, "short"); + fclose(fout); +endfunction diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m index ae77ae66..b9de56bd 100644 --- a/codec2-dev/octave/fdmdv_ut.m +++ b/codec2-dev/octave/fdmdv_ut.m @@ -7,17 +7,15 @@ % Version 2 % -clear all; - fdmdv; % load modem code % Simulation Parameters -------------------------------------- -frames = 100; -EbNo_dB = 7.3; -Foff_hz = -100; +frames = 50; +EbNo_dB = 73; +Foff_hz = 0; modulation = 'dqpsk'; -hpa_clip = 100; +hpa_clip = 150; % ------------------------------------------------------------ @@ -32,9 +30,8 @@ total_bits = 0; rx_fdm_log = []; rx_baseband_log = []; rx_bits_offset = zeros(Nc*Nb*2); -prev_tx_symbols = sqrt(2)*ones(Nc+1,1); -prev_tx_symbols(Nc+1) = 1; -prev_rx_symbols = sqrt(2)*ones(Nc+1,1); +prev_tx_symbols = ones(Nc+1,1); +prev_rx_symbols = ones(Nc+1,1); foff_log = []; tx_baseband_log = []; tx_fdm_log = []; @@ -92,7 +89,7 @@ t = 0; % Main loop % --------------------------------------------------------------------- -for i=1:frames +for f=1:frames % ------------------- % Modulator @@ -103,21 +100,13 @@ for i=1:frames prev_tx_symbols = tx_symbols; tx_baseband = tx_filter(tx_symbols); tx_baseband_log = [tx_baseband_log tx_baseband]; - [tx_fdm pilot_tx] = fdm_upconvert(tx_baseband); + tx_fdm = fdm_upconvert(tx_baseband); tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M); % ------------------- % Channel simulation % ------------------- - % HPA non-linearity - - i = find(abs(tx_fdm) > hpa_clip); - tx_fdm(i) = hpa_clip*exp(j*angle(tx_fdm(i))); - tx_fdm_log = [tx_fdm_log tx_fdm]; - - rx_fdm = tx_fdm; - % frequency offset for i=1:M @@ -127,9 +116,18 @@ for i=1:frames Foff = Foff_hz; freq_offset = exp(j*2*pi*Foff/Fs); phase_offset *= freq_offset; - rx_fdm(i) = phase_offset*real(rx_fdm(i)); + tx_fdm(i) = phase_offset*tx_fdm(i); end + tx_fdm = real(tx_fdm); + + % HPA non-linearity + + tx_fdm(find(abs(tx_fdm) > hpa_clip)) = hpa_clip; + tx_fdm_log = [tx_fdm_log tx_fdm]; + + rx_fdm = tx_fdm; + % AWGN noise noise = Ngain*randn(1,M); @@ -137,10 +135,11 @@ for i=1:frames rx_fdm += noise; rx_fdm_log = [rx_fdm_log rx_fdm]; - % delay + % Delay - rx_fdm_delay(1:Ndelay-M) = rx_fdm_delay(M+1:Ndelay); - rx_fdm_delay(Ndelay-M+1:Ndelay) = rx_fdm; + %rx_fdm_delay(1:Ndelay-M) = rx_fdm_delay(M+1:Ndelay); + %rx_fdm_delay(Ndelay-M+1:Ndelay) = rx_fdm; + rx_fdm_delay = rx_fdm; % ------------------- % Demodulator @@ -149,7 +148,7 @@ for i=1:frames % frequency offset estimation and correction [pilot pilot_rx_bit pilot_symbol pilot_filter_mem pilot_phase] = generate_pilot_fdm(pilot_rx_bit, pilot_symbol, pilot_filter_mem, pilot_phase, pilot_freq); - foff = rx_est_freq_offset(rx_fdm, pilot); + foff = rx_est_freq_offset(rx_fdm_delay, pilot); foff_log = [ foff_log foff ]; %foff = 0; foff_rect = exp(j*2*pi*foff/Fs); @@ -182,9 +181,10 @@ for i=1:frames sync_log = [sync_log sync]; % count bit errors if we find a test frame + % Allow 15 frames for filter memories to fill and time est to settle [test_frame_sync bit_errors] = put_test_bits(rx_bits); - if (test_frame_sync == 1) + if ((test_frame_sync == 1) && (f > 15)) total_bit_errors = total_bit_errors + bit_errors; total_bits = total_bits + Ntest_bits; bit_errors_log = [bit_errors_log bit_errors]; @@ -216,7 +216,8 @@ printf("Eb/No (meas): %2.2f (%2.2f) dB %d bits %d errors BER: (%1.4f) PAPR: % figure(1) clf; [n m] = size(rx_symbols_log); -plot(real(rx_symbols_log(1:Nc,20:m)),imag(rx_symbols_log(1:Nc,20:m)),'+') +plot(real(rx_symbols_log(1:Nc+1,10:m)),imag(rx_symbols_log(1:Nc+1,10:m)),'+') +axis([-2 2 -2 2]); title('Scatter Diagram'); figure(2) @@ -238,7 +239,7 @@ Nfft=Fs; S=fft(rx_fdm_log,Nfft); SdB=20*log10(abs(S)); plot(SdB(1:Fs/4)) -title('FDM Tx Spectrum'); +title('FDM Rx Spectrum'); figure(4) clf; -- 2.25.1