--- /dev/null
+
+
+Modeling sample clock errors using sox:
+
+sox -r 8000 -s -2 mod_dqpsk.raw -s -2 mod_dqpsk_8008hz.raw rate -h 8008
% Version 2
%
+% reqd to mak sure we get same random bits at mod and demod
+
+rand('state',1);
+randn('state',1);
+
% Constants
global Fs = 8000; % sample rate in Hz
msb = tx_bits_matrix(c,1); lsb = tx_bits_matrix(c,2);
if ((msb == 0) && (lsb == 0))
- tx_symbols(c) = prev_tx_symbols(c);
+ tx_symbols(c) = prev_tx_symbols(c);
endif
if ((msb == 0) && (lsb == 1))
tx_symbols(c) = j*prev_tx_symbols(c);
end
else
% QPSK mapping
- tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2);
+ tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2);
endif
endfunction
tx_fdm(i) = tx_fdm(i) + tx_filt(c,i)*phase_tx(c);
end
end
-
+
% Nc/2 tones above centre freq
for c=Nc/2+1:Nc
end
end
- tx_fdm = real(tx_fdm);
+ % add centre pilot tone at twice amplitude of other tones
+
+ c = Nc+1;
+ for i=1:M
+ phase_tx(c) = phase_tx(c) * freq(c);
+ tx_fdm(i) = tx_fdm(i) + 1/sqrt(2)*phase_tx(c);
+ end
+
+ % Scale such that total Carrier power C of real(tx_fdm) = Nc
+ % we return the complex (single sided) signal to make frequency
+ % shifting for the purpose of testing easier
+
+ tx_fdm = sqrt(2)*tx_fdm;
endfunction
endfunction
-% Estimate optimum timing offset, and symbol receive symbols
+% Estimate frequency offset of FDM signal by peak picking pilot
+% tone at Fcentre
+function foff = rx_est_freq_offset(rx_fdm)
+ global M;
+ global Nc;
+ global Fs;
+ global freq;
+ global phase_rx;
+ global Npilotbaseband;
+ global Npilotlpf;
+ global Npilotcoeff;
+ global pilot_baseband;
+ global pilot_lpf;
+ global pilot_coeff;
+
+ % down convert latest M samples of pilot
+
+ pilot_baseband(1:Npilotbaseband-M) = pilot_baseband(M+1:Npilotbaseband);
+ c = Nc+1;
+ for i=1:M
+ phase_rx(c) = phase_rx(c) * freq(c);
+ pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i) * phase_rx(c)';
+ end
+
+ % LPF cutoff 200Hz
+
+ pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf);
+ j = 1;
+ for i = Npilotlpf-M+1:Npilotlpf
+ pilot_lpf(i) = pilot_baseband(j:j+Npilotcoeff) * pilot_coeff';
+ j++;
+ end
+
+ % decimate to improve DFT resolution, window and DFT
+
+ Mpilot = Fs/(2*200); % calc decimation rate given new sample rate is twice LPF freq
+ Mpilotfft = 256;
+ s = pilot_lpf(1:Mpilot:Npilotlpf) .* hanning(Npilotlpf/Mpilot)';
+ S = abs(fft(s, Mpilotfft));
+
+ % peak pick and convert to Hz
+
+ [x ix] = max(S);
+ r = 2*200/Mpilotfft; % maps FFT bin to frequency in Hz
+
+ if ix > Mpilotfft/2
+ foff = (ix - Mpilotfft - 1)*r;
+ else
+ foff = (ix - 1)*r;
+ endif
+
+endfunction
+
+
+% Estimate optimum timing offset, and symbol receive symbols
function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
global M;
% map phase to estimated optimum timing instant at rate M
% the M/2 + 1 part was adjusted by experment, I know not why....
- rx_timing = angle(x)*M/pi + M/2 + 1;
+ rx_timing = angle(x)*M/(2*pi) + M/4;
if (rx_timing > M)
rx_timing -= M;
end
% Phase estimation function - probably won't work over a HF channel
% Tries to operate over a sinle symbol but uses phase information from
% all Nc carriers which should increase the SNR of phase estimate.
-% Maybe phase is coherent over a coupel of symbols in HF channel,not
+% Maybe phase is coherent over a couple of symbols in HF channel,not
% sure but it's worth 3dB so worth experimenting or using coherent as
% an option.
% if less than a thresh we are aligned and in sync
- ber = bit_errors/Ntest_bits;
+ ber = bit_errors/Ntest_bits;
sync = 0;
- if (ber < 0.2)
+ if (ber < 0.1)
sync = 1;
endif
endfunction
% phasors used for up and down converters
-global freq = zeros(Nc,1);;
+global 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);
freq(c) = exp(j*2*pi*carrier_freq/Fs);
end
-global phase_tx = ones(Nc,1);
-global phase_rx = ones(Nc,1);
+freq(Nc+1) = exp(j*2*pi*Fcentre/Fs);
+
+global phase_tx = ones(Nc+1,1);
+global phase_rx = ones(Nc+1,1);
+
+% Freq offset estimator states
+
+global Npilotcoeff = 30; % number of pilot LPF coeffs
+global pilot_coeff = fir1(Npilotcoeff, 200/(Fs/2))'; % 200Hz LPF
+global Npilotbaseband = Npilotcoeff + 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
-% Timing estimateor states
+% Timing estimator states
global rx_filter_mem_timing = zeros(Nc, Nt*P);
global rx_baseband_mem_timing = zeros(Nc, Nfiltertiming);
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;
if (nargin == 1)
frames = floor(length(rx_fdm)/M);
else
- frames = nbits/(Nc*Nb);
+ frames = nbits/(Nc*Nb);
endif
+prev_rx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
+
total_bit_errors = 0;
total_bits = 0;
rx_timing_log = [];
rx_symbols_log = [];
rx_phase_log = [];
-prev_rx_symbols = ones(Nc,1)*exp(j*pi/4);
-modulation = 'dqpsk';
+sync_log = [];
% Main loop ----------------------------------------------------
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);
%rx_phase_log = [rx_phase_log rx_phase];
%rx_symbols = rx_symbols*exp(j*rx_phase);
+
if strcmp(modulation,'dqpsk')
rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)];
else
prev_rx_symbols = rx_symbols;
[sync bit_errors] = put_test_bits(rx_bits);
+ sync_log = [sync_log sync];
if sync == 1
total_bit_errors = total_bit_errors + bit_errors;
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(311)
+%subplot(211)
plot(rx_timing_log)
-subplot(312)
-Nfft=Fs;
-S=fft(rx_fdm,Nfft);
-SdB=20*log10(abs(S));
-plot(SdB(1:Fs/4))
-subplot(313)
-%plot(rx_phase_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
frames = floor(nbits/(Nc*Nb));
tx_fdm = [];
gain = 1000; % Scale up to 16 bit shorts
-prev_tx_symbols = ones(Nc,1)*exp(j*pi/4);
+prev_tx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
for i=1:frames
tx_bits = get_test_bits(Nc*Nb);
- tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits,'dqpsk');
+ 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 fdm_upconvert(tx_baseband)];
+tx_fdm = [tx_fdm real(fdm_upconvert(tx_baseband))];
end
tx_fdm *= gain;
clear all;
fdmdv; % load modem code
-rand('state',1);
-randn('state',1);
% Simulation Parameters --------------------------------------
-frames = 50;
-EbNo_dB = 40;
-Foff_hz = 1;
+frames = 100;
+EbNo_dB = 7.3;
+Foff_hz = 124;
modulation = 'dqpsk';
% ------------------------------------------------------------
total_bit_errors = 0;
total_bits = 0;
rx_fdm_log = [];
+rx_baseband_log = [];
rx_bits_offset = zeros(Nc*Nb*2);
-prev_tx_symbols = ones(Nc,1)*exp(j*pi/4);
-prev_rx_symbols = ones(Nc,1)*exp(j*pi/4);
+prev_tx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
+prev_rx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
+foff_log = [];
+
+Ndelay = M+20;
+rx_fdm_delay = zeros(Ndelay,1);
% Eb/No calculations. We need to work out Eb/No for each FDM carrier.
% Total power is sum of power in all FDM carriers
-C = 1; % power of each FDM carrier (energy/sample)
-N = 1; % total noise power (energy/sample) of noise source before scaling
- % by Ngain
+C = 1; % power of each FDM carrier (energy/sample). Total Carrier power should = Nc*C = Nc
+N = 1; % total noise power (energy/sample) of noise source across entire bandwidth
% Eb = Carrier power * symbol time / (bits/symbol)
-% = C *(Rs/Fs) / 2
-Eb_dB = 10*log10(C) + 10*log10(Rs) - 10*log10(Fs) - 10*log10(2);
+% = C *(1/Rs) / 2
+Eb_dB = 10*log10(C) - 10*log10(Rs) - 10*log10(2);
No_dBHz = Eb_dB - EbNo_dB;
-% Noise power = Noise spectral density * bandwidth;
-N_dB = No_dBHz + 10*log10(Fs);
+% Noise power = Noise spectral density * bandwidth
+% Noise power = Noise spectral density * Fs/2 for real signals
+N_dB = No_dBHz + 10*log10(Fs/2);
Ngain_dB = N_dB - 10*log10(N);
Ngain = 10^(Ngain_dB/20);
-% C/No = Carrier Power/noise spectral denity
-% = power per carrier*number of carriers / noise spectral denity
+% C/No = Carrier Power/noise spectral density
+% = power per carrier*number of carriers / noise spectral density
CNo_dB = 10*log10(C) + 10*log10(Nc) - No_dBHz;
-% SNR in equivalent 2400 Hz SSB channel
+% SNR in equivalent 3000 Hz SSB channel
-B = 2400;
+B = 3000;
SNR = CNo_dB - 10*log10(B);
phase_offset = 1;
freq_offset = exp(j*2*pi*Foff_hz/Fs);
+foff_phase = 1;
% Main loop ----------------------------------------------------
for i=1:frames
+
+ % -------------------
+ % Modulator
+ % -------------------
+
tx_bits = get_test_bits(Nc*Nb);
- %tx_bits = zeros(1,Nc*Nb);
- %tx_bits = ones(1,Nc*Nb);
- %tx_bits = [0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
- %tx_bits = [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0];
tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation);
prev_tx_symbols = tx_symbols;
tx_baseband = tx_filter(tx_symbols);
- tx_fdm = fdm_upconvert(tx_baseband);
- tx_pwr = 0.9*tx_pwr + 0.1*tx_fdm*tx_fdm'/(M);
+ tx_fdm = fdm_upconvert(tx_baseband) + 2*cos(2*pi*(0:M-1)*Fcentre/Fs);
+ tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M);
+
+ % -------------------
+ % Channel simulation
+ % -------------------
+
+ % frequency offset
- noise = Ngain/sqrt(2)*[randn(1,M) + j*randn(1,M)];
- noise_pwr = 0.9*noise_pwr + 0.1*noise*noise'/M;
for i=1:M
phase_offset *= freq_offset;
- rx_fdm(i) = phase_offset*tx_fdm(i);
+ rx_fdm(i) = phase_offset*real(tx_fdm(i));
end
+
+ % AWGN noise
+
+ noise = Ngain*randn(1,M);
+ noise_pwr = 0.9*noise_pwr + 0.1*noise*noise'/M;
rx_fdm += noise;
rx_fdm_log = [rx_fdm_log rx_fdm];
- rx_baseband = fdm_downconvert(rx_fdm);
+ % delay
+
+ rx_fdm_delay(1:Ndelay-M) = rx_fdm_delay(M+1:Ndelay);
+ rx_fdm_delay(Ndelay-M+1:Ndelay) = rx_fdm;
+
+ % -------------------
+ % Demodulator
+ % -------------------
+
+ % frequency offset estimation and correction
+
+ foff = rx_est_freq_offset(rx_fdm);
+ foff_log = [ foff_log foff ];
+ foff_rect = exp(j*2*pi*foff/Fs);
+
+ for i=1:M
+ foff_phase *= foff_rect';
+ rx_fdm_delay(i) = rx_fdm_delay(i)*foff_phase;
+ end
+
+ % baseband processing
+
+ rx_baseband = fdm_downconvert(rx_fdm_delay(1:M));
+ rx_baseband_log = [rx_baseband_log rx_baseband];
rx_filt = rx_filter(rx_baseband);
[rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband);
rx_bits = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
prev_rx_symbols = rx_symbols;
+ % count bit errors
+
[sync bit_errors] = put_test_bits(rx_bits);
if (sync == 1)
total_bit_errors = total_bit_errors + bit_errors;
end
ber = total_bit_errors/total_bits;
-printf("Eb/No: %2.2f dB %d bits %d errors Meas BER: %1.4f Theor BER: %1.4f\n", EbNo_dB,
- total_bits, total_bit_errors, ber, 0.5*erfc(sqrt(10.^(EbNo_dB/10))));
+printf("Eb/No (meas): %2.2f (%2.2f) dB %d bits %d errors QPSK BER (meas): %1.4f (%1.4f)\n",
+ EbNo_dB, 10*log10(0.25*tx_pwr*Fs/(Rs*Nc*noise_pwr)),
+ total_bits, total_bit_errors, 0.5*erfc(sqrt(10.^(EbNo_dB/10))), ber );
figure(1)
clf;
[n m] = size(rx_symbols_log);
plot(real(rx_symbols_log(:,20:m)),imag(rx_symbols_log(:,20:m)),'+')
+
figure(2)
clf;
subplot(211)
plot(rx_timing_log)
subplot(212)
-Nfft=Fs;
-S=fft(rx_fdm_log,Nfft);
-SdB=20*log10(abs(S));
+%Nfft=Fs;
+%S=fft(rx_fdm_log,Nfft);
+%SdB=20*log10(abs(S));
%plot(-Fs/2+1:Fs/2,fftshift(SdB))
-plot(SdB(1:Fs/4))
+%plot(SdB(1:Fs/4))
+plot(foff_log)
+
+
% TODO
% + handling sample slips, extra plus/minus samples
% Sensitivity to Fs
% Try BPSK
% second term of QPSK eqn
+
+% Faster sync:
+%
+% 1/ Maybe Ask Bill how we can sync in less than 10 symbols? What to
+% put in filter memories? If high SNR should sync very quickly
+% Maybe start feeding in symbols to middle of filter memory? Or do timing
+% sync before rx filtering at start?
+
projects / freetel-svn-tracking.git / commitdiff