From 0fca5bad1e9ceea9a89f78b63711c16513f7c205 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Thu, 12 Mar 2015 06:34:52 +0000 Subject: [PATCH] first test of coh pdsk modem passes, tets framework started git-svn-id: https://svn.code.sf.net/p/freetel/code@2069 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/autotest.m | 68 ++ codec2-dev/octave/{test_ldpc.m => cohpsk.m} | 1101 +++++-------------- codec2-dev/octave/tcohpsk.m | 30 + codec2-dev/octave/test_cohpsk.m | 578 ++++++++++ codec2-dev/octave/tfdmdv.m | 65 +- codec2-dev/src/cohpsk.c | 8 +- codec2-dev/unittest/tcohpsk.c | 4 +- 7 files changed, 956 insertions(+), 898 deletions(-) create mode 100644 codec2-dev/octave/autotest.m rename codec2-dev/octave/{test_ldpc.m => cohpsk.m} (60%) create mode 100644 codec2-dev/octave/tcohpsk.m create mode 100644 codec2-dev/octave/test_cohpsk.m diff --git a/codec2-dev/octave/autotest.m b/codec2-dev/octave/autotest.m new file mode 100644 index 00000000..9901535b --- /dev/null +++ b/codec2-dev/octave/autotest.m @@ -0,0 +1,68 @@ +% autotest.m +% David Rowe Mar 2015 +% +% Helper functions to plot output of C verson and difference between Octave and C versions + +1; + +function stem_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec) + global no_plot_list; + + if find(no_plot_list == plotnum) + return; + end + figure(plotnum) + subplot(subplotnum) + stem(sig,'g;Octave version;'); + hold on; + stem(error,'r;Octave - C version (hopefully 0);'); + hold off; + if nargin == 6 + axis(axisvec); + end + title(titlestr); +endfunction + + +function plot_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec) + global no_plot_list; + + if find(no_plot_list == plotnum) + return; + end + + figure(plotnum) + subplot(subplotnum) + plot(sig,'g;Octave version;'); + hold on; + plot(error,'r;Octave - C version (hopefully 0);'); + hold off; + if nargin == 6 + axis(axisvec); + end + title(titlestr); +endfunction + + +function check(a, b, test_name) + global passes; + global fails; + + [m n] = size(a); + printf("%s", test_name); + for i=1:(25-length(test_name)) + printf("."); + end + printf(": "); + + e = sum(abs(a - b))/n; + if e < 1E-3 + printf("OK\n"); + passes++; + else + printf("FAIL\n"); + fails++; + end +endfunction + + diff --git a/codec2-dev/octave/test_ldpc.m b/codec2-dev/octave/cohpsk.m similarity index 60% rename from codec2-dev/octave/test_ldpc.m rename to codec2-dev/octave/cohpsk.m index 603c0096..20f4a5d7 100644 --- a/codec2-dev/octave/test_ldpc.m +++ b/codec2-dev/octave/cohpsk.m @@ -1,225 +1,38 @@ -% test_ldpc.m -% David Rowe Oct 2014 +% cohpsk.m +% David Rowe Mar 2015 % - -% Simulation to test FDM QPSK with pilot based coherent detection, -% DSSS, and rate 1/2 LDPC -% -% TODO -% [X] Nc carriers, 588 bit frames -% [X] FEC -% [X] pilot insertion and removal -% [ ] delay on parity carriers -% [X] pilot based phase est -% [ ] uncoded and coded frame sync -% [X] timing estimation, RN filtering, carrier FDM -% [ ] this file getting too big - refactor -% [ ] robust coarse timing -% [ ] Np knob on GUI -% [ ] experimental tuning on EnNo_, fading[], to optimise LDPC dec perf - -% reqd to make sure we can repeat tests exactly - -rand('state',1); -randn('state',1); -graphics_toolkit ("gnuplot"); +% Coherent PSK modem functions, with support for LDPC and DSSS +% (diversity). 1; -% Symbol rate processing for tx side (modulator) - -function [tx_symb tx_bits prev_sym_tx] = symbol_rate_tx(sim_in, tx_bits, code_param, prev_sym_tx) - ldpc_code = sim_in.ldpc_code; - rate = sim_in.ldpc_code_rate; - framesize = sim_in.framesize; - Nsymbrow = sim_in.Nsymbrow; - Nsymbrowpilot = sim_in.Nsymbrowpilot; - Nc = sim_in.Nc; - Npilotsframe = sim_in.Npilotsframe; - Ns = sim_in.Ns; - Nchip = sim_in.Nchip; - modulation = sim_in.modulation; - pilot = sim_in.pilot; - - if ldpc_code - [tx_bits, tmp] = ldpc_enc(tx_bits, code_param); - end - - % modulate -------------------------------------------- - - % organise symbols into a Nsymbrow rows by Nc cols - % data and parity bits are on separate carriers - - tx_symb = zeros(Nsymbrow,Nc); - - for c=1:Nc - for r=1:Nsymbrow - i = (c-1)*Nsymbrow + r; - tx_symb(r,c) = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); - end - end - %tx_symb = zeros(Nsymbrow,Nc); - - % insert pilots, one every Ns data symbols - - tx_symb_pilot = zeros(Nsymbrowpilot, Nc); - - for p=1:Npilotsframe - tx_symb_pilot((p-1)*(Ns+1)+1,:) = pilot(p,:); % row of pilots - %printf("%d %d %d %d\n", (p-1)*(Ns+1)+2, p*(Ns+1), (p-1)*Ns+1, p*Ns); - tx_symb_pilot((p-1)*(Ns+1)+2:p*(Ns+1),:) = tx_symb((p-1)*Ns+1:p*Ns,:); % payload symbols - end - tx_symb = tx_symb_pilot; - - % Optionally copy to other carriers (spreading) - - for c=Nc+1:Nc:Nc*Nchip - tx_symb(:,c:c+Nc-1) = tx_symb(:,1:Nc); - end - - % Optionally DQPSK encode - - if strcmp(modulation,'dqpsk') - for c=1:Nc*Nchip - for r=1:Nsymbrowpilot - tx_symb(r,c) *= prev_sym_tx(c); - prev_sym_tx(c) = tx_symb(r,c); - end - end - end - - % ensures energy/symbol is normalised when spreading - - tx_symb = tx_symb/sqrt(Nchip); -end - - -% compression, John Gibbs pointed out it's best to perform -% non-linear operations on an oversampled signals as they -% tend to generate broadband noise that will be aliased into -% passband if bandwidth is too low - -function y = compress(x, power) - - % oversample by a factor of M - - M = 4; - Ntap = 47; - n = length(x); - - b = fir1(Ntap,1/M); - xM = zeros(1,M*n); - for i=1:n - xM(i*M) = M*x(i); - end - - xM = filter(b,1,xM); - - % non linearity - - yM = sign(xM).*(abs(xM) .^ power); - - % decimate by a factor of M - - yM = filter(b,1,yM); - y = yM(1:M:n*M); - -endfunction - - -% Init HF channel model from stored sample files of spreading signal ---------------------------------- - -function [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, nsam) - - % convert "spreading" samples from 1kHz carrier at Fs to complex - % baseband, generated by passing a 1kHz sine wave through PathSim - % with the ccir-poor model, enabling one path at a time. - - Fc = 1000; M = Fs/Rs; - fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb"); - spread1k = fread(fspread, "int16")/10000; - fclose(fspread); - fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb"); - spread1k_2ms = fread(fspread, "int16")/10000; - fclose(fspread); - - % down convert to complex baseband - spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))'); - spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))'); - - % remove -2000 Hz image - b = fir1(50, 5/Fs); - spread = filter(b,1,spreadbb); - spread_2ms = filter(b,1,spreadbb_2ms); - - % discard first 1000 samples as these were near 0, probably as - % PathSim states were ramping up - - spread = spread(1000:length(spread)); - spread_2ms = spread_2ms(1000:length(spread_2ms)); - - % decimate down to Rs - - spread = spread(1:M:length(spread)); - spread_2ms = spread_2ms(1:M:length(spread_2ms)); - - % Determine "gain" of HF channel model, so we can normalise - % carrier power during HF channel sim to calibrate SNR. I imagine - % different implementations of ccir-poor would do this in - % different ways, leading to different BER results. Oh Well! +% Gray coded QPSK modulation function - hf_gain = 1.0/sqrt(var(spread(1:nsam))+var(spread_2ms(1:nsam))); +function symbol = qpsk_mod(two_bits) + two_bits_decimal = sum(two_bits .* [2 1]); + switch(two_bits_decimal) + case (0) symbol = 1; + case (1) symbol = j; + case (2) symbol = -j; + case (3) symbol = -1; + endswitch endfunction -function write_pilot_file(pilot, Nsymbrowpilot, Ns, Nsymrow, Npilotsframe, Nc); - - filename = sprintf("../src/cohpsk_defs.h", Npilotsframe, Nc); - f=fopen(filename,"wt"); - fprintf(f,"/* Generated by write_pilot_file() Octave function */\n\n"); - fprintf(f,"#define NSYMROW %d /* number of data symbols for each row (i.e. each carrier) */\n", Nsymrow); - fprintf(f,"#define NS %d /* number of data symbols between pilots */\n", Ns); - fprintf(f,"#define NPILOTSFRAME %d /* number of pilot symbols of each row */\n", Npilotsframe); - fprintf(f,"#define PILOTS_NC %d /* number of carriers in coh modem */\n\n", Nc); - fprintf(f,"#define NSYMROWPILOT %d /* length of row after pilots inserted */\n\n", Nsymbrowpilot); - fclose(f); +% Gray coded QPSK demodulation function - filename = sprintf("../src/pilots_coh.h", Npilotsframe, Nc); - f=fopen(filename,"wt"); - fprintf(f,"/* Generated by write_pilot_file() Octave function */\n\n"); - fprintf(f,"float pilots_coh[][PILOTS_NC]={\n"); - for r=1:Npilotsframe - fprintf(f, " {"); - for c=1:Nc-1 - fprintf(f, " %f,", pilot(r, c)); - end - if r < Npilotsframe - fprintf(f, " %f},\n", pilot(r, Nc)); - else - fprintf(f, " %f}\n};", pilot(r, Nc)); +function two_bits = qpsk_demod(symbol) + if isscalar(symbol) == 0 + printf("only works with scalars\n"); + return; end - end - fclose(f); -endfunction - - -% Save test bits frame to a text file in the form of a C array - -function test_bits_coh_file(test_bits_coh) - - f=fopen("test_bits_coh.h","wt"); - fprintf(f,"/* Generated by test_bits_coh_file() Octave function */\n\n"); - fprintf(f,"const int test_bits_coh[]={\n"); - for m=1:length(test_bits_coh)-1 - fprintf(f," %d,\n",test_bits_coh(m)); - endfor - fprintf(f," %d\n};\n",test_bits_coh(length(test_bits_coh))); - fclose(f); - + bit0 = real(symbol*exp(j*pi/4)) < 0; + bit1 = imag(symbol*exp(j*pi/4)) < 0; + two_bits = [bit1 bit0]; endfunction -% init function for symbol rate processing +% init function for symbol rate processing -------------------------------------------------------- function sim_in = symbol_rate_init(sim_in) sim_in.Fs = Fs = 8000; @@ -271,7 +84,7 @@ function sim_in = symbol_rate_init(sim_in) % pilot sequence is used for phase and amplitude estimation, and frame sync - pilot = 1 - 2*(rand(Npilotsframe,Nc) > 0.5) + pilot = 1 - 2*(rand(Npilotsframe,Nc) > 0.5); sim_in.pilot = pilot; sim_in.tx_pilot_buf = [pilot; pilot; pilot]; if sim_in.do_write_pilot_file @@ -312,6 +125,76 @@ function sim_in = symbol_rate_init(sim_in) endfunction +% Symbol rate processing for tx side (modulator) ------------------------------------------------------- + +function [tx_symb tx_bits prev_sym_tx] = symbol_rate_tx(sim_in, tx_bits, code_param, prev_sym_tx) + ldpc_code = sim_in.ldpc_code; + rate = sim_in.ldpc_code_rate; + framesize = sim_in.framesize; + Nsymbrow = sim_in.Nsymbrow; + Nsymbrowpilot = sim_in.Nsymbrowpilot; + Nc = sim_in.Nc; + Npilotsframe = sim_in.Npilotsframe; + Ns = sim_in.Ns; + Nchip = sim_in.Nchip; + modulation = sim_in.modulation; + pilot = sim_in.pilot; + + if ldpc_code + [tx_bits, tmp] = ldpc_enc(tx_bits, code_param); + end + + % modulate -------------------------------------------- + + % organise symbols into a Nsymbrow rows by Nc cols + % data and parity bits are on separate carriers + + tx_symb = zeros(Nsymbrow,Nc); + + for c=1:Nc + for r=1:Nsymbrow + i = (c-1)*Nsymbrow + r; + tx_symb(r,c) = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); + end + end + %tx_symb = zeros(Nsymbrow,Nc); + + % insert pilots, one every Ns data symbols + + tx_symb_pilot = zeros(Nsymbrowpilot, Nc); + + for p=1:Npilotsframe + tx_symb_pilot((p-1)*(Ns+1)+1,:) = pilot(p,:); % row of pilots + %printf("%d %d %d %d\n", (p-1)*(Ns+1)+2, p*(Ns+1), (p-1)*Ns+1, p*Ns); + tx_symb_pilot((p-1)*(Ns+1)+2:p*(Ns+1),:) = tx_symb((p-1)*Ns+1:p*Ns,:); % payload symbols + end + tx_symb = tx_symb_pilot; + + % Optionally copy to other carriers (spreading) + + for c=Nc+1:Nc:Nc*Nchip + tx_symb(:,c:c+Nc-1) = tx_symb(:,1:Nc); + end + + % Optionally DQPSK encode + + if strcmp(modulation,'dqpsk') + for c=1:Nc*Nchip + for r=1:Nsymbrowpilot + tx_symb(r,c) *= prev_sym_tx(c); + prev_sym_tx(c) = tx_symb(r,c); + end + end + end + + % ensures energy/symbol is normalised when spreading + + tx_symb = tx_symb/sqrt(Nchip); +end + + +% Symbol rate processing for rx side (demodulator) ------------------------------------------------------- + function [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx sim_in] = symbol_rate_rx(sim_in, s_ch, prev_sym_rx) framesize = sim_in.framesize; Nsymb = sim_in.Nsymb; @@ -416,65 +299,192 @@ function [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx amp_(r,1:Nc) = amp_(r,1:Nc) + amp_(r,c:c+Nc-1); end end - - % demodulate stage 2 - - rx_symb_linear = zeros(1,Nsymb); - amp_linear = zeros(1,Nsymb); - rx_bits = zeros(1, framesize); - for c=1:Nc - for r=1:Nsymbrow - i = (c-1)*Nsymbrow + r; - rx_symb_linear(i) = rx_symb(r,c); - amp_linear(i) = amp_(r,c); - rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb(r,c)); - end + + % demodulate stage 2 + + rx_symb_linear = zeros(1,Nsymb); + amp_linear = zeros(1,Nsymb); + rx_bits = zeros(1, framesize); + for c=1:Nc + for r=1:Nsymbrow + i = (c-1)*Nsymbrow + r; + rx_symb_linear(i) = rx_symb(r,c); + amp_linear(i) = amp_(r,c); + rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb(r,c)); + end + end + + % Estimate noise power from demodulated symbols. One method is to + % calculate the distance of each symbol from the average symbol + % position. However this is complicated by fading, which means the + % amplitude of the symbols is constantly changing. + + % Now the scatter diagram in a fading channel is a X shape. The + % noise can be resolved into two components at right angles to + % each other. The component along the the "thickness" of the arms + % is proportional to the noise power and not affected by fading. + + v = zeros(1,Nsymb); + for i=1:Nsymb + s = rx_symb_linear(i); + if abs(real(s)) > abs(imag(s)) + v(i) = imag(s); + else + v(i) = real(s); + end + %printf("s: %f %f v: %f\n", real(s), imag(s), v(i)); + end + + % Note we are only measuring variance in one axis, as other axis is obscured by fading. We assume + % that total noise power is shared between both axis so multiply by sqrt(2) to get an estimate of + % total noise pwr. Small constant prevents divide by zero errors on start up. + + No_ = var(v)*sqrt(2) + 1E-6; + + % Estimate signal power + + Es_ = mean(amp_linear .^ 2); + + EsNo_ = Es_/No_; + %printf("Es_: %f No_: %f Es/No: %f Es/No dB: %f\n", Es_, No_, Es_/No_, 10*log10(EsNo_)); + + % LDPC decoder requires some amplitude normalisation + % (AGC), was found to break ow. So we adjust the symbol + % amplitudes so that they are an averge of 1 + + rx_symb_linear /= mean(amp_linear); + amp_linear /= mean(amp_linear); + +endfunction + + +% Compression, John Gibbs pointed out it's best to perform non-linear +% operations on an oversampled signals as they tend to generate +% broadband noise that will be aliased into passband if bandwidth is +% too low + +function y = compress(x, power) + + % oversample by a factor of M + + M = 4; + Ntap = 47; + n = length(x); + + b = fir1(Ntap,1/M); + xM = zeros(1,M*n); + for i=1:n + xM(i*M) = M*x(i); + end + + xM = filter(b,1,xM); + + % non linearity + + yM = sign(xM).*(abs(xM) .^ power); + + % decimate by a factor of M + + yM = filter(b,1,yM); + y = yM(1:M:n*M); + +endfunction + + +% Init HF channel model from stored sample files of spreading signal ---------------------------------- + +function [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, nsam) + + % convert "spreading" samples from 1kHz carrier at Fs to complex + % baseband, generated by passing a 1kHz sine wave through PathSim + % with the ccir-poor model, enabling one path at a time. + + Fc = 1000; M = Fs/Rs; + fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb"); + spread1k = fread(fspread, "int16")/10000; + fclose(fspread); + fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb"); + spread1k_2ms = fread(fspread, "int16")/10000; + fclose(fspread); + + % down convert to complex baseband + spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))'); + spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))'); + + % remove -2000 Hz image + b = fir1(50, 5/Fs); + spread = filter(b,1,spreadbb); + spread_2ms = filter(b,1,spreadbb_2ms); + + % discard first 1000 samples as these were near 0, probably as + % PathSim states were ramping up + + spread = spread(1000:length(spread)); + spread_2ms = spread_2ms(1000:length(spread_2ms)); + + % decimate down to Rs + + spread = spread(1:M:length(spread)); + spread_2ms = spread_2ms(1:M:length(spread_2ms)); + + % Determine "gain" of HF channel model, so we can normalise + % carrier power during HF channel sim to calibrate SNR. I imagine + % different implementations of ccir-poor would do this in + % different ways, leading to different BER results. Oh Well! + + hf_gain = 1.0/sqrt(var(spread(1:nsam))+var(spread_2ms(1:nsam))); +endfunction + + +function write_pilot_file(pilot, Nsymbrowpilot, Ns, Nsymrow, Npilotsframe, Nc); + + filename = sprintf("../src/cohpsk_defs.h", Npilotsframe, Nc); + f=fopen(filename,"wt"); + fprintf(f,"/* Generated by write_pilot_file() Octave function */\n\n"); + fprintf(f,"#define NSYMROW %d /* number of data symbols for each row (i.e. each carrier) */\n", Nsymrow); + fprintf(f,"#define NS %d /* number of data symbols between pilots */\n", Ns); + fprintf(f,"#define NPILOTSFRAME %d /* number of pilot symbols of each row */\n", Npilotsframe); + fprintf(f,"#define PILOTS_NC %d /* number of carriers in coh modem */\n\n", Nc); + fprintf(f,"#define NSYMROWPILOT %d /* length of row after pilots inserted */\n\n", Nsymbrowpilot); + fclose(f); + + filename = sprintf("../src/pilots_coh.h", Npilotsframe, Nc); + f=fopen(filename,"wt"); + fprintf(f,"/* Generated by write_pilot_file() Octave function */\n\n"); + fprintf(f,"float pilots_coh[][PILOTS_NC]={\n"); + for r=1:Npilotsframe + fprintf(f, " {"); + for c=1:Nc-1 + fprintf(f, " %f,", pilot(r, c)); + end + if r < Npilotsframe + fprintf(f, " %f},\n", pilot(r, Nc)); + else + fprintf(f, " %f}\n};", pilot(r, Nc)); end + end + fclose(f); +endfunction - % Estimate noise power from demodulated symbols. One method is to - % calculate the distance of each symbol from the average symbol - % position. However this is complicated by fading, which means the - % amplitude of the symbols is constantly changing. - - % Now the scatter diagram in a fading channel is a X shape. The - % noise can be resolved into two components at right angles to - % each other. The component along the the "thickness" of the arms - % is proportional to the noise power and not affected by fading. - - v = zeros(1,Nsymb); - for i=1:Nsymb - s = rx_symb_linear(i); - if abs(real(s)) > abs(imag(s)) - v(i) = imag(s); - else - v(i) = real(s); - end - %printf("s: %f %f v: %f\n", real(s), imag(s), v(i)); - end - % Note we are only measuring variance in one axis, as other axis is obscured by fading. We assume - % that total noise power is shared between both axis so multiply by sqrt(2) to get an estimate of - % total noise pwr. Small constant prevents divide by zero errors on start up. +% Save test bits frame to a text file in the form of a C array - No_ = var(v)*sqrt(2) + 1E-6; +function test_bits_coh_file(test_bits_coh) - % Estimate signal power - - Es_ = mean(amp_linear .^ 2); - - EsNo_ = Es_/No_; - %printf("Es_: %f No_: %f Es/No: %f Es/No dB: %f\n", Es_, No_, Es_/No_, 10*log10(EsNo_)); - - % LDPC decoder requires some amplitude normalisation - % (AGC), was found to break ow. So we adjust the symbol - % amplitudes so that they are an averge of 1 + f=fopen("test_bits_coh.h","wt"); + fprintf(f,"/* Generated by test_bits_coh_file() Octave function */\n\n"); + fprintf(f,"const int test_bits_coh[]={\n"); + for m=1:length(test_bits_coh)-1 + fprintf(f," %d,\n",test_bits_coh(m)); + endfor + fprintf(f," %d\n};\n",test_bits_coh(length(test_bits_coh))); + fclose(f); - rx_symb_linear /= mean(amp_linear); - amp_linear /= mean(amp_linear); - endfunction +% Rate Rs BER tests ------------------------------------------------------------------------------ + function sim_out = ber_test(sim_in) sim_in = symbol_rate_init(sim_in); @@ -751,29 +761,7 @@ function sim_out = ber_test(sim_in) endfunction -% Gray coded QPSK modulation function - -function symbol = qpsk_mod(two_bits) - two_bits_decimal = sum(two_bits .* [2 1]); - switch(two_bits_decimal) - case (0) symbol = 1; - case (1) symbol = j; - case (2) symbol = -j; - case (3) symbol = -1; - endswitch -endfunction - -% Gray coded QPSK demodulation function -function two_bits = qpsk_demod(symbol) - if isscalar(symbol) == 0 - printf("only works with scalars\n"); - return; - end - bit0 = real(symbol*exp(j*pi/4)) < 0; - bit1 = imag(symbol*exp(j*pi/4)) < 0; - two_bits = [bit1 bit0]; -endfunction function sim_in = standard_init sim_in.verbose = 1; @@ -794,552 +782,3 @@ function sim_in = standard_init sim_in.Nchip = 1; endfunction - -function test_curves - - sim_in = standard_init(); - - sim_in.verbose = 1; - sim_in.plot_scatter = 1; - - sim_in.Esvec = 10; - sim_in.hf_sim = 1; - sim_in.Ntrials = 1000; - sim_in.Rs = 50; - sim_in.Nc = 9; - sim_in.Np = 4; - sim_in.Ns = 8; - sim_in.Nchip = 1; - sim_in.modulation = 'qpsk'; - sim_in.ldpc_code_rate = 0.5; - sim_in.ldpc_code = 0; - - sim_qpsk = ber_test(sim_in); - - sim_in.hf_sim = 0; - sim_in.plot_scatter = 0; - sim_in.Esvec = 10:20; - Ebvec = sim_in.Esvec - 10*log10(2); - BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); - - sim_in.Np = 0; - sim_in.Nchip = 1; - - sim_in.modulation = 'dqpsk'; - sim_dqpsk = ber_test(sim_in, 'dqpsk'); - sim_in.hf_sim = 1; - sim_in.hf_mag_only = 1; - sim_in.modulation = 'qpsk'; - sim_qpsk_hf_ideal = ber_test(sim_in, 'qpsk'); - sim_in.modulation = 'dqpsk'; - sim_in.hf_mag_only = 0; - sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); - sim_in.modulation = 'qpsk'; - sim_in.Ns = 4 - sim_in.Np = 2; - sim_qpsk_hf_pilot = ber_test(sim_in, 'qpsk'); - - figure(1); - clf; - semilogy(Ebvec, BER_theory,'r;QPSK theory;') - hold on; - semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') - semilogy(sim_qpsk_hf_ideal.Ebvec, sim_qpsk_hf_ideal.BERvec,'b;QPSK HF ideal;') - semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'k;DQPSK HF;') - semilogy(sim_qpsk_hf_pilot.Ebvec, sim_qpsk_hf_pilot.BERvec,'r;QPSK Np=2 Ns=4 HF;') - hold off; - - xlabel('Eb/N0') - ylabel('BER') - grid("minor") - axis([min(Ebvec) max(Ebvec) 1E-3 1]) -endfunction - -function test_single - - sim_in = standard_init(); - - sim_in.verbose = 1; - sim_in.plot_scatter = 1; - - sim_in.framesize = 160; - sim_in.Nc = 4; - sim_in.Rs = 50; - sim_in.Ns = 4; - sim_in.Np = 2; - sim_in.Nchip = 1; -% sim_in.ldpc_code_rate = 0.5; -% sim_in.ldpc_code = 1; - sim_in.ldpc_code_rate = 1; - sim_in.ldpc_code = 0; - - sim_in.Ntrials = 20; - sim_in.Esvec = 100; - sim_in.hf_sim = 0; - sim_in.hf_mag_only = 0; - sim_in.modulation = 'qpsk'; - - sim_in.do_write_pilot_file = 1; - - sim_qpsk_hf = ber_test(sim_in); - - fep=fopen("errors_450.bin","wb"); fwrite(fep, sim_qpsk_hf.ldpc_errors_log, "short"); fclose(fep); -endfunction - -% Rate Fs test funcs ----------------------------------------------------------- - -function rate_Fs_tx(tx_filename) - sim_in = standard_init(); - - sim_in.verbose = 1; - sim_in.plot_scatter = 1; - - sim_in.framesize = 160; - sim_in.Nc = 4; - sim_in.Rs = 50; - sim_in.Ns = 4; - sim_in.Np = 2; - sim_in.Nchip = 1; - sim_in.ldpc_code_rate = 1; - sim_in.ldpc_code = 0; - - sim_in.Ntrials = 20; - sim_in.Esvec = 7; - sim_in.hf_sim = 1; - sim_in.hf_mag_only = 0; - sim_in.modulation = 'qpsk'; - - sim_in = symbol_rate_init(sim_in); - - prev_sym_tx = sim_in.prev_sym_tx; - code_param = sim_in.code_param; - tx_bits_buf = sim_in.tx_bits_buf; - framesize = sim_in.framesize; - rate = sim_in.ldpc_code_rate; - Ntrials = sim_in.Ntrials; - Rs = sim_in.Rs; - Fs = sim_in.Fs; - Nc = sim_in.Nc; - - M = Fs/Rs; - - EsNodB = sim_in.Esvec(1); - EsNo = 10^(EsNodB/10); - - rx_symb_log = []; av_tx_pwr = []; - - rn_coeff = gen_rn_coeffs(0.5, 1/Fs, Rs, 6, M); - tx_symb_buf = []; - - tx_bits = round(rand(1,framesize*rate)); - - for nn=1:Ntrials+2 - - [tx_symb tx_bits_out prev_sym_tx] = symbol_rate_tx(sim_in, tx_bits, code_param, prev_sym_tx); - tx_bits_buf(1:framesize) = tx_bits_buf(framesize+1:2*framesize); - tx_bits_buf(framesize+1:2*framesize) = tx_bits_out; - tx_symb_buf = [tx_symb_buf; tx_symb]; - end - - % zero pad and tx filter - - [m n] = size(tx_symb_buf); - zp = []; - for i=1:m - zrow = M*tx_symb_buf(i,:); - zp = [zp; zrow; zeros(M-1,Nc)]; - end - - for c=1:Nc - tx_filt(:,c) = filter(rn_coeff, 1, zp(:,c)); - end - - % upconvert to real IF and save to disk - - [m n] = size(tx_filt); - tx_fdm = zeros(1,m); - Fc = 1500; - for c=1:Nc - freq(c) = exp(j*2*pi*(Fc - c*Rs*1.5)/Fs); - end - phase_tx = ones(1,Nc); - %phase_tx = exp(j*2*pi*(0:Nc)/(Nc+1)); - - for c=1:Nc - for i=1:m - phase_tx(c) = phase_tx(c) * freq(c); - tx_fdm(i) = tx_fdm(i) + tx_filt(i,c)*phase_tx(c); - end - end - - tx_fdm = real(tx_fdm); - - %tx_fdm = compress(tx_fdm, 0.4); - %tx_fdm = sign(tx_fdm) .* (abs(tx_fdm) .^ 0.4); - %hpa_clip = max(abs(tx_fdm))*0.8 - %tx_fdm(find(abs(tx_fdm) > hpa_clip)) = hpa_clip; - - papr = max(tx_fdm.*conj(tx_fdm)) / mean(tx_fdm.*conj(tx_fdm)); - papr_dB = 10*log10(papr); - printf("PAPR: %4.2f dB\n", papr_dB); - - Ascale = 2000; - figure(1); - clf; - plot(Ascale*tx_fdm(1:8000)) - - ftx=fopen(tx_filename,"wb"); fwrite(ftx, Ascale*real(tx_fdm), "short"); fclose(ftx); - -endfunction - - -function rate_Fs_rx(rx_filename) - sim_in = standard_init(); - - sim_in.verbose = 1; - sim_in.plot_scatter = 1; - - sim_in.framesize = 160; - sim_in.Nc = 4; - sim_in.Rs = 50; - sim_in.Ns = 4; - sim_in.Np = 4; - sim_in.Nchip = 1; - sim_in.ldpc_code_rate = 1; - sim_in.ldpc_code = 0; - - sim_in.Ntrials = 10; - sim_in.Esvec = 40; - sim_in.hf_sim = 1; - sim_in.hf_mag_only = 0; - sim_in.modulation = 'qpsk'; - - sim_in = symbol_rate_init(sim_in); - - prev_sym_tx = sim_in.prev_sym_tx; - prev_sym_rx = sim_in.prev_sym_rx; - code_param = sim_in.code_param; - tx_bits_buf = sim_in.tx_bits_buf; - framesize = sim_in.framesize; - rate = sim_in.ldpc_code_rate; - Ntrials = sim_in.Ntrials; - Rs = sim_in.Rs; - Fs = sim_in.Fs; - Nc = sim_in.Nc; - Nsymbrowpilot = sim_in.Nsymbrowpilot; - pilot = sim_in.pilot; - Ns = sim_in.Ns; - Npilotsframe = sim_in.Npilotsframe; - - M = Fs/Rs; - - EsNodB = sim_in.Esvec(1); - EsNo = 10^(EsNodB/10); - - phi_log = []; amp_log = []; EsNo__log = []; - rx_symb_log = []; av_tx_pwr = []; - Terrs = Tbits = 0; - errors_log = []; Nerrs_log = []; - ldpc_Nerrs_log = []; ldpc_errors_log = []; - Ferrsldpc = Terrsldpc = Tbitsldpc = 0; - - rn_coeff = gen_rn_coeffs(0.5, 1/Fs, Rs, 6, M); - - tx_bits = round(rand(1,framesize*rate)); - - % read from disk - - Ascale = 2000; - frx=fopen(rx_filename,"rb"); rx_fdm = fread(frx, "short")/Ascale; fclose(frx); - - rx_fdm=sqrt(2)*rx_fdm; - - if 0 - % optionally add AWGN noise for testing calibration of Es//No measurement - - snr_3000Hz_dB = -9; - snr_4000Hz_lin = 10 ^ (snr_3000Hz_dB/10); - snr_4000Hz_lin *= (3000/4000); - variance = var(rx_fdm)/snr_4000Hz_lin; - rx_fdm += sqrt(variance)*randn(length(rx_fdm),1); - end - - figure(1) - plot(rx_fdm(800:1200)); - - % freq offset estimation - - printf("Freq offset and coarse timing est...\n"); - [f_max max_s_Fs] = test_freq_off_est(rx_filename, 1,5*6400); - f_max = 0; max_s_Fs = 4; - max_s = floor(max_s_Fs/M + 6); - - printf("Downconverting...\n"); - - [m n] = size(rx_fdm); - rx_symb = zeros(m,Nc); - Fc = 1500; - for c=1:Nc - freq(c) = exp(-j*2*pi*(f_max + Fc - c*Rs*1.5)/Fs); - end - phase_rx = ones(1,Nc); - rx_bb = zeros(m,Nc); - - for c=1:Nc - for i=1:m - phase_rx(c) = phase_rx(c) * freq(c); - rx_bb(i,c) = rx_fdm(i)*phase_rx(c); - end - end - - printf("Filtering...\n"); - for c=1:Nc - rx_filt(:,c) = filter(rn_coeff, 1, rx_bb(:,c)); - end - - %subplot(211); - %plot(real(rx_filt(1:10*M,9))); - %subplot(212); - %plot(imag(rx_filt(1:10*M,9))); - - % Fine timing estimation and decimation to symbol rate Rs. Break rx - % signal into ft sample blocks. If clock offset is 1000ppm, - % that's one more/less sample over Ft samples at Fs=8000 Hz. - - printf("Fine timing estimation....\n"); - ft = M*10; - [nsam m] = size(rx_filt); - rx_symb_buf = []; rx_timing_log = []; - - for st=1:ft:floor(nsam/ft - 1)*ft - % fine timing and decimation - - env = zeros(ft,1); - for c=1:Nc - env = env + abs(rx_filt(st:st+ft-1,c)); - end - - % The envelope has a frequency component at the symbol rate. The - % phase of this frequency component indicates the timing. So work out - % single DFT at frequency 2*pi/M - - x = exp(-j*2*pi*(0:ft-1)/M) * env; - - norm_rx_timing = angle(x)/(2*pi); - %norm_rx_timing = -0.4; - if norm_rx_timing < 0 - rx_timing = -floor(norm_rx_timing*M+0.5) + M; - else - rx_timing = -floor(norm_rx_timing*M+0.5) + 2*M; - end - rx_timing_log = [rx_timing_log norm_rx_timing]; - - % printf("%d %d\n", st+rx_timing, st+rx_timing+ft-1); - rx_symb_buf = [rx_symb_buf; rx_filt(st+rx_timing:M:st+rx_timing+ft-1,:)]; - end - - figure(2) - clf; - plot(rx_timing_log) - axis([1 length(rx_timing_log) -0.5 0.5 ]) - title('fine timing') - - printf("Symbol rate demodulation....\n"); - phase_off = 1; - Ntrials = floor((nsam/M)/Nsymbrowpilot) - 2; - %max_s = 6; - - for nn=1:Ntrials - - s_ch = rx_symb_buf((nn-1)*Nsymbrowpilot+max_s:nn*Nsymbrowpilot+max_s-1,:); - [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx sim_in] = symbol_rate_rx(sim_in, s_ch, prev_sym_rx); - - rx_symb_log = [rx_symb_log rx_symb_linear]; - phi_log = [phi_log; phi_]; - amp_log = [amp_log; amp_]; - - if nn > 1 - EsNo__log = [EsNo__log EsNo_]; - - % Measure BER - - error_positions = xor(rx_bits(1:framesize*rate), tx_bits(1:framesize*rate)); - Nerrs = sum(error_positions); - Terrs += Nerrs; - Tbits += framesize*rate; - errors_log = [errors_log error_positions]; - Nerrs_log = [Nerrs_log Nerrs]; - - if sim_in.ldpc_code - % LDPC decode - - detected_data = ldpc_dec(code_param, sim_in.max_iterations, sim_in.demod_type, sim_in.decoder_type, ... - rx_symb_linear, min(100,EsNo_), amp_linear); - error_positions = xor(detected_data(1:framesize*rate), tx_bits(1:framesize*rate) ); - Nerrs = sum(error_positions); - ldpc_Nerrs_log = [ldpc_Nerrs_log Nerrs]; - ldpc_errors_log = [ldpc_errors_log error_positions]; - if Nerrs - Ferrsldpc++; - end - Terrsldpc += Nerrs; - Tbitsldpc += framesize*rate; - end - end - end - - EsNo_av = mean(10*log10(EsNo__log)); - printf("EsNo est (dB): %3.1f SNR est: %3.2f Terrs: %d BER %4.2f QPSK BER theory %4.2f av_tx_pwr: %3.2f", - EsNo_av, mean(EsNo_to_SNR(10*log10(EsNo__log))), - Terrs, - Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), av_tx_pwr); - if sim_in.ldpc_code - printf("\n LDPC: Terrs: %d BER: %4.2f Ferrs: %d FER: %4.2f\n", - Terrsldpc, Terrsldpc/Tbitsldpc, Ferrsldpc, Ferrsldpc/Ntrials); - end - - figure(3); - clf; - scat = rx_symb_log .* exp(j*pi/4); - scat = scat((framesize):length(scat)); - plot(real(scat), imag(scat),'+'); - title('Scatter plot'); - - figure(4) - clf - subplot(211) - stem(Nerrs_log) - axis([0 Ntrials+1 0 max(Nerrs_log)+1]) - title('Uncoded Errors/Frame'); - if sim_in.ldpc_code - subplot(212) - stem(ldpc_Nerrs_log) - axis([0 Ntrials+1 0 max(ldpc_Nerrs_log)+1]) - title('Coded Errors/Frame'); - end - - figure(5) - clf - [m1 n1] = size(phi_log); - phi_x = []; - phi_x_counter = 1; - p = Ns; - for r=1:m1 - if p == Ns - phi_x_counter++; - p = 0; - end - p++; - phi_x = [phi_x phi_x_counter++]; - end - - subplot(211) - plot(phi_x, phi_log(:,1),'r+;Estimated HF channel phase;') - ylabel('Phase (rads)'); - - subplot(212) - plot(phi_x, amp_log(:,1),'r+;Estimated HF channel amp;') - ylabel('Amplitude'); - xlabel('Time (symbols)'); - - figure(6); - clf - plot(EsNo_to_SNR(10*log10(EsNo__log))); - title('SNR est (dB)') - - fep=fopen("errors_450.bin","wb"); fwrite(fep, ldpc_errors_log, "short"); fclose(fep); - -endfunction - -% ideas: cld estimate timing with freq offset and decimate to save cpu load -% fft to do cross correlation - -function [f_max s_max] = test_freq_off_est(rx_filename, offset, n) - fpilot = fopen("tx_zero.raw","rb"); tx_pilot = fread(fpilot, "short"); fclose(fpilot); - frx=fopen(rx_filename,"rb"); rx_fdm = fread(frx, "short"); fclose(frx); - - Fs = 8000; - nc = 1800; % portion we wish to correlate over (first 2 rows on pilots) - - % downconvert to complex baseband to remove images - - f = 1000; - foff_rect = exp(j*2*pi*f*(1:2*n)/Fs); - tx_pilot_bb = tx_pilot(1:n) .* foff_rect(1:n)'; - rx_fdm_bb = rx_fdm(offset:offset+2*n-1) .* foff_rect'; - - % remove -2000 Hz image - - b = fir1(50, 1000/Fs); - tx_pilot_bb_lpf = filter(b,1,tx_pilot_bb); - rx_fdm_bb_lpf = filter(b,1,rx_fdm_bb); - - % decimate by M - - M = 4; - tx_pilot_bb_lpf = tx_pilot_bb_lpf(1:M:length(tx_pilot_bb_lpf)); - rx_fdm_bb_lpf = rx_fdm_bb_lpf(1:M:length(rx_fdm_bb_lpf)); - n /= M; - nc /= M; - - % correlate over a range of frequency offsets and delays - - c_max = 0; - f_n = 1; - f_range = -75:2.5:75; - c_log=zeros(n, length(f_range)); - - for f=f_range - foff_rect = exp(j*2*pi*(f*M)*(1:nc)/Fs); - for s=1:n - - c = abs(tx_pilot_bb_lpf(1:nc)' * (rx_fdm_bb_lpf(s:s+nc-1) .* foff_rect')); - c_log(s,f_n) = c; - if c > c_max - c_max = c; - f_max = f; - s_max = s; - end - end - f_n++; - %printf("f: %f c_max: %f f_max: %f s_max: %d\n", f, c_max, f_max, s_max); - end - - figure(1); - y = f_range; - x = max(s_max-25,1):min(s_max+25, n); - mesh(y,x, c_log(x,:)); - grid - - s_max *= M; - s_max -= floor(s_max/6400)*6400; - printf("f_max: %f s_max: %d\n", f_max, s_max); - - % decimated position at sample rate. need to relate this to symbol - % rate position. - -endfunction - -function snr = EsNo_to_SNR(EsNo) - snr = interp1([20 11.8 9.0 6.7 4.9 3.3 -10], [ 3 3 0 -3 -6 -9 -9], EsNo); -endfunction - -function gen_test_bits() - sim_in = standard_init(); - framesize = 160; - tx_bits = round(rand(1,framesize)); - test_bits_coh_file(tx_bits); -endfunction - -% Start simulations --------------------------------------- - -more off; -%test_curves(); -test_single(); -%rate_Fs_tx("tx_zero.raw"); -%rate_Fs_tx("tx.raw"); -%rate_Fs_rx("tx_-4dB.wav") -%rate_Fs_rx("tx.raw") -%test_freq_off_est("tx.raw",40,6400) -%gen_test_bits(); - diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m new file mode 100644 index 00000000..2c6e9527 --- /dev/null +++ b/codec2-dev/octave/tcohpsk.m @@ -0,0 +1,30 @@ +% tcohpsk.m +% David Rowe Oct 2014 +% +% Octave script that tests the C port of the coherent PSK modem. This +% script loads the output of unittest/tcohpsk.c and compares it to the +% output of the reference versions of the same functions written in +% Octave. + +rand('state',1); +randn('state',1); +graphics_toolkit ("gnuplot"); + +cohpsk; +autotest; + +n = 160; +frames = 35; +framesize = 160; + +load ../build_linux/unittest/tcohpsk_out.txt + +sim_in = standard_init(); +tx_bits = round(rand(1,framesize)); +tx_bits_log = []; +for i=1:frames + tx_bits_log = [tx_bits_log tx_bits]; +end + +stem_sig_and_error(1, 211, tx_bits_log_c(1:n), tx_bits_log(1:n) - tx_bits_log_c(1:n), 'tx bits', [1 n -1.5 1.5]) +check(tx_bits_log, tx_bits_log_c, 'tx_bits'); diff --git a/codec2-dev/octave/test_cohpsk.m b/codec2-dev/octave/test_cohpsk.m new file mode 100644 index 00000000..f03ed725 --- /dev/null +++ b/codec2-dev/octave/test_cohpsk.m @@ -0,0 +1,578 @@ +% test_cohpsk.m +% David Rowe Oct 2014 +% + +% Bunch of simulations to test coherent FDM QPSK with pilot based +% coherent detection, DSSS (diversity), and rate 1/2 LDPC + +% TODO +% [X] Nc carriers, 588 bit frames +% [X] FEC +% [X] pilot insertion and removal +% [ ] delay on parity carriers +% [X] pilot based phase est +% [ ] uncoded and coded frame sync +% [X] timing estimation, RN filtering, carrier FDM +% [ ] this file getting too big - refactor +% [ ] robust coarse timing +% [ ] Np knob on GUI +% [ ] experimental tuning on EnNo_, fading[], to optimise LDPC dec perf + +% reqd to make sure we can repeat tests exactly + +rand('state',1); +randn('state',1); +graphics_toolkit ("gnuplot"); + +cohpsk; + +function test_curves + + sim_in = standard_init(); + + sim_in.verbose = 1; + sim_in.plot_scatter = 1; + + sim_in.Esvec = 10; + sim_in.hf_sim = 1; + sim_in.Ntrials = 1000; + sim_in.Rs = 50; + sim_in.Nc = 9; + sim_in.Np = 4; + sim_in.Ns = 8; + sim_in.Nchip = 1; + sim_in.modulation = 'qpsk'; + sim_in.ldpc_code_rate = 0.5; + sim_in.ldpc_code = 0; + + sim_qpsk = ber_test(sim_in); + + sim_in.hf_sim = 0; + sim_in.plot_scatter = 0; + sim_in.Esvec = 10:20; + Ebvec = sim_in.Esvec - 10*log10(2); + BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10))); + + sim_in.Np = 0; + sim_in.Nchip = 1; + + sim_in.modulation = 'dqpsk'; + sim_dqpsk = ber_test(sim_in, 'dqpsk'); + sim_in.hf_sim = 1; + sim_in.hf_mag_only = 1; + sim_in.modulation = 'qpsk'; + sim_qpsk_hf_ideal = ber_test(sim_in, 'qpsk'); + sim_in.modulation = 'dqpsk'; + sim_in.hf_mag_only = 0; + sim_dqpsk_hf = ber_test(sim_in, 'dqpsk'); + sim_in.modulation = 'qpsk'; + sim_in.Ns = 4 + sim_in.Np = 2; + sim_qpsk_hf_pilot = ber_test(sim_in, 'qpsk'); + + figure(1); + clf; + semilogy(Ebvec, BER_theory,'r;QPSK theory;') + hold on; + semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;') + semilogy(sim_qpsk_hf_ideal.Ebvec, sim_qpsk_hf_ideal.BERvec,'b;QPSK HF ideal;') + semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'k;DQPSK HF;') + semilogy(sim_qpsk_hf_pilot.Ebvec, sim_qpsk_hf_pilot.BERvec,'r;QPSK Np=2 Ns=4 HF;') + hold off; + + xlabel('Eb/N0') + ylabel('BER') + grid("minor") + axis([min(Ebvec) max(Ebvec) 1E-3 1]) +endfunction + + +function test_single + + sim_in = standard_init(); + + sim_in.verbose = 1; + sim_in.plot_scatter = 1; + + sim_in.framesize = 160; + sim_in.Nc = 4; + sim_in.Rs = 50; + sim_in.Ns = 4; + sim_in.Np = 2; + sim_in.Nchip = 1; +% sim_in.ldpc_code_rate = 0.5; +% sim_in.ldpc_code = 1; + sim_in.ldpc_code_rate = 1; + sim_in.ldpc_code = 0; + + sim_in.Ntrials = 20; + sim_in.Esvec = 100; + sim_in.hf_sim = 0; + sim_in.hf_mag_only = 0; + sim_in.modulation = 'qpsk'; + + sim_in.do_write_pilot_file = 1; + + sim_qpsk_hf = ber_test(sim_in); + + fep=fopen("errors_450.bin","wb"); fwrite(fep, sim_qpsk_hf.ldpc_errors_log, "short"); fclose(fep); +endfunction + + +% Rate Fs test funcs ----------------------------------------------------------- + +function rate_Fs_tx(tx_filename) + sim_in = standard_init(); + + sim_in.verbose = 1; + sim_in.plot_scatter = 1; + + sim_in.framesize = 160; + sim_in.Nc = 4; + sim_in.Rs = 50; + sim_in.Ns = 4; + sim_in.Np = 2; + sim_in.Nchip = 1; + sim_in.ldpc_code_rate = 1; + sim_in.ldpc_code = 0; + + sim_in.Ntrials = 20; + sim_in.Esvec = 7; + sim_in.hf_sim = 1; + sim_in.hf_mag_only = 0; + sim_in.modulation = 'qpsk'; + + sim_in = symbol_rate_init(sim_in); + + prev_sym_tx = sim_in.prev_sym_tx; + code_param = sim_in.code_param; + tx_bits_buf = sim_in.tx_bits_buf; + framesize = sim_in.framesize; + rate = sim_in.ldpc_code_rate; + Ntrials = sim_in.Ntrials; + Rs = sim_in.Rs; + Fs = sim_in.Fs; + Nc = sim_in.Nc; + + M = Fs/Rs; + + EsNodB = sim_in.Esvec(1); + EsNo = 10^(EsNodB/10); + + rx_symb_log = []; av_tx_pwr = []; + + rn_coeff = gen_rn_coeffs(0.5, 1/Fs, Rs, 6, M); + tx_symb_buf = []; + + tx_bits = round(rand(1,framesize*rate)); + + for nn=1:Ntrials+2 + + [tx_symb tx_bits_out prev_sym_tx] = symbol_rate_tx(sim_in, tx_bits, code_param, prev_sym_tx); + tx_bits_buf(1:framesize) = tx_bits_buf(framesize+1:2*framesize); + tx_bits_buf(framesize+1:2*framesize) = tx_bits_out; + tx_symb_buf = [tx_symb_buf; tx_symb]; + end + + % zero pad and tx filter + + [m n] = size(tx_symb_buf); + zp = []; + for i=1:m + zrow = M*tx_symb_buf(i,:); + zp = [zp; zrow; zeros(M-1,Nc)]; + end + + for c=1:Nc + tx_filt(:,c) = filter(rn_coeff, 1, zp(:,c)); + end + + % upconvert to real IF and save to disk + + [m n] = size(tx_filt); + tx_fdm = zeros(1,m); + Fc = 1500; + for c=1:Nc + freq(c) = exp(j*2*pi*(Fc - c*Rs*1.5)/Fs); + end + phase_tx = ones(1,Nc); + %phase_tx = exp(j*2*pi*(0:Nc)/(Nc+1)); + + for c=1:Nc + for i=1:m + phase_tx(c) = phase_tx(c) * freq(c); + tx_fdm(i) = tx_fdm(i) + tx_filt(i,c)*phase_tx(c); + end + end + + tx_fdm = real(tx_fdm); + + %tx_fdm = compress(tx_fdm, 0.4); + %tx_fdm = sign(tx_fdm) .* (abs(tx_fdm) .^ 0.4); + %hpa_clip = max(abs(tx_fdm))*0.8 + %tx_fdm(find(abs(tx_fdm) > hpa_clip)) = hpa_clip; + + papr = max(tx_fdm.*conj(tx_fdm)) / mean(tx_fdm.*conj(tx_fdm)); + papr_dB = 10*log10(papr); + printf("PAPR: %4.2f dB\n", papr_dB); + + Ascale = 2000; + figure(1); + clf; + plot(Ascale*tx_fdm(1:8000)) + + ftx=fopen(tx_filename,"wb"); fwrite(ftx, Ascale*real(tx_fdm), "short"); fclose(ftx); + +endfunction + + +function rate_Fs_rx(rx_filename) + sim_in = standard_init(); + + sim_in.verbose = 1; + sim_in.plot_scatter = 1; + + sim_in.framesize = 160; + sim_in.Nc = 4; + sim_in.Rs = 50; + sim_in.Ns = 4; + sim_in.Np = 4; + sim_in.Nchip = 1; + sim_in.ldpc_code_rate = 1; + sim_in.ldpc_code = 0; + + sim_in.Ntrials = 10; + sim_in.Esvec = 40; + sim_in.hf_sim = 1; + sim_in.hf_mag_only = 0; + sim_in.modulation = 'qpsk'; + + sim_in = symbol_rate_init(sim_in); + + prev_sym_tx = sim_in.prev_sym_tx; + prev_sym_rx = sim_in.prev_sym_rx; + code_param = sim_in.code_param; + tx_bits_buf = sim_in.tx_bits_buf; + framesize = sim_in.framesize; + rate = sim_in.ldpc_code_rate; + Ntrials = sim_in.Ntrials; + Rs = sim_in.Rs; + Fs = sim_in.Fs; + Nc = sim_in.Nc; + Nsymbrowpilot = sim_in.Nsymbrowpilot; + pilot = sim_in.pilot; + Ns = sim_in.Ns; + Npilotsframe = sim_in.Npilotsframe; + + M = Fs/Rs; + + EsNodB = sim_in.Esvec(1); + EsNo = 10^(EsNodB/10); + + phi_log = []; amp_log = []; EsNo__log = []; + rx_symb_log = []; av_tx_pwr = []; + Terrs = Tbits = 0; + errors_log = []; Nerrs_log = []; + ldpc_Nerrs_log = []; ldpc_errors_log = []; + Ferrsldpc = Terrsldpc = Tbitsldpc = 0; + + rn_coeff = gen_rn_coeffs(0.5, 1/Fs, Rs, 6, M); + + tx_bits = round(rand(1,framesize*rate)); + + % read from disk + + Ascale = 2000; + frx=fopen(rx_filename,"rb"); rx_fdm = fread(frx, "short")/Ascale; fclose(frx); + + rx_fdm=sqrt(2)*rx_fdm; + + if 0 + % optionally add AWGN noise for testing calibration of Es//No measurement + + snr_3000Hz_dB = -9; + snr_4000Hz_lin = 10 ^ (snr_3000Hz_dB/10); + snr_4000Hz_lin *= (3000/4000); + variance = var(rx_fdm)/snr_4000Hz_lin; + rx_fdm += sqrt(variance)*randn(length(rx_fdm),1); + end + + figure(1) + plot(rx_fdm(800:1200)); + + % freq offset estimation + + printf("Freq offset and coarse timing est...\n"); + [f_max max_s_Fs] = test_freq_off_est(rx_filename, 1,5*6400); + f_max = 0; max_s_Fs = 4; + max_s = floor(max_s_Fs/M + 6); + + printf("Downconverting...\n"); + + [m n] = size(rx_fdm); + rx_symb = zeros(m,Nc); + Fc = 1500; + for c=1:Nc + freq(c) = exp(-j*2*pi*(f_max + Fc - c*Rs*1.5)/Fs); + end + phase_rx = ones(1,Nc); + rx_bb = zeros(m,Nc); + + for c=1:Nc + for i=1:m + phase_rx(c) = phase_rx(c) * freq(c); + rx_bb(i,c) = rx_fdm(i)*phase_rx(c); + end + end + + printf("Filtering...\n"); + for c=1:Nc + rx_filt(:,c) = filter(rn_coeff, 1, rx_bb(:,c)); + end + + %subplot(211); + %plot(real(rx_filt(1:10*M,9))); + %subplot(212); + %plot(imag(rx_filt(1:10*M,9))); + + % Fine timing estimation and decimation to symbol rate Rs. Break rx + % signal into ft sample blocks. If clock offset is 1000ppm, + % that's one more/less sample over Ft samples at Fs=8000 Hz. + + printf("Fine timing estimation....\n"); + ft = M*10; + [nsam m] = size(rx_filt); + rx_symb_buf = []; rx_timing_log = []; + + for st=1:ft:floor(nsam/ft - 1)*ft + % fine timing and decimation + + env = zeros(ft,1); + for c=1:Nc + env = env + abs(rx_filt(st:st+ft-1,c)); + end + + % The envelope has a frequency component at the symbol rate. The + % phase of this frequency component indicates the timing. So work out + % single DFT at frequency 2*pi/M + + x = exp(-j*2*pi*(0:ft-1)/M) * env; + + norm_rx_timing = angle(x)/(2*pi); + %norm_rx_timing = -0.4; + if norm_rx_timing < 0 + rx_timing = -floor(norm_rx_timing*M+0.5) + M; + else + rx_timing = -floor(norm_rx_timing*M+0.5) + 2*M; + end + rx_timing_log = [rx_timing_log norm_rx_timing]; + + % printf("%d %d\n", st+rx_timing, st+rx_timing+ft-1); + rx_symb_buf = [rx_symb_buf; rx_filt(st+rx_timing:M:st+rx_timing+ft-1,:)]; + end + + figure(2) + clf; + plot(rx_timing_log) + axis([1 length(rx_timing_log) -0.5 0.5 ]) + title('fine timing') + + printf("Symbol rate demodulation....\n"); + phase_off = 1; + Ntrials = floor((nsam/M)/Nsymbrowpilot) - 2; + %max_s = 6; + + for nn=1:Ntrials + + s_ch = rx_symb_buf((nn-1)*Nsymbrowpilot+max_s:nn*Nsymbrowpilot+max_s-1,:); + [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx sim_in] = symbol_rate_rx(sim_in, s_ch, prev_sym_rx); + + rx_symb_log = [rx_symb_log rx_symb_linear]; + phi_log = [phi_log; phi_]; + amp_log = [amp_log; amp_]; + + if nn > 1 + EsNo__log = [EsNo__log EsNo_]; + + % Measure BER + + error_positions = xor(rx_bits(1:framesize*rate), tx_bits(1:framesize*rate)); + Nerrs = sum(error_positions); + Terrs += Nerrs; + Tbits += framesize*rate; + errors_log = [errors_log error_positions]; + Nerrs_log = [Nerrs_log Nerrs]; + + if sim_in.ldpc_code + % LDPC decode + + detected_data = ldpc_dec(code_param, sim_in.max_iterations, sim_in.demod_type, sim_in.decoder_type, ... + rx_symb_linear, min(100,EsNo_), amp_linear); + error_positions = xor(detected_data(1:framesize*rate), tx_bits(1:framesize*rate) ); + Nerrs = sum(error_positions); + ldpc_Nerrs_log = [ldpc_Nerrs_log Nerrs]; + ldpc_errors_log = [ldpc_errors_log error_positions]; + if Nerrs + Ferrsldpc++; + end + Terrsldpc += Nerrs; + Tbitsldpc += framesize*rate; + end + end + end + + EsNo_av = mean(10*log10(EsNo__log)); + printf("EsNo est (dB): %3.1f SNR est: %3.2f Terrs: %d BER %4.2f QPSK BER theory %4.2f av_tx_pwr: %3.2f", + EsNo_av, mean(EsNo_to_SNR(10*log10(EsNo__log))), + Terrs, + Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), av_tx_pwr); + if sim_in.ldpc_code + printf("\n LDPC: Terrs: %d BER: %4.2f Ferrs: %d FER: %4.2f\n", + Terrsldpc, Terrsldpc/Tbitsldpc, Ferrsldpc, Ferrsldpc/Ntrials); + end + + figure(3); + clf; + scat = rx_symb_log .* exp(j*pi/4); + scat = scat((framesize):length(scat)); + plot(real(scat), imag(scat),'+'); + title('Scatter plot'); + + figure(4) + clf + subplot(211) + stem(Nerrs_log) + axis([0 Ntrials+1 0 max(Nerrs_log)+1]) + title('Uncoded Errors/Frame'); + if sim_in.ldpc_code + subplot(212) + stem(ldpc_Nerrs_log) + axis([0 Ntrials+1 0 max(ldpc_Nerrs_log)+1]) + title('Coded Errors/Frame'); + end + + figure(5) + clf + [m1 n1] = size(phi_log); + phi_x = []; + phi_x_counter = 1; + p = Ns; + for r=1:m1 + if p == Ns + phi_x_counter++; + p = 0; + end + p++; + phi_x = [phi_x phi_x_counter++]; + end + + subplot(211) + plot(phi_x, phi_log(:,1),'r+;Estimated HF channel phase;') + ylabel('Phase (rads)'); + + subplot(212) + plot(phi_x, amp_log(:,1),'r+;Estimated HF channel amp;') + ylabel('Amplitude'); + xlabel('Time (symbols)'); + + figure(6); + clf + plot(EsNo_to_SNR(10*log10(EsNo__log))); + title('SNR est (dB)') + + fep=fopen("errors_450.bin","wb"); fwrite(fep, ldpc_errors_log, "short"); fclose(fep); + +endfunction + +% ideas: cld estimate timing with freq offset and decimate to save cpu load +% fft to do cross correlation + +function [f_max s_max] = test_freq_off_est(rx_filename, offset, n) + fpilot = fopen("tx_zero.raw","rb"); tx_pilot = fread(fpilot, "short"); fclose(fpilot); + frx=fopen(rx_filename,"rb"); rx_fdm = fread(frx, "short"); fclose(frx); + + Fs = 8000; + nc = 1800; % portion we wish to correlate over (first 2 rows on pilots) + + % downconvert to complex baseband to remove images + + f = 1000; + foff_rect = exp(j*2*pi*f*(1:2*n)/Fs); + tx_pilot_bb = tx_pilot(1:n) .* foff_rect(1:n)'; + rx_fdm_bb = rx_fdm(offset:offset+2*n-1) .* foff_rect'; + + % remove -2000 Hz image + + b = fir1(50, 1000/Fs); + tx_pilot_bb_lpf = filter(b,1,tx_pilot_bb); + rx_fdm_bb_lpf = filter(b,1,rx_fdm_bb); + + % decimate by M + + M = 4; + tx_pilot_bb_lpf = tx_pilot_bb_lpf(1:M:length(tx_pilot_bb_lpf)); + rx_fdm_bb_lpf = rx_fdm_bb_lpf(1:M:length(rx_fdm_bb_lpf)); + n /= M; + nc /= M; + + % correlate over a range of frequency offsets and delays + + c_max = 0; + f_n = 1; + f_range = -75:2.5:75; + c_log=zeros(n, length(f_range)); + + for f=f_range + foff_rect = exp(j*2*pi*(f*M)*(1:nc)/Fs); + for s=1:n + + c = abs(tx_pilot_bb_lpf(1:nc)' * (rx_fdm_bb_lpf(s:s+nc-1) .* foff_rect')); + c_log(s,f_n) = c; + if c > c_max + c_max = c; + f_max = f; + s_max = s; + end + end + f_n++; + %printf("f: %f c_max: %f f_max: %f s_max: %d\n", f, c_max, f_max, s_max); + end + + figure(1); + y = f_range; + x = max(s_max-25,1):min(s_max+25, n); + mesh(y,x, c_log(x,:)); + grid + + s_max *= M; + s_max -= floor(s_max/6400)*6400; + printf("f_max: %f s_max: %d\n", f_max, s_max); + + % decimated position at sample rate. need to relate this to symbol + % rate position. + +endfunction + +function snr = EsNo_to_SNR(EsNo) + snr = interp1([20 11.8 9.0 6.7 4.9 3.3 -10], [ 3 3 0 -3 -6 -9 -9], EsNo); +endfunction + +function gen_test_bits() + sim_in = standard_init(); + framesize = 160; + tx_bits = round(rand(1,framesize)); + test_bits_coh_file(tx_bits); +endfunction + +% Start simulations --------------------------------------- + +more off; +%test_curves(); +test_single(); +%rate_Fs_tx("tx_zero.raw"); +%rate_Fs_tx("tx.raw"); +%rate_Fs_rx("tx_-4dB.wav") +%rate_Fs_rx("tx.raw") +%test_freq_off_est("tx.raw",40,6400) +%gen_test_bits(); + diff --git a/codec2-dev/octave/tfdmdv.m b/codec2-dev/octave/tfdmdv.m index 1ea15a46..ab20e21e 100644 --- a/codec2-dev/octave/tfdmdv.m +++ b/codec2-dev/octave/tfdmdv.m @@ -12,7 +12,8 @@ more off NumCarriers = 14; fdmdv; % load modem code - +autotest; + % Generate reference vectors using Octave implementation of FDMDV modem global passes; @@ -171,44 +172,6 @@ end load ../unittest/tfdmdv_out.txt -% Helper functions to plot output of C verson and difference between Octave and C versions - -function stem_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec) - global no_plot_list; - - if find(no_plot_list == plotnum) - return; - end - figure(plotnum) - subplot(subplotnum) - stem(sig,'g;Octave version;'); - hold on; - stem(error,'r;Octave - C version (hopefully 0);'); - hold off; - if nargin == 6 - axis(axisvec); - end - title(titlestr); -endfunction - -function plot_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec) - global no_plot_list; - - if find(no_plot_list == plotnum) - return; - end - - figure(plotnum) - subplot(subplotnum) - plot(sig,'g;Octave version;'); - hold on; - plot(error,'r;Octave - C version (hopefully 0);'); - hold off; - if nargin == 6 - axis(axisvec); - end - title(titlestr); -endfunction % --------------------------------------------------------------------------------------- % Plot output and test each C function @@ -281,30 +244,6 @@ stem_sig_and_error(13, 212, imag(rx_symbols_log(:,f)), imag(rx_symbols_log(:,f) stem_sig_and_error(17, 211, real(phase_difference_log(:,f)), real(phase_difference_log(:,f) - phase_difference_log_c(:,f)), 'phase difference real' ) stem_sig_and_error(17, 212, imag(phase_difference_log(:,f)), imag(phase_difference_log(:,f) - phase_difference_log_c(:,f)), 'phase difference imag' ) -% --------------------------------------------------------------------------------------- -% AUTOMATED CHECKS ------------------------------------------ -% --------------------------------------------------------------------------------------- - -function check(a, b, test_name) - global passes; - global fails; - - [m n] = size(a); - printf("%s", test_name); - for i=1:(25-length(test_name)) - printf("."); - end - printf(": "); - - e = sum(abs(a - b))/n; - if e < 1E-3 - printf("OK\n"); - passes++; - else - printf("FAIL\n"); - fails++; - end -endfunction check(tx_bits_log, tx_bits_log_c, 'tx_bits'); check(tx_symbols_log, tx_symbols_log_c, 'tx_symbols'); diff --git a/codec2-dev/src/cohpsk.c b/codec2-dev/src/cohpsk.c index 46530955..81e8e9e2 100644 --- a/codec2-dev/src/cohpsk.c +++ b/codec2-dev/src/cohpsk.c @@ -77,7 +77,7 @@ void bits_to_qpsk_symbols(COMP tx_symb[][PILOTS_NC], int tx_bits[], int framesiz short bits; assert(COHPSK_NC == PILOTS_NC); - assert((NSYMROW*PILOTS_NC)/2 == framesize); + assert((NSYMROW*PILOTS_NC)*2 == framesize); /* Organise QPSK symbols into a NSYMBROWS rows by PILOTS_NC cols matrix, @@ -109,9 +109,11 @@ void bits_to_qpsk_symbols(COMP tx_symb[][PILOTS_NC], int tx_bits[], int framesiz /* NS rows of data symbols */ - for(i=0; i