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
+
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 ---------------
% 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;
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
% 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
global Fcentre;
global freq;
global freq_rx_pilot;
- global phase_rx_pilot;
global Npilotbaseband;
global Npilotlpf;
global Npilotcoeff;
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
% 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);
% 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
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
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
% 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;
% ------------------------------------------------------------
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 = [];
% Main loop
% ---------------------------------------------------------------------
-for i=1:frames
+for f=1:frames
% -------------------
% Modulator
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
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);
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
% 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);
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];
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)
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;
projects / freetel-svn-tracking.git / commitdiff