+ e.g. hang on thru the fades?
[ ] PAPR idea
+ automatically tweak phases to reduce PAPR, e.g. slow variations in freq...
-
+[ ] why is pilot noise_est twice as big as other carriers
global Nfilter = Nsym*M;
global Nfiltertiming = M+Nfilter+M;
alpha = 0.5;
+global snr_coeff = 0.9 % SNR est averaging filter coeff
% root raised cosine (Root Nyquist) filter
% convert symbols back to an array of bits
-function [rx_bits sync_bit f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
+function [rx_bits sync_bit f_err phase_difference] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
global Nc;
global Nb;
global Nb;
% extra 45 degree clockwise lets us use real and imag axis as
% decision boundaries
+ phase_difference = zeros(Nc+1,1);
phase_difference(1:Nc) = rx_symbols(1:Nc) .* conj(prev_rx_symbols(1:Nc)) * exp(j*pi/4);
% map (Nc,1) DQPSK symbols back into an (1,Nc*Nb) array of bits
% Extract DBPSK encoded Sync bit
- phase_difference(Nc+1) = rx_symbols(Nc+1) .* conj(prev_rx_symbols(Nc+1));
+ phase_difference(Nc+1,1) = rx_symbols(Nc+1) .* conj(prev_rx_symbols(Nc+1));
if (real(phase_difference(Nc+1)) < 0)
sync_bit = 1;
f_err = imag(phase_difference(Nc+1));
f_err = -imag(phase_difference(Nc+1));
end
+ % pilot carrier gets an extra pi/4 rotation to make it consistent with
+ % other carriers, as we need it for snr_update and scatter diagram
+
+ phase_difference(Nc+1) *= exp(j*pi/4);
else
% map (Nc,1) QPSK symbols back into an (1,Nc*Nb) array of bits
endfunction
+% given phase differences update estimates of signal and noise levels
+
+function [sig_est noise_est] = snr_update(sig_est, noise_est, phase_difference)
+ global snr_coeff;
+ global Nc;
+
+ % mag of each symbol is distance from origin, this gives us a
+ % vector of mags, one for each carrier.
+
+ s = abs(phase_difference);
+
+ % signal mag estimate for each carrier is a smoothed version
+ % of instantaneous magntitude, this gives us a vector of smoothed
+ % mag estimates, one for each carrier.
+
+ sig_est = snr_coeff*sig_est + (1 - snr_coeff)*s;
+
+ % noise mag estimate is distance of current symbol from average
+ % location of that symbol. We reflect all symbols into the first
+ % quadrant for convenience.
+
+ refl_symbols = abs(real(phase_difference)) + j*abs(imag(phase_difference));
+ n = abs(exp(j*pi/4)*sig_est - refl_symbols);
+
+ % noise mag estimate for each carrier is a smoothed version of
+ % instantaneous noise mag, this gives us a vector of smoothed
+ % noise power estimates, one for each carrier.
+
+ noise_est = snr_coeff*noise_est + (1 - snr_coeff)*n;
+endfunction
+
+
+% calculate current SNR estimate (3000Hz noise BW)
+
+function snr_dB = calc_snr(sig_est, noise_est)
+ global Rs;
+
+ % find total signal power by summing power in all carriers
+
+ S = sum(sig_est .^2);
+ SdB = 10*log10(S);
+
+ % Average noise mag across all carriers and square to get an average
+ % noise power. This is an estimate of the noise power in Rs = 50Hz of
+ % BW (note for raised root cosine filters Rs is the noise BW of the
+ % filter)
+
+ N50 = mean(noise_est).^2;
+ N50dB = 10*log10(N50);
+
+ % Now multiply by (3000 Hz)/(50 Hz) to find the total noise power in
+ % 3000 Hz
+
+ N3000dB = N50dB + 10*log10(3000/Rs);
+
+ snr_dB = SdB - N3000dB;
+
+endfunction
+
% returns nbits from a repeating sequence of random data
function bits = get_test_bits(nbits)
global rx_test_bits_mem;
rx_test_bits_mem = zeros(1,Ntest_bits);
-
test_frame_sync_log = [];
test_frame_sync_state = 0;
+ % SNR states
+
+ sig_est = zeros(Nc+1,1);
+ noise_est = zeros(Nc+1,1);
+
+ % logs of various states for plotting
+
rx_symbols_log = [];
rx_timing_log = [];
foff_log = [];
rx_fdm_log = [];
+ snr_est_log = [];
% misc states
else
rx_symbols_log = [rx_symbols_log rx_symbols];
endif
- [rx_bits sync f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+ [rx_bits sync f_err pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+ [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
+ snr_est = calc_snr(sig_est, noise_est);
+ snr_est_log = [snr_est_log snr_est];
foff -= 0.5*f_err;
prev_rx_symbols = rx_symbols;
sync_log = [sync_log sync];
axis([0 secs 0 1.5]);
title('Test Frame Sync')
-
-
+ figure(5)
+ clf;
+ plot(xt, snr_est_log);
+ title('SNR Estimates')
+
endfunction
axis([0 secs 0 1.5]);
title('Test Frame Sync')
+ figure(5)
+ clf;
+ plot(xt, snr_est_log_c(1:frames));
+ title('SNR Estimates')
endfunction
% Simulation Parameters --------------------------------------
-frames = 100;
+frames = 25;
EbNo_dB = 7.3;
Foff_hz = 0;
modulation = 'dqpsk';
test_frame_sync_log = [];
test_frame_sync_state = 0;
+% SNR estimation states
+
+sig_est = zeros(Nc+1,1);
+noise_est = zeros(Nc+1,1);
+
% fixed delay simuation
Ndelay = M+20;
track = 0;
track_log = [];
+
% ---------------------------------------------------------------------
% Main loop
% ---------------------------------------------------------------------
%rx_phase_log = [rx_phase_log rx_phase];
%rx_symbols = rx_symbols*exp(j*rx_phase);
+ [rx_bits sync foff_fine pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
if strcmp(modulation,'dqpsk')
- rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)];
+ %rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)];
+ rx_symbols_log = [rx_symbols_log pd];
else
rx_symbols_log = [rx_symbols_log rx_symbols];
endif
- [rx_bits sync foff_fine] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
foff -= 0.5*ferr;
prev_rx_symbols = rx_symbols;
sync_log = [sync_log sync];
[track fest_state] = freq_state(sync, fest_state);
track_log = [track_log track];
+ % Update SNR est
+
+ [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
+
% count bit errors if we find a test frame
% Allow 15 frames for filter memories to fill and time est to settle
ber = total_bit_errors / total_bits;
-% Peak to Average Power Ratio from http://www.dsplog.com
+% Peak to Average Power Ratio calcs from http://www.dsplog.com
papr = max(tx_fdm_log.*conj(tx_fdm_log)) / mean(tx_fdm_log.*conj(tx_fdm_log));
papr_dB = 10*log10(papr);
-printf("Eb/No (meas): %2.2f (%2.2f) dB\nbits........: %d\nerrors......: %d\nBER.........: %1.4f\nPAPR........: %1.2f dB\nSNR.........: %2.1f dB\n",
- EbNo_dB, 10*log10(0.25*tx_pwr*Fs/(Rs*Nc*noise_pwr)),
- total_bits, total_bit_errors, ber, papr_dB, SNR );
+% Note Eb/No set point is for Nc data carriers only, exclduing pilot.
+% This is convenient for testing BER versus Eb/No. Measured Eb/No
+% includes power of pilot. Similar for SNR, first number is SNR excluding
+% pilot pwr for Eb/No set point, 2nd value is measured SNR which will be a little
+% higher as pilot power is included.
+
+printf("Eb/No (meas): %2.2f (%2.2f) dB\n", EbNo_dB, 10*log10(0.25*tx_pwr*Fs/(Rs*Nc*noise_pwr)));
+printf("bits........: %d\n", total_bits);
+printf("errors......: %d\n", total_bit_errors);
+printf("BER.........: %1.4f\n", ber);
+printf("PAPR........: %1.2f dB\n", papr_dB);
+printf("SNR...(meas): %2.2f (%2.2f) dB\n", SNR, calc_snr(sig_est, noise_est));
% ---------------------------------------------------------------------
% Plots
function phase2(samname, png)
N = 16000;
- f=45;
+ f=43;
model = load("../src/hts1a_phase_model.txt");
phase = load("../src/hts1a_phase_phase.txt");
Wo = model(f,1);
L = model(f,2);
A = model(f,3:(L+2));
phi = phase(f,1:L);
+ phi = zeros(1,L);
+ phi(3) = -pi/2;
+ phi(4) = -pi/4;
+ phi(5) = pi/2;
s = zeros(1,N);
- for m=1:L/2
+ for m=3:5
s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
s = s + s_m;
endfor
clf;
plot(s(1:250));
+ figure(2);
+ clf;
+ subplot(211)
+ plot((1:L)*Wo*4000/pi, 20*log10(A),'+');
+ subplot(212)
+ plot((1:L)*Wo*4000/pi, phi,'+');
+
fs=fopen(samname,"wb");
fwrite(fs,s,"short");
fclose(fs);
clf;
% s = [ Sn(2*(f-2)-1,:) Sn(2*(f-2),:) ];
s = [ Sn(2*f-1,:) Sn(2*f,:) ];
+ size(s)
plot(s);
axis([1 length(s) -20000 20000]);
channel = [];
channel_count = 0;
next_nin = M;
+sig_est = zeros(Nc+1,1);
+noise_est = zeros(Nc+1,1);
% Octave outputs we want to collect for comparison to C version
sync_bit_log = [];
coarse_fine_log = [];
nin_log = [];
+sig_est_log = [];
+noise_est_log = [];
for f=1:frames
end
nin_log = [nin_log nin];
- [rx_bits sync_bit foff_fine] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
+ [rx_bits sync_bit foff_fine pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
+
+ [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
+ sig_est_log = [sig_est_log sig_est];
+ noise_est_log = [noise_est_log noise_est];
+
prev_rx_symbols = rx_symbols;
rx_bits_log = [rx_bits_log rx_bits];
foff_fine_log = [foff_fine_log foff_fine];
stem_sig_and_error(15, 211, rx_bits_log_c(st:en), rx_bits_log(st:en) - rx_bits_log_c(st:en), 'RX bits', [1 en-st -1.5 1.5])
stem_sig_and_error(15, 212, nin_log_c, nin_log - nin_log_c, 'nin')
+c = 1;
+plot_sig_and_error(16, 211, sig_est_log(c,:), sig_est_log(c,:) - sig_est_log_c(c,:), 'sig est for SNR' )
+plot_sig_and_error(16, 212, noise_est_log(c,:), noise_est_log(c,:) - noise_est_log_c(c,:), 'noise est for SNR' )
+
% ---------------------------------------------------------------------------------------
% AUTOMATED CHECKS ------------------------------------------
% ---------------------------------------------------------------------------------------
check(sync_bit_log, sync_bit_log_c, 'sync bit');
check(coarse_fine_log, coarse_fine_log_c, 'coarse_fine');
check(nin_log, nin_log_c, 'nin');
+check(sig_est_log, sig_est_log_c, 'sig_est');
+check(noise_est_log, noise_est_log_c, 'noise_est');
printf("\npasses: %d fails: %d\n", passes, fails);
int vector_quant_Wo_e = 0;
int dump_pitch_e = 0;
FILE *fjvm = NULL;
+ #ifdef DUMP
+ int dump;
+ #endif
char* opt_string = "ho:";
struct option long_options[] = {
return res;
}
+static float cabsolute(COMP a)
+{
+ return sqrt(pow(a.real, 2.0) + pow(a.imag, 2.0));
+}
+
/*---------------------------------------------------------------------------*\
FUNCTION....: fdmdv_create
f->fest_state = 0;
f->coarse_fine = COARSE;
+
+ for(c=0; c<NC+1; c++) {
+ f->sig_est[c] = 0.0;
+ f->noise_est[c] = 0.0;
+ }
return f;
}
for(i=0; i<NT*P; i++) {
env[i] = 0.0;
for(c=0; c<NC+1; c++)
- env[i] += sqrt(pow(rx_filter_mem_timing[c][i].real,2.0) + pow(rx_filter_mem_timing[c][i].imag,2.0));
+ env[i] += cabsolute(rx_filter_mem_timing[c][i]);
}
/* The envelope has a frequency component at the symbol rate. The
*sync_bit = 0;
ferr = -phase_difference[NC].imag;
}
+
+ /* pilot carrier gets an extra pi/4 rotation to make it consistent
+ with other carriers, as we need it for snr_update and scatter
+ diagram */
+
+ phase_difference[NC] = cmult(phase_difference[NC], pi_on_4);
return ferr;
}
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: snr_update()
+ AUTHOR......: David Rowe
+ DATE CREATED: 17 May 2012
+
+ Given phase differences update estimates of signal and noise levels.
+
+\*---------------------------------------------------------------------------*/
+
+void snr_update(float sig_est[], float noise_est[], COMP phase_difference[])
+{
+ float s[NC+1];
+ COMP refl_symbols[NC+1];
+ float n[NC+1];
+ COMP pi_on_4;
+ int c;
+
+ pi_on_4.real = cos(PI/4.0);
+ pi_on_4.imag = sin(PI/4.0);
+
+ /* mag of each symbol is distance from origin, this gives us a
+ vector of mags, one for each carrier. */
+
+ for(c=0; c<NC+1; c++)
+ s[c] = cabsolute(phase_difference[c]);
+
+ /* signal mag estimate for each carrier is a smoothed version of
+ instantaneous magntitude, this gives us a vector of smoothed
+ mag estimates, one for each carrier. */
+
+ for(c=0; c<NC+1; c++)
+ sig_est[c] = SNR_COEFF*sig_est[c] + (1.0 - SNR_COEFF)*s[c];
+
+ /* noise mag estimate is distance of current symbol from average
+ location of that symbol. We reflect all symbols into the first
+ quadrant for convenience. */
+
+ for(c=0; c<NC+1; c++) {
+ refl_symbols[c].real = fabs(phase_difference[c].real);
+ refl_symbols[c].imag = fabs(phase_difference[c].imag);
+ n[c] = cabsolute(cadd(fcmult(sig_est[c], pi_on_4), cneg(refl_symbols[c])));
+ }
+
+ /* noise mag estimate for each carrier is a smoothed version of
+ instantaneous noise mag, this gives us a vector of smoothed
+ noise power estimates, one for each carrier. */
+
+ for(c=0; c<NC+1; c++)
+ noise_est[c] = SNR_COEFF*noise_est[c] + (1 - SNR_COEFF)*n[c];
+}
+
/*---------------------------------------------------------------------------*\
FUNCTION....: fdmdv_put_test_bits()
foff_fine = qpsk_to_bits(rx_bits, sync_bit, fdmdv->phase_difference, fdmdv->prev_rx_symbols, rx_symbols);
memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP)*(NC+1));
+ snr_update(fdmdv->sig_est, fdmdv->noise_est, fdmdv->phase_difference);
/* freq offset estimation state machine */
fdmdv->foff -= TRACK_COEFF*foff_fine;
}
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: calc_snr()
+ AUTHOR......: David Rowe
+ DATE CREATED: 17 May 2012
+
+ Calculate current SNR estimate (3000Hz noise BW)
+
+\*---------------------------------------------------------------------------*/
+
+float calc_snr(float sig_est[], float noise_est[])
+{
+ float S, SdB;
+ float mean, N50, N50dB, N3000dB;
+ float snr_dB;
+ int c;
+
+ S = 0.0;
+ for(c=0; c<NC+1; c++)
+ S += pow(sig_est[c], 2.0);
+ SdB = 10.0*log10(S);
+
+ /* Average noise mag across all carriers and square to get an
+ average noise power. This is an estimate of the noise power in
+ Rs = 50Hz of BW (note for raised root cosine filters Rs is the
+ noise BW of the filter) */
+
+ mean = 0.0;
+ for(c=0; c<NC+1; c++)
+ mean += noise_est[c];
+ mean /= (NC+1);
+ N50 = pow(mean, 2.0);
+ N50dB = 10.0*log10(N50);
+
+ /* Now multiply by (3000 Hz)/(50 Hz) to find the total noise power
+ in 3000 Hz */
+
+ N3000dB = N50dB + 10.0*log10(3000.0/RS);
+
+ snr_dB = SdB - N3000dB;
+
+ return snr_dB;
+}
+
/*---------------------------------------------------------------------------*\
FUNCTION....: fdmdv_get_demod_stats()
void fdmdv_get_demod_stats(struct FDMDV *fdmdv, struct FDMDV_STATS *fdmdv_stats)
{
int c;
- COMP pi_on_4;
- pi_on_4.real = cos(PI/4.0);
- pi_on_4.imag = sin(PI/4.0);
-
- fdmdv_stats->snr = 0.0; /* TODO - implement SNR estimation */
+ fdmdv_stats->snr_est = calc_snr(fdmdv->sig_est, fdmdv->noise_est);
fdmdv_stats->fest_coarse_fine = fdmdv->coarse_fine;
fdmdv_stats->foff = fdmdv->foff;
fdmdv_stats->rx_timing = fdmdv->rx_timing;
assert((NC+1) == FDMDV_NSYM);
- for(c=0; c<NC; c++) {
+ for(c=0; c<NC+1; c++) {
fdmdv_stats->rx_symbols[c] = fdmdv->phase_difference[c];
}
-
- /* place pilots somewhere convenient on scatter diagram */
-
- fdmdv_stats->rx_symbols[NC] = cmult(fdmdv->phase_difference[NC], pi_on_4);
-
}
/*---------------------------------------------------------------------------*\
struct FDMDV;
struct FDMDV_STATS {
- float snr; /* estimated SNR of rx signal in dB */
+ float snr_est; /* estimated SNR of rx signal in dB (3 kHz noise BW) */
COMP rx_symbols[FDMDV_NSYM]; /* latest received symbols, for scatter plot */
int fest_coarse_fine; /* freq est state, 0-coarse 1-fine */
float foff; /* estimated freq offset in Hz */
float foff_log[MAX_FRAMES];
int sync_bit_log[MAX_FRAMES];
int rx_bits_log[FDMDV_BITS_PER_FRAME*MAX_FRAMES];
+ float snr_est_log[MAX_FRAMES];
if (argc < 3) {
printf("usage: %s InputModemRawFile OutputBitFile [OctaveDumpFile]\n", argv[0]);
coarse_fine_log[f] = stats.fest_coarse_fine;
sync_bit_log[f] = sync_bit;
memcpy(&rx_bits_log[FDMDV_BITS_PER_FRAME*f], rx_bits, sizeof(int)*FDMDV_BITS_PER_FRAME);
-
+ snr_est_log[f] = stats.snr_est;
f++;
}
else
octave_save_int(foct, "coarse_fine_log_c", coarse_fine_log, 1, f);
octave_save_int(foct, "rx_bits_log_c", rx_bits_log, 1, FDMDV_BITS_PER_FRAME*f);
octave_save_int(foct, "sync_bit_log_c", sync_bit_log, 1, f);
+ octave_save_float(foct, "snr_est_log_c", snr_est_log, 1, f, MAX_FRAMES);
fclose(foct);
}
}
#define FINE 1
#define TRACK_COEFF 0.5
+#define SNR_COEFF 0.9 /* SNR est averaging filter coeff */
+
/*---------------------------------------------------------------------------*\
STRUCT for States
int fest_state;
int coarse_fine;
+
+ /* SNR estimation states */
+
+ float sig_est[NC+1];
+ float noise_est[NC+1];
};
/*---------------------------------------------------------------------------*\
COMP rx_baseband_mem_timing[NC+1][NFILTERTIMING],
int nin);
float qpsk_to_bits(int rx_bits[], int *sync_bit, COMP phase_difference[], COMP prev_rx_symbols[], COMP rx_symbols[]);
+void snr_update(float sig_est[], float noise_est[], COMP phase_difference[]);
int freq_state(int sync_bit, int *state);
+float calc_snr(float sig_est[], float noise_est[]);
#endif
float env_log[NT*P*FRAMES];
float rx_timing_log[FRAMES];
COMP rx_symbols_log[NC+1][FRAMES];
+ float sig_est_log[NC+1][FRAMES];
+ float noise_est_log[NC+1][FRAMES];
int rx_bits_log[FDMDV_BITS_PER_FRAME*FRAMES];
float foff_fine_log[FRAMES];
int sync_bit_log[FRAMES];
rx_filter(rx_filt, rx_baseband, fdmdv->rx_filter_memory, nin);
rx_timing = rx_est_timing(rx_symbols, rx_filt, rx_baseband, fdmdv->rx_filter_mem_timing, env, fdmdv->rx_baseband_mem_timing, nin);
foff_fine = qpsk_to_bits(rx_bits, &sync_bit, fdmdv->phase_difference, fdmdv->prev_rx_symbols, rx_symbols);
+ snr_update(fdmdv->sig_est, fdmdv->noise_est, fdmdv->phase_difference);
memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP)*(NC+1));
next_nin = M;
/* qpsk_to_bits() */
memcpy(&rx_bits_log[FDMDV_BITS_PER_FRAME*f], rx_bits, sizeof(int)*FDMDV_BITS_PER_FRAME);
+ for(c=0; c<NC+1; c++) {
+ sig_est_log[c][f] = fdmdv->sig_est[c];
+ noise_est_log[c][f] = fdmdv->noise_est[c];
+ }
foff_fine_log[f] = foff_fine;
sync_bit_log[f] = sync_bit;
octave_save_float(fout, "env_log_c", env_log, 1, NT*P*FRAMES, NT*P*FRAMES);
octave_save_float(fout, "rx_timing_log_c", rx_timing_log, 1, FRAMES, FRAMES);
octave_save_complex(fout, "rx_symbols_log_c", (COMP*)rx_symbols_log, (NC+1), FRAMES, FRAMES);
+ octave_save_float(fout, "sig_est_log_c", (float*)sig_est_log, (NC+1), FRAMES, FRAMES);
+ octave_save_float(fout, "noise_est_log_c", (float*)noise_est_log, (NC+1), FRAMES, FRAMES);
octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1, FDMDV_BITS_PER_FRAME*FRAMES);
octave_save_float(fout, "foff_fine_log_c", foff_fine_log, 1, FRAMES, FRAMES);
octave_save_int(fout, "sync_bit_log_c", sync_bit_log, 1, FRAMES);
projects / freetel-svn-tracking.git / commitdiff