pkg load signal;  % Octave needs this
Fs = 20000;  % sample rate
fc = Fs/4 * 1; % optional IF frequency
Rs = 2000;   % symbol rate
sps = Fs/Rs; % samples per symbol

%
% Root raised cosine pulse filter
% https://www.michael-joost.de/rrcfilter.pdf
%
r = 0.22; % bandwidth factor

ntaps = 8 * sps + 1;
st = [-floor(ntaps/2):floor(ntaps/2)] / sps;
hpulse = 1/sqrt(sps) * (sin ((1-r)*pi*st) + 4*r*st.*cos((1+r)*pi*st)) ./ (pi*st.*(1-(4*r*st).^2));

% fix the removable singularities
hpulse(ceil(ntaps/2)) = 1/sqrt(sps) * (1 - r + 4*r/pi); % t = 0 singulatiry
sing_idx = find(abs(1-(4*r*st).^2) < 0.000001);
for k = [1:length(sing_idx)]
    hpulse(sing_idx) = 1/sqrt(sps) * r/sqrt(2) * ((1+2/pi)*sin(pi/(4*r))+(1-2/pi)*cos(pi/(4*r)));
endfor

% normalize to 0 dB gain
hpulse = hpulse / sum(hpulse);
delay_hpulse = (ntaps-1)/2;

%
% Create a correlation filter for a QPSK preamble
%

% 
preamble_symbols = [3 3 3 0 0 3 0];
constellation_map = 1/sqrt(2) * [-1-1j -1+1j 1-1j 1+1j];

% Map the preamble symbols to the constellation points
preamble_qpsk = constellation_map(preamble_symbols+1);

% Upsample and pulse shape the preamble
x = sps*[upsample(preamble_qpsk, sps), zeros(1, delay_hpulse)];
preamble_bb_init = filter(hpulse, [1], x);

% Discard the most of the transient
preamble_bb = preamble_bb_init(floor(delay_hpulse*0.70)+1:end);
% Set the correlation preak ref to the middle of the last symbol
preamble_bb = preamble_bb(1:(end-ceil(sps*0.75)));

% Modulate the preamble up to IF and then create the correlation filter
k = [0 : length(preamble_bb)-1];
preamble_if = preamble_bb .* exp(1j*2*pi*fc/Fs*k);

% Use preamble_if for the "Symbols" input of the corr_est block
preamble_if

% Use this for the actual correlation filter taps in Octave/Matlab
hcorr = fliplr(conj(preamble_if));

figure(1);
t1 = [0:length(preamble_bb_init)-1]/Fs;
plot(t1, real(preamble_bb_init), 'o-', t1, imag(preamble_bb_init), 'x-');
xlabel("Time (s)");
ylabel("Amplitude");
title("Pulse Shaped Preamble Symbols - Baseband, Untrimmed");

figure(2);
t2 = [0:length(preamble_bb)-1]/Fs;
plot(t2, real(preamble_bb), 'o-', t2, imag(preamble_bb), 'x-');
xlabel("Time (s)");
ylabel("Amplitude");
title("Pulse Shaped Preamble Symbols - Baseband, Trimmed");

figure(3);
t3 = [0:length(preamble_if)-1]/Fs;
plot(t3, real(preamble_if), '.-', t3, imag(preamble_if), '.-');
xlabel("Time (s)");
ylabel("Amplitude");
title("Pulse Shaped Preamble Symbols - Modulated to IF");

figure(4);
t4 = [0:length(hcorr)-1]/Fs;
plot(t4, real(hcorr), '.-', t4, imag(hcorr), '.-');
xlabel("Time (s)");
ylabel("Amplitude");
title("Preamble Correlation Filter Taps at IF");