pkg load signal;  % Octave needs this

% One of the possible encodings of the AIS preamble symbols
h_bits = [-1,-1, 1, 1, ...
          -1,-1, 1, 1, ...
          -1,-1, 1, 1, ...
          -1,-1, 1, 1, ...
          -1,-1, 1, 1, ...
          -1,-1, 1, 1, ...
          -1,-1,-1,-1,-1,-1,-1, 1].';

% AIS specific GMSK parameters
Rb = 9600;
BT = 0.4;
L = 3;
modulation_index = 1/2;

% Working sample rate
sps = 5;
Fs = Rb * sps;

% Build Gaussian pulse filter
Ls = round(L*sps);
alpha = sqrt(2/log(2)) * pi * BT;
k = [(-Ls/2+1):1:(Ls/2-1)];
taps = (erf(alpha*(k/sps + 0.5)) - erf(alpha*(k/sps - 0.5)))*0.5/sps;
K = length(taps);
if (mod(K,2) == 0)
   delay = K/2;
else
   delay = (K-1)/2;
end

% Upsample and pulse shape the correlation sequence
x = sps*[upsample(h_bits,sps); zeros(delay, size(h_bits)(2))];
h_ibaseband = filter(taps, [1], x);
% Trim the pulse shaped correlation sequence
h_baseband = h_ibaseband((delay):(end-(sps-1)),:);

% modulate the correlation sequence
fm_gain = pi/(Fs/2) * Rb/2 * modulation_index;
x = h_baseband * fm_gain;
phase = cumsum(x); % phase is integral of frequency
h_iq = exp(1i*mod(phase, 2*pi));

% create the correlation filter
h = conj(h_iq(end:-1:1,:));

figure(1);
t1 = [1:size(h_ibaseband)(1)];
plot(t1, h_ibaseband(:,1), 'x-');
title('Premable Baseband Gaussian Pulses');

figure(2);
t2 = [1:size(h_baseband)(1)];
plot(t2, h_baseband(:,1), 'x-');
title('Trimmed Premable Baseband Gaussian Pulses');

figure(3);
t3 = [1:size(h)(1)];
plot(t3, real(h(:,1)), 'x-', t3, imag(h(:,1)), 'x-');
title('Correlation Filter Taps');
h(:,1)