% fdmdv.m
%
% Functions that implement a Frequency Divison Multiplexed Modem for
% Digital Voice (FDMDV) over HF channels.
%
% Copyright David Rowe 2012
% This program is distributed under the terms of the GNU General Public License
% 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
global T = 1/Fs; % sample period in seconds
global Rs = 50; % symbol rate in Hz
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 = 6; % number of symbols to filter over
global Fsep = 75; % Separation between carriers (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;
alpha = 0.5;
global snr_coeff = 0.9 % SNR est averaging filter coeff
% root raised cosine (Root Nyquist) filter
global gt_alpha5_root;
gt_alpha5_root = gen_rn_coeffs(alpha, T, Rs, Nsym, M);
% Functions ----------------------------------------------------
% generate Nc+1 QPSK symbols from vector of (1,Nc*Nb) input bits. The
% Nc+1 symbol is the +1 -1 +1 .... BPSK sync carrier
function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
global Nc;
global Nb;
global pilot_bit;
% re-arrange as (Nc,Nb) matrix
tx_bits_matrix = zeros(Nc,Nb);
tx_bits_matrix(1:Nc,1) = tx_bits(1:Nb:Nb*Nc);
tx_bits_matrix(1:Nc,2) = tx_bits(2:Nb:Nb*Nc);
if (strcmp(modulation,'dqpsk'))
% map to (Nc,1) DQPSK symbols
for c=1:Nc
msb = tx_bits_matrix(c,1); lsb = tx_bits_matrix(c,2);
if ((msb == 0) && (lsb == 0))
tx_symbols(c) = prev_tx_symbols(c);
endif
if ((msb == 0) && (lsb == 1))
tx_symbols(c) = j*prev_tx_symbols(c);
endif
if ((msb == 1) && (lsb == 0))
tx_symbols(c) = -prev_tx_symbols(c);
endif
if ((msb == 1) && (lsb == 1))
tx_symbols(c) = -j*prev_tx_symbols(c);
endif
end
else
% QPSK mapping
tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2);
endif
% +1 -1 +1 -1 BPSK sync carrier, once filtered becomes two spectral
% lines at +/- Rs/2
if pilot_bit
tx_symbols(Nc+1) = -prev_tx_symbols(Nc+1);
else
tx_symbols(Nc+1) = prev_tx_symbols(Nc+1);
end
if pilot_bit
pilot_bit = 0;
else
pilot_bit = 1;
end
endfunction
% Given Nc*Nb bits construct M samples (1 symbol) of Nc filtered
% symbols streams
function tx_baseband = tx_filter(tx_symbols)
global Nc;
global M;
global tx_filter_memory;
global Nfilter;
global gt_alpha5_root;
tx_baseband = zeros(Nc+1,M);
% tx filter each symbol, generate M filtered output samples for each symbol.
% Efficient polyphase filter techniques used as tx_filter_memory is sparse
tx_filter_memory(:,Nfilter) = sqrt(2)/2*tx_symbols;
for i=1:M
tx_baseband(:,i) = M*tx_filter_memory(:,M:M:Nfilter) * gt_alpha5_root(M-i+1:M:Nfilter)';
end
tx_filter_memory(:,1:Nfilter-M) = tx_filter_memory(:,M+1:Nfilter);
tx_filter_memory(:,Nfilter-M+1:Nfilter) = zeros(Nc+1,M);
endfunction
% Construct FDM signal by frequency shifting each filtered symbol
% stream. Returns complex signal so we can apply frequency offsets
% easily.
function tx_fdm = fdm_upconvert(tx_filt)
global Fs;
global M;
global Nc;
global Fsep;
global phase_tx;
global freq;
tx_fdm = zeros(1,M);
% Nc/2 tones below centre freq
for c=1:Nc/2
for i=1:M
phase_tx(c) = phase_tx(c) * freq(c);
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
for i=1:M
phase_tx(c) = phase_tx(c) * freq(c);
tx_fdm(i) = tx_fdm(i) + tx_filt(c,i)*phase_tx(c);
end
end
% add centre pilot tone
c = Nc+1;
for i=1:M
phase_tx(c) = phase_tx(c) * freq(c);
pilot(i) = 2*tx_filt(c,i)*phase_tx(c);
tx_fdm(i) = tx_fdm(i) + pilot(i);
end
% Scale such that total Carrier power C of real(tx_fdm) = Nc. This
% excludes the power of the pilot tone.
% We return the complex (single sided) signal to make frequency
% shifting for the purpose of testing easier
tx_fdm = 2*tx_fdm;
endfunction
% Frequency shift each modem carrier down to Nc+1 baseband signals
function rx_baseband = fdm_downconvert(rx_fdm, nin)
global Fs;
global M;
global Nc;
global Fsep;
global phase_rx;
global freq;
rx_baseband = zeros(1,nin);
% Nc/2 tones below centre freq
for c=1:Nc/2
for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
end
% Nc/2 tones above centre freq
for c=Nc/2+1:Nc
for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
end
% Pilot
c = Nc+1;
for i=1:nin
phase_rx(c) = phase_rx(c) * freq(c);
rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
end
endfunction
% Receive filter each baseband signal at oversample rate P
function rx_filt = rx_filter(rx_baseband, nin)
global Nc;
global M;
global P;
global rx_filter_memory;
global Nfilter;
global gt_alpha5_root;
global Fsep;
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: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);
j+=1;
end
endfunction
% LPF and peak pick part of freq est, put in a function as we call it twice
function [foff imax pilot_lpf S] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin)
global M;
global Npilotlpf;
global Npilotcoeff;
global Fs;
global Mpilotfft;
global pilot_coeff;
% LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset
pilot_lpf(1:Npilotlpf-nin) = pilot_lpf(nin+1:Npilotlpf);
j = 1;
for i = Npilotlpf-nin+1:Npilotlpf
pilot_lpf(i) = pilot_baseband(j:j+Npilotcoeff-1) * 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
h = hanning(Npilotlpf);
s = pilot_lpf(1:Mpilot:Npilotlpf) .* h(1:Mpilot:Npilotlpf)';
s = [s zeros(1,Mpilotfft-Npilotlpf/Mpilot)];
S = fft(s, Mpilotfft);
% peak pick and convert to Hz
[imax 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 frequency offset of FDM signal using BPSK pilot. This is quite
% sensitive to pilot tone level wrt other carriers
function [foff S1 S2] = rx_est_freq_offset(rx_fdm, pilot, pilot_prev, nin)
global M;
global Npilotbaseband;
global pilot_baseband1;
global pilot_baseband2;
global pilot_lpf1;
global pilot_lpf2;
% down convert latest nin samples of pilot by multiplying by
% ideal BPSK pilot signal we have generated locally. This
% peak of the resulting signal is sensitive to the time shift between
% the received and local version of the pilot, so we do it twice at
% different time shifts and choose the maximum.
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 S1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1, nin);
[foff2 max2 pilot_lpf2 S2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2, nin);
if max1 > max2
foff = foff1;
else
foff = foff2;
end
endfunction
% Estimate optimum timing offset, re-filter receive symbols
function [rx_symbols rx_timing env] = rx_est_timing(rx_filt, rx_baseband, nin)
global M;
global Nt;
global Nc;
global rx_filter_mem_timing;
global rx_baseband_mem_timing;
global P;
global Nfilter;
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+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P);
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(:,:))); % 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
% phase of this frequency component indicates the timing. So work out
% 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/4 part was adjusted by experiment, I know not why....
rx_timing = angle(x)*M/(2*pi) + M/4;
if (rx_timing > M)
rx_timing -= M;
end
if (rx_timing < -M)
rx_timing += M;
end
% rx_baseband_mem_timing contains M + Nfilter + M samples of the
% 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-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
s = round(rx_timing) + M;
rx_symbols = rx_baseband_mem_timing(:,s+1:s+Nfilter) * gt_alpha5_root';
endfunction
% Phase estimation function - probably won't work over a HF channel
% Tries to operate over a single symbol but uses phase information from
% all Nc carriers which should increase the SNR of phase estimate.
% 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.
function rx_phase = rx_est_phase(prev_rx_symbols, rx_symbols)
% modulation strip
rx_phase = angle(sum(rx_symbols .^ 4))/4;
endfunction
% convert symbols back to an array of bits
function [rx_bits sync_bit f_err phase_difference] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
global Nc;
global Nb;
global Nb;
if (strcmp(modulation,'dqpsk'))
% 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
for c=1:Nc
d = phase_difference(c);
if ((real(d) >= 0) && (imag(d) >= 0))
msb = 0; lsb = 0;
endif
if ((real(d) < 0) && (imag(d) >= 0))
msb = 0; lsb = 1;
endif
if ((real(d) < 0) && (imag(d) < 0))
msb = 1; lsb = 0;
endif
if ((real(d) >= 0) && (imag(d) < 0))
msb = 1; lsb = 1;
endif
rx_bits(2*(c-1)+1) = msb;
rx_bits(2*(c-1)+2) = lsb;
end
% Extract DBPSK encoded Sync bit
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));
else
sync_bit = 0;
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
rx_bits(1:Nb:Nc*Nb) = real(rx_symbols) > 0;
rx_bits(2:Nb:Nc*Nb) = imag(rx_symbols) > 0;
endif
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 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
% Accepts nbits from rx and attempts to sync with test_bits sequence.
% if sync OK measures bit errors
function [sync bit_errors] = put_test_bits(rx_bits)
global Ntest_bits; % length of test sequence
global test_bits;
global rx_test_bits_mem;
% Append to our memory
[m n] = size(rx_bits);
rx_test_bits_mem(1:Ntest_bits-n) = rx_test_bits_mem(n+1:Ntest_bits);
rx_test_bits_mem(Ntest_bits-n+1:Ntest_bits) = rx_bits;
% see how many bit errors we get when checked against test sequence
bit_errors = sum(xor(test_bits,rx_test_bits_mem));
% if less than a thresh we are aligned and in sync with test sequence
ber = bit_errors/Ntest_bits;
sync = 0;
if (ber < 0.2)
sync = 1;
endif
endfunction
% Generate M samples of DBPSK pilot signal for Freq offset estimation
function [pilot_fdm bit symbol filter_mem phase] = generate_pilot_fdm(bit, symbol, filter_mem, phase, freq)
global M;
global Nfilter;
global gt_alpha5_root;
% +1 -1 +1 -1 DBPSK sync carrier, once filtered becomes two spectral
% lines at +/- Rs/2
if bit
symbol = -symbol;
else
symbol = symbol;
end
if bit
bit = 0;
else
bit = 1;
end
% filter DPSK symbol to create M baseband samples
filter_mem(Nfilter) = (sqrt(2)/2)*symbol;
for i=1:M
tx_baseband(i) = M*filter_mem(M:M:Nfilter) * gt_alpha5_root(M-i+1:M:Nfilter)';
end
filter_mem(1:Nfilter-M) = filter_mem(M+1:Nfilter);
filter_mem(Nfilter-M+1:Nfilter) = zeros(1,M);
% upconvert
for i=1:M
phase = phase * freq;
pilot_fdm(i) = sqrt(2)*2*tx_baseband(i)*phase;
end
endfunction
% Generate a 4M sample vector of DBPSK pilot signal. As the pilot signal
% is periodic in 4M samples we can then use this vector as a look up table
% for pilot signal generation in the demod.
function pilot_lut = generate_pilot_lut()
global Nc;
global Nfilter;
global M;
global freq;
% pilot states
pilot_rx_bit = 0;
pilot_symbol = sqrt(2);
pilot_freq = freq(Nc+1);
pilot_phase = 1;
pilot_filter_mem = zeros(1, Nfilter);
%prev_pilot = zeros(M,1);
pilot_lut = [];
F=8;
for f=1:F
[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);
%prev_pilot = pilot;
pilot_lut = [pilot_lut pilot];
end
% discard first 4 symbols as filter memory is filling, just keep last
% four symbols
pilot_lut = pilot_lut(4*M+1:M*F);
endfunction
% grab next pilot samples for freq offset estimation at demod
function [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, nin)
global M;
global pilot_lut;
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
endfunction
% Change the sample rate by a small amount, for example 1000ppm (ratio
% = 1.001). Always returns nout samples in buf_out, but uses a
% variable number of input samples nin to accomodate the change in
% sample rate. nin is nominally set to nout, but may use nout +/- 2
% samples to accomodate the different sample rates. buf_in should be
% of length nout+6 samples to accomodate this, and buf_in should be
% updated externally based on the nin returned each time. "ratio" is
% Fs_in/Fs_out, for example 48048/48000 = 1.001 (+1000ppm) or
% 47952/48000 = 0.999 (-1000ppm). Uses linear interpolation to
% perform the resampling. This requires a highly over-sampled signal,
% for example 48000Hz sample rate for the modem signal centred on
% 1kHz, otherwise linear interpolation will have a low pass filter effect
% (for example an 8000Hz sample rate for modem signal centred on 1kHz
% would cause problems).
function [buf_out t nin] = resample(buf_in, t, ratio, nout)
for i=1:nout
c = floor(t);
a = t - c;
b = 1 - a;
buf_out(i) = buf_in(c)*b + buf_in(c+1)*a;
t += ratio;
end
t -= nout;
% adjust nin and t so that on next call we start with 3 < t < 4,
% this gives us +/- 2 samples room to move before we hit start or
% end of buf_in
delta = floor(t - 3);
nin = nout + delta;
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 until 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
% Save test bits to a text file in the form of a C array
function test_bits_file(filename)
global test_bits;
global Ntest_bits;
f=fopen(filename,"wt");
fprintf(f,"/* Generated by test_bits_file() Octave function */\n\n");
fprintf(f,"const int test_bits[]={\n");
for m=1:Ntest_bits-1
fprintf(f," %d,\n",test_bits(m));
endfor
fprintf(f," %d\n};\n",test_bits(Ntest_bits));
fclose(f);
endfunction
% Saves RN filter coeffs to a text file in the form of a C array
function rn_file(filename)
global gt_alpha5_root;
global Nfilter;
f=fopen(filename,"wt");
fprintf(f,"/* Generated by rn_file() Octave function */\n\n");
fprintf(f,"const float gt_alpha5_root[]={\n");
for m=1:Nfilter-1
fprintf(f," %g,\n",gt_alpha5_root(m));
endfor
fprintf(f," %g\n};\n",gt_alpha5_root(Nfilter));
fclose(f);
endfunction
function pilot_coeff_file(filename)
global pilot_coeff;
global Npilotcoeff;
f=fopen(filename,"wt");
fprintf(f,"/* Generated by pilot_coeff_file() Octave function */\n\n");
fprintf(f,"const float pilot_coeff[]={\n");
for m=1:Npilotcoeff-1
fprintf(f," %g,\n",pilot_coeff(m));
endfor
fprintf(f," %g\n};\n",pilot_coeff(Npilotcoeff));
fclose(f);
endfunction
% Saves hanning window coeffs to a text file in the form of a C array
function hanning_file(filename)
global Npilotlpf;
h = hanning(Npilotlpf);
f=fopen(filename,"wt");
fprintf(f,"/* Generated by hanning_file() Octave function */\n\n");
fprintf(f,"const float hanning[]={\n");
for m=1:Npilotlpf-1
fprintf(f," %g,\n", h(m));
endfor
fprintf(f," %g\n};\n", h(Npilotlpf));
fclose(f);
endfunction
function png_file(fig, pngfilename)
figure(fig);
pngname = sprintf("%s.png",pngfilename);
print(pngname, '-dpng', "-S500,500")
pngname = sprintf("%s_large.png",pngfilename);
print(pngname, '-dpng', "-S800,600")
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-1, 200/(Fs/2))'; % 200Hz LPF
global Npilotbaseband = Npilotcoeff + M + M/P; % number of pilot baseband samples reqd for pilot LPF
global Npilotlpf = 4*M; % number of samples we DFT pilot over, pilot est window
% pilot LUT, used for copy of pilot at rx
global pilot_lut;
pilot_lut = generate_pilot_lut();
pilot_lut_index = 1;
prev_pilot_lut_index = 3*M+1;
% 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);