global Fs = 8000; % sample rate in Hz
global T = 1/Fs; % sample period in seconds
global Rs = 50; % symbol rate in Hz
-global Ts = 1/Rs; % symbol period in seconds
global Nc = 14; % number of carriers
global Nb = 2; % Bits/symbol for QPSK modulation
global Rb = Nc*Rs*Nb; % bit rate
global M = Fs/Rs; % oversampling factor
global Nsym = 4; % number of symbols to filter over
global Fsep = 75; % Separation between carriers (Hz)
-global Fcentre = 1200; % Centre frequency, Nc/2 below this, N/c above (Hz)
+global Fcentre = 1200; % Centre frequency, Nc/2 carriers below this, N/c carriers above (Hz)
global Nt = 5; % number of symbols we estimate timing over
global P = 4; % oversample factor used for rx symbol filtering
+global Nfilter = Nsym*M;
+global Nfiltertiming = M+Nfilter+M;
-% Generate root raised cosine (Root Nyquist) filter ---------------
+% root raised cosine (Root Nyquist) filter
-% thanks http://www.dsplog.com/db-install/wp-content/uploads/2008/05/raised_cosine_filter.m
+global gt_alpha5_root = gen_rn_coeffs(alpha, Rs, Nsym, M);
-alpha = 0.5; % excess bandwidth
-n = -Nsym*Ts/2:T:Nsym*Ts/2;
-global Nfilter = Nsym*M;
-global Nfiltertiming = M+Nfilter+M;
+% Initialise ----------------------------------------------------
-sincNum = sin(pi*n/Ts); % numerator of the sinc function
-sincDen = (pi*n/Ts); % denominator of the sinc function
-sincDenZero = find(abs(sincDen) < 10^-10);
-sincOp = sincNum./sincDen;
-sincOp(sincDenZero) = 1; % sin(pix/(pix) =1 for x =0
-
-cosNum = cos(alpha*pi*n/Ts);
-cosDen = (1-(2*alpha*n/Ts).^2);
-cosDenZero = find(abs(cosDen)<10^-10);
-cosOp = cosNum./cosDen;
-cosOp(cosDenZero) = pi/4;
-gt_alpha5 = sincOp.*cosOp;
-Nfft = 4096;
-GF_alpha5 = fft(gt_alpha5,Nfft)/M;
-% sqrt causes stop band to be amplifed, this hack pushes it down again
-for i=1:Nfft
- if (abs(GF_alpha5(i)) < 0.02)
- GF_alpha5(i) *= 0.001;
- endif
+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;
+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);
end
-GF_alpha5_root = sqrt(abs(GF_alpha5)) .* exp(j*angle(GF_alpha5));
-ifft_GF_alpha5_root = ifft(GF_alpha5_root);
-global gt_alpha5_root = real((ifft_GF_alpha5_root(1:Nfilter)));
+for c=Nc/2+1:Nc
+ carrier_freq = (-Nc/2 + c)*Fsep + Fcentre;
+ freq(c) = exp(j*2*pi*carrier_freq/Fs);
+end
+
+freq(Nc+1) = exp(j*2*pi*Fcentre/Fs);
+
+% Spread initial FDM carrier phase out as far as possible. This
+% helped PAPR for a few dB. We don't need to adjust rx phase as DQPSK
+% takes care of that.
+
+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 constants
+
+global Mpilotfft = 256;
+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
+
+% Freq offset estimator states
+
+global pilot_baseband1;
+global pilot_baseband2;
+pilot_baseband1 = zeros(1, Npilotbaseband); % pilot baseband samples
+pilot_baseband2 = zeros(1, Npilotbaseband); % pilot baseband samples
+global pilot_lpf1
+global pilot_lpf2
+pilot_lpf1 = zeros(1, Npilotlpf); % LPF pilot samples
+pilot_lpf2 = zeros(1, Npilotlpf); % LPF pilot samples
+
+% Timing estimator states
+
+global rx_filter_mem_timing;
+rx_filter_mem_timing = zeros(Nc+1, Nt*P);
+global rx_baseband_mem_timing;
+rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming);
+
+% Test bit stream constants
+
+global Ntest_bits = Nc*Nb*4; % length of test sequence
+global test_bits = rand(1,Ntest_bits) > 0.5;
+
+% 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);
% Functions ----------------------------------------------------
% 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, nin)
+function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin)
global M;
global Npilotlpf;
global Npilotcoeff;
endfunction
-% Estimate optimum timing offset, and symbol receive symbols
+% Estimate optimum timing offset, re-filter receive symbols
function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin)
global M;
% update buffer of Nt rate P filtered symbols
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,:);
+ rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(:,:);
% sum envelopes of all carriers
- env = sum(abs(rx_filter_mem_timing(:,:)));
+ env = sum(abs(rx_filter_mem_timing(:,:))); % use all Nc+1 carriers for timing
+ %env = abs(rx_filter_mem_timing(Nc+1,:)); % just use BPSK pilot
[n m] = size(env);
% The envelope has a frequency component at the symbol rate. The
% single DFT at frequency 2*pi/P
x = env * exp(-j*2*pi*(0:m-1)/P)';
-
+
% map phase to estimated optimum timing instant at rate M
- % the M/2 + 1 part was adjusted by experment, I know not why....
+ % the M/4 part was adjusted by experiment, I know not why....
rx_timing = angle(x)*M/(2*pi) + M/4;
if (rx_timing > M)
endfunction
-% Initialise ----------------------------------------------------
-
-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;
-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);
-end
-for c=Nc/2+1:Nc
- carrier_freq = (-Nc/2 + c)*Fsep + Fcentre;
- freq(c) = exp(j*2*pi*carrier_freq/Fs);
-end
-
-freq(Nc+1) = exp(j*2*pi*Fcentre/Fs);
-
-% Spread initial FDM carrier phase out as far as possible. This
-% helped PAPR for a few dB. We don't need to adjust rx phase as DQPSK
-% takes care of that.
-
-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 constants
-
-global Mpilotfft = 256;
-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
-
-% Freq offset estimator states constants
-
-global pilot_baseband1;
-global pilot_baseband2;
-pilot_baseband1 = zeros(1, Npilotbaseband); % pilot baseband samples
-pilot_baseband2 = zeros(1, Npilotbaseband); % pilot baseband samples
-global pilot_lpf1
-global pilot_lpf2
-pilot_lpf1 = zeros(1, Npilotlpf); % LPF pilot samples
-pilot_lpf2 = zeros(1, Npilotlpf); % LPF pilot samples
-
-% Timing estimator states
-
-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 constants
-
-global Ntest_bits = Nc*Nb*4; % length of test sequence
-global test_bits = rand(1,Ntest_bits) > 0.5;
-
-% 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
pilot_lut = [pilot_lut pilot];
end
- % discard first 4 symbols as filer memory is filling, just keep last
+ % discard first 4 symbols as filter memory is filling, just keep last
% four symbols
pilot_lut = pilot_lut(4*M+1:M*F);
frames = 100;
EbNo_dB = 7.3;
-Foff_hz = 10;
+Foff_hz = 0;
modulation = 'dqpsk';
hpa_clip = 150;
% 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 = 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
figure(2)
clf;
subplot(211)
-plot(rx_timing_log(15:m))
+plot(rx_timing_log)
title('timing offset (samples)');
subplot(212)
-plot(foff_log(15:m))
+plot(foff_log)
title('Freq offset (Hz)');
figure(3)
axis([0 frames 0 1.5]);
title('Test Frame Sync')
-% TODO
-% + handling sample slips, extra plus/minus samples
-% + simulating sample clock offsets
-% + timing, frequency offset get function
-% + memory of recent tx and rx signal for spectrogram
-% + scatter diagram get function
-
-% DPSK
-% add phase offset, frequency offset, timing offset
-% fading simulator
-% file I/O to test with Codec
-% code to measure sync recovery time
-% dump file type plotting & instrumentation
-% determine if error pattern is bursty
-% HF channel simulation
-% Offset or pi/4 QPSK and tests with real tx HPA
-% real time SNR get function
-%
-% phase estimator not working too well and would need a UW
-% to resolve ambiguity. But this is probably worth it for
-% 3dB. Test with small freq offset
-
-% Implementation loss BER issues:
-% QPSK mapping
-% Single sided noise issues
-% interference between carriers due to excess BW
-% Crappy RN coeffs
-% timing recovery off by one
-% Use a true PR sequence
-% 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?
-