+ does modem fall over?
[ ] try non-flat channel, e.g. 3dB difference between hi and low tones
+ make sure all estimators keep working
+[ ] test rx level sensitivity, i.e. 0 to 20dB attenuation
[ ] try to run from shell script
[ ] run a few tests
[ ] start coding in C and repeat tests
if (strcmp(modulation,'dqpsk'))
-
- if 1
% map to (Nc,1) DQPSK symbols
for c=1:Nc
tx_symbols(c) = -j*prev_tx_symbols(c);
endif
end
- else
- % map to pi/4 DQPSK from spra341 Eq (6) & (7)
-
- for c=1:Nc
-
- msb = tx_bits_matrix(c,1); lsb = tx_bits_matrix(c,2);
- a = 2*msb - 1;
- b = 2*lsb - 1;
- p = prev_tx_symbols(c);
- inphase = (real(p)*a - imag(p)*b)*0.707;
- quadrature = (imag(p)*a + real(p)*b)*0.707;
- tx_symbols(c) = inphase + j*quadrature;
- end
- end
-
else
% QPSK mapping
tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2);
% Frequency shift each modem carrier down to Nc baseband signals
-function rx_baseband = fdm_downconvert(rx_fdm)
+function rx_baseband = fdm_downconvert(rx_fdm, nin)
global Fs;
global M;
global Nc;
global phase_rx;
global freq;
- rx_baseband = zeros(1,M);
+ rx_baseband = zeros(1,nin);
% Nc/2 tones below centre freq
for c=1:Nc/2
- for i=1:M
+ for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
% Nc/2 tones above centre freq
for c=Nc/2+1:Nc
- for i=1:M
+ for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
% Pilot
c = Nc+1;
- for i=1:M
+ for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
% Receive filter each baseband signal at oversample rate P
-function rx_filt = rx_filter(rx_baseband)
+function rx_filt = rx_filter(rx_baseband, nin)
global Nc;
global M;
global P;
global gt_alpha5_root;
global Fsep;
- rx_filt = zeros(Nc+1,P);
+ rx_filt = zeros(Nc+1,nin*P/M);
% rx filter each symbol, generate P filtered output samples for each symbol.
% Note we keep memory at rate M, it's just the filter output at rate P
N=M/P;
j=1;
- for i=1:N:M
+ for i=1:N:nin
rx_filter_memory(:,Nfilter-N+1:Nfilter) = rx_baseband(:,i:i-1+N);
rx_filt(:,j) = rx_filter_memory * gt_alpha5_root';
rx_filter_memory(:,1:Nfilter-N) = rx_filter_memory(:,1+N:Nfilter);
% Estimate frequency offset of FDM signal using BPSK pilot. This is quite
% sensitive to pilot tone level wrt other carriers
-function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev)
+function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev, nin)
global M;
global Npilotbaseband;
global pilot_baseband1;
% down convert latest M samples of pilot by multiplying by
% ideal BPSK pilot signal we have generated locally
- pilot_baseband1(1:Npilotbaseband-M) = pilot_baseband1(M+1:Npilotbaseband);
- pilot_baseband2(1:Npilotbaseband-M) = pilot_baseband2(M+1:Npilotbaseband);
- for i=1:M
- pilot_baseband1(Npilotbaseband-M+i) = rx_fdm(i) * conj(pilot(i));
- pilot_baseband2(Npilotbaseband-M+i) = rx_fdm(i) * conj(pilot_prev(i));
+ pilot_baseband1(1:Npilotbaseband-nin) = pilot_baseband1(nin+1:Npilotbaseband);
+ pilot_baseband2(1:Npilotbaseband-nin) = pilot_baseband2(nin+1:Npilotbaseband);
+ for i=1:nin
+ pilot_baseband1(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot(i));
+ pilot_baseband2(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot_prev(i));
end
- [foff1 max1 pilot_lpf1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1);
- [foff2 max2 pilot_lpf2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2);
+ [foff1 max1 pilot_lpf1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1, nin);
+ [foff2 max2 pilot_lpf2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2, nin);
if max1 > max2
foff = foff1;
% LPF and peak pick part of freq est, put in a function as we call it twice
-function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf)
+function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin)
global M;
global Npilotlpf;
global Npilotcoeff;
% LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset
- pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf);
+ pilot_lpf(1:Npilotlpf-nin) = pilot_lpf(nin+1:Npilotlpf);
j = 1;
- for i = Npilotlpf-M+1:Npilotlpf
+ for i = Npilotlpf-nin+1:Npilotlpf
pilot_lpf(i) = pilot_baseband(j:j+Npilotcoeff) * pilot_coeff';
j++;
end
% Estimate optimum timing offset, and symbol receive symbols
-function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
+function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin)
global M;
global Nt;
global Nc;
global Nfiltertiming;
global gt_alpha5_root;
+ % nin adjust
+ % --------------------------------
+ % 120 -1 (one less rate P sample)
+ % 160 0 (nominal)
+ % 200 1 (one more rate P sample)
+
+ adjust = P - nin*P/M;
+
% update buffer of Nt rate P filtered symbols
- rx_filter_mem_timing(:,1:(Nt-1)*P) = rx_filter_mem_timing(:,P+1:Nt*P);
- rx_filter_mem_timing(:,(Nt-1)*P+1:Nt*P) = rx_filt(1:Nc,:);
+ rx_filter_mem_timing(:,1:(Nt-1)*P+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P);
+ %size((Nt-1)*P+1+adjust:Nt*P)
+ %size(rx_filt)
+ %adjust
+ %nin
+ rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(1:Nc,:);
% sum envelopes of all carriers
% baseband signal at rate M this enables us to resample the filtered
% rx symbol with M sample precision once we have rx_timing
- rx_baseband_mem_timing(:,1:Nfiltertiming-M) = rx_baseband_mem_timing(:,M+1:Nfiltertiming);
- rx_baseband_mem_timing(:,Nfiltertiming-M+1:Nfiltertiming) = rx_baseband;
+ rx_baseband_mem_timing(:,1:Nfiltertiming-nin) = rx_baseband_mem_timing(:,nin+1:Nfiltertiming);
+ rx_baseband_mem_timing(:,Nfiltertiming-nin+1:Nfiltertiming) = rx_baseband;
% sample right in the middle of the timing estimator window, by filtering
% at rate M
% convert symbols back to an array of bits
-function [rx_bits sync_bit] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
+function [rx_bits sync_bit f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
global Nc;
global Nb;
global Nb;
phase_difference(Nc+1) = rx_symbols(Nc+1) .* conj(prev_rx_symbols(Nc+1));
if (real(phase_difference(Nc+1)) < 0)
sync_bit = 0;
+ f_err = imag(phase_difference(Nc+1));
else
sync_bit = 1;
+ f_err = -imag(phase_difference(Nc+1));
end
else
t -= delta;
endfunction
+
+
+% freq offset state machine. Moves between acquire and track states based
+% on BPSK pilot sequence. Freq offset estimator occasionally makes mistakes
+% when used continuously. So we use it unit we have acquired the BPSK pilot,
+% then switch to a more robust tracking algorithm. If we lose sync we switch
+% back to acquire mode for fast-requisition.
+
+function [track state] = freq_state(sync_bit, state)
+
+ % acquire state, look for 6 symbol 010101 sequence from sync bit
+
+ next_state = state;
+ if state == 0
+ if sync_bit == 0
+ next_state = 1;
+ end
+ end
+ if state == 1
+ if sync_bit == 1
+ next_state = 2;
+ else
+ next_state = 0;
+ end
+ end
+ if state == 2
+ if sync_bit == 0
+ next_state = 3;
+ else
+ next_state = 0;
+ end
+ end
+ if state == 3
+ if sync_bit == 1
+ next_state = 4;
+ else
+ next_state = 0;
+ end
+ end
+ if state == 4
+ if sync_bit == 0
+ next_state = 5;
+ else
+ next_state = 0;
+ end
+ end
+ if state == 5
+ if sync_bit == 1
+ next_state = 6;
+ else
+ next_state = 0;
+ end
+ end
+
+ % states 6 and above are track mode, make sure we keep getting 0101 sync bit sequence
+
+ if state == 6
+ if sync_bit == 0
+ next_state = 7;
+ else
+ next_state = 0;
+ end
+ end
+ if state == 7
+ if sync_bit == 1
+ next_state = 6;
+ else
+ next_state = 0;
+ end
+ end
+
+ state = next_state;
+ if state >= 6
+ track = 1;
+ else
+ track = 0;
+ end
+endfunction
+
% fdmdv_demod.m
%
-% Demodulator function for FDMDV modem. Requires 48kHz sample rate raw files
+% Demodulator function for FDMDV modem. Requires 8kHz sample rate raw files
% as input
%
% Copyright David Rowe 2012
pilot_lut = generate_pilot_lut;
pilot_lut_index = 1;
- prev_pilot = zeros(1,M);
+ prev_pilot_lut_index = 3*M+1;
% BER stats
rx_symbols_log = [];
rx_timing_log = [];
foff_log = [];
+ rx_fdm_log = [];
- % resampler states
+ % misc states
- t = 3;
- ratio = 1.000;
- F=6;
- MF=M*F;
- nin = MF;
- nin_size = MF+6;
- buf_in = zeros(1,nin_size);
- rx_fdm_buf = [];
- ratio_log = [];
+ nin = M; % timing correction for sample rate differences
+ foff = 0;
+ track_log = [];
+ track = 1;
+ fest_state = 0;
% Main loop ----------------------------------------------------
for f=1:frames
- % update buf_in memory
-
- m = nin_size - nin;
- for i=1:m
- buf_in(i) = buf_in(i+nin);
- end
-
% obtain nin samples of the test input signal
-
- for i=m+1:nin_size
- buf_in(i) = fread(fin, 1, "short")/gain;
+
+ for i=1:nin
+ rx_fdm(i) = fread(fin, 1, "short")/gain;
end
- % resample at 48kHz and decimate to 8kHz
-
- [rx_fdm_mf t nin] = resample(buf_in, t, 1.0, MF);
- rx_fdm = rx_fdm_mf(1:F:MF);
-
- %rx_fdm = buf_in(m+1:m+n);
-
- %for i=1:M
- % rx_fdm(i) = fread(fin, 1, "short")/gain;
- %end
- rx_fdm_buf = [rx_fdm_buf rx_fdm];
+ rx_fdm_log = [rx_fdm_log rx_fdm(1:nin)];
% frequency offset estimation and correction
- for i=1:M
+ for i=1:nin
pilot(i) = pilot_lut(pilot_lut_index);
pilot_lut_index++;
if pilot_lut_index > 4*M
pilot_lut_index = 1;
end
+ prev_pilot(i) = pilot_lut(prev_pilot_lut_index);
+ prev_pilot_lut_index++;
+ if prev_pilot_lut_index > 4*M
+ prev_pilot_lut_index = 1;
+ end
+ end
+ foff_coarse = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin);
+ if track == 0
+ foff = foff_coarse;
end
- foff = rx_est_freq_offset(rx_fdm, pilot, prev_pilot);
- prev_pilot = pilot;
foff_log = [ foff_log foff ];
- foff = 0;
foff_rect = exp(j*2*pi*foff/Fs);
- for i=1:M
+ for i=1:nin
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_baseband = fdm_downconvert(rx_fdm, nin);
+ rx_filt = rx_filter(rx_baseband, nin);
- [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband);
+ [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin);
rx_timing_log = [rx_timing_log rx_timing];
- beta = 1E-7;
- ratio += beta*rx_timing;
- if (ratio > 1.002)
- ratio = 1.002;
+ nin = M;
+ if rx_timing > 2*M/P
+ nin += M/P;
end
- if (ratio < 0.998)
- ratio = 0.998;
+ if rx_timing < 0;
+ nin -= M/P;
end
- ratio_log = [ratio_log ratio];
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);
+ [rx_bits sync f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+ foff -= 0.5*f_err;
prev_rx_symbols = rx_symbols;
sync_log = [sync_log sync];
+ % freq est state machine
+
+ [track fest_state] = freq_state(sync, fest_state);
+ track_log = [track_log track];
+
% 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 && f > 15)
+ if (test_frame_sync == 1)
total_bit_errors = total_bit_errors + bit_errors;
total_bits = total_bits + Ntest_bits;
- bit_errors_log = [bit_errors_log bit_errors];
+ bit_errors_log = [bit_errors_log bit_errors/Ntest_bits];
+ else
+ bit_errors_log = [bit_errors_log 0];
end
% test frame sync state machine, just for more informative plots
% Plots
% ---------------------------------------------------------------------
+ xt = (1:frames)/Rs;
+ secs = frames/Rs;
+
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)),'+')
+ plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+')
axis([-2 2 -2 2]);
title('Scatter Diagram');
figure(2)
clf;
subplot(211)
- plot(rx_timing_log(15:m))
+ plot(xt, rx_timing_log)
title('timing offset (samples)');
subplot(212)
- plot(foff_log)
+ plot(xt, foff_log)
title('Freq offset (Hz)');
grid
figure(3)
clf;
subplot(211)
- %plot(rx_fdm_buf);
- %title('FDM Rx Signal');
- plot(ratio_log-1);
- title('Sampling Clock error (ppm)');
+ [a b] = size(rx_fdm_log);
+ xt1 = (1:b)/Fs;
+ plot(xt1, rx_fdm_log);
+ title('Rx FDM Signal');
subplot(212)
Nfft=Fs;
- S=fft(rx_fdm_buf,Nfft);
+ 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;
subplot(311)
- stem(sync_log)
- axis([0 frames 0 1.5]);
+ stem(xt, sync_log)
+ axis([0 secs 0 1.5]);
title('BPSK Sync')
subplot(312)
- stem(bit_errors_log);
+ stem(xt, bit_errors_log);
title('Bit Errors for test frames')
subplot(313)
- plot(test_frame_sync_log);
- axis([0 frames 0 1.5]);
+ plot(xt, test_frame_sync_log);
+ axis([0 secs 0 1.5]);
title('Test Frame Sync')
- mean(ratio_log)
+ figure(5)
+ clf;
+ plot(xt, track_log);
+ axis([0 secs 0 1.5]);
endfunction
% Simulation Parameters --------------------------------------
-frames = 50;
+frames = 100;
EbNo_dB = 7.3;
-Foff_hz = 0;
+Foff_hz = 10;
modulation = 'dqpsk';
hpa_clip = 150;
rx_bits_offset = zeros(Nc*Nb*2);
prev_tx_symbols = ones(Nc+1,1);
prev_rx_symbols = ones(Nc+1,1);
+f_err = 0;
+foff = 0;
foff_log = [];
tx_baseband_log = [];
tx_fdm_log = [];
pilot_lut = generate_pilot_lut;
pilot_lut_index = 1;
+prev_pilot_lut_index = 3*M+1;
% fixed delay simuation
% frequency offset
+ %Foff_hz += 1/Rs;
+ Foff = Foff_hz;
for i=1:M
% Time varying freq offset
- % Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs));
- % t++;
- Foff = Foff_hz;
+ %Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs));
+ %t++;
freq_offset = exp(j*2*pi*Foff/Fs);
phase_offset *= freq_offset;
tx_fdm(i) = phase_offset*tx_fdm(i);
if pilot_lut_index > 4*M
pilot_lut_index = 1;
end
+ prev_pilot(i) = pilot_lut(prev_pilot_lut_index);
+ prev_pilot_lut_index++;
+ if prev_pilot_lut_index > 4*M
+ prev_pilot_lut_index = 1;
+ end
end
- foff = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot);
- prev_pilot = pilot;
+ %foff = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M);
foff_log = [ foff_log foff ];
%foff = 0;
foff_rect = exp(j*2*pi*foff/Fs);
% baseband processing
- rx_baseband = fdm_downconvert(rx_fdm_delay(1:M));
+ rx_baseband = fdm_downconvert(rx_fdm_delay(1:M), M);
rx_baseband_log = [rx_baseband_log rx_baseband];
- rx_filt = rx_filter(rx_baseband);
+ rx_filt = rx_filter(rx_baseband, M);
- [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband);
+ [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, M);
rx_timing_log = [rx_timing_log rx_timing];
%rx_phase = rx_est_phase(rx_symbols);
else
rx_symbols_log = [rx_symbols_log rx_symbols];
endif
- [rx_bits sync] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+ [rx_bits sync f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+ foff -= 0.5*f_err;
prev_rx_symbols = rx_symbols;
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
clf;
[n m] = size(rx_symbols_log);
plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+')
-axis([-2 2 -2 2]);
+axis([-3 3 -3 3]);
title('Scatter Diagram');
figure(2)
clf;
subplot(211)
-plot(rx_timing_log)
+plot(rx_timing_log(15:m))
title('timing offset (samples)');
subplot(212)
-plot(foff_log)
+plot(foff_log(15:m))
title('Freq offset (Hz)');
figure(3)