%
% MatLab program to generate filter coefficients for WWV/H clock driver
%
fn = 4000; nb = (2700 - 50) / fn;
fprintf(1, 'sample frequency %5.0f, noise bandwidth %.0f\n', 2 * fn, fn * nb)
fid = fopen('table.c', 'w');
%
% Construct sin table 100 Hz at 125-us increments
%
fprintf(fid, 'sin table 100 Hz\n');
t = linspace(0, 7999, 8000) / 8000;
b = sin(100 * t(1:1 + 10 * 8) * (2 * pi)) / pi;
for i = 1:length(b)
        fprintf(fid, ' %13.6e,', b(i));
end
%
% Construct noise sample bandlimited in 2400-Hz
%
n = 60; f = [0 300 600 2700 3000 4000] / fn; m = [0 0 1 1 0 0];
b = remez(n, f, m);
[h, w] = freqz(b, 1, 256); w = w * fn / pi; h = 20 * log10(abs(h));
r = filter(b, 1, randn(1, 8000)); r = r / std(r);
%
% The WWV sync pulse is five cycles at 1000 Hz; The WWVH sync pulse is
% six cycles at 1200 Hz.
%
sync1 = cos(1000 * t(1:1 + 5 * 8) * 2 * pi);
sig1 = linspace(0, 0, 8000);
sig1(1001:1000 + length(sync1)) = sync1;
sync2 = cos(1200 * t(1:1 + 5 * 8) * 2 * pi);
sig2 = linspace(0, 0, 8000);
sig2(1001:1000 + length(sync2)) = sync2;
sync3(1:1 + 5 * 8) = 2;
%
% Construct the modulated waveform
%
mod = sqrt(2) * sin(100 * t(1:1 + 170 * 8) * (2 * pi));
sig3 = linspace(0, 0, 8000);
sig3(1001 + (200 * 8): 1000 + (200 * 8) + length(mod)) = mod;
%
% Construct 600-Hz bandpass filter at 1100 Hz for WWV/WWVH
%
fprintf(fid, '600-Hz bandpass (b, a)\n');
n = 4; f = [800 1400] / fn;
[b a] = ellip(n, .2, 50, f);
for i = 1:length(b)
        fprintf(fid, '%d %13.6e %13.6e\n', i - 1, b(i), a(i));
end
[h, w] = freqz(b, a, 256);
w = w * fn / pi;
fprintf(1, 'bpf processing gain %5.1f\n', 20 * log10(nb * 256 / sum(abs(h))))
h = 20 * log10(abs(h));
plot(w, h);
set(gca, 'FontSize', 16, 'Xlim', [0 2500], 'Ylim', [-60 10])
xlabel('Frequency (Hz)');
ylabel('Response (dB)');
print -deps bpf1100;
%
% Construct matched filter for five cycles at 1000 Hz
%
fprintf(fid, '5-ms matched filter at 1000 Hz\n');
c1 = sqrt(3) * sync1 / length(sync1);
for i = 1:length(c1)
        fprintf(fid, '%d %13.6e\n ', i - 1, c1(i));
end
y1 = xcorr(sync1, sync1) * pi / 4;
plot(y1)
set(gca, 'FontSize', 16, 'xlim', [1 80])
xlabel('A/D Samples');
ylabel('Amplitude');
print -deps wwv;
x = 20 * log10(max(y1));
fprintf(1, 'sync proc gain %.1f\n', x)
%
% Construct matched filter for six cycles at 1200 Hz
%
fprintf(fid, '5-ms matched filter at 1200 Hz\n');
c2 = sqrt(3) * sync2 / length(sync2);
for i = 1:length(c2)
        fprintf(fid, '%d %13.6e\n', i-1, c2(i));
end
y2 = xcorr(sync1, sync2) * pi / 4;
x = 20 * log10(max(y1) / max(abs(y2)));
fprintf(1, 'cross proc gain %.1f\n', x)
plot(y2)
set(gca, 'FontSize', 16, 'xlim', [1 80])
xlabel('A/D Samples');
ylabel('Amplitude');
print -deps wwvh;
%
% The data bits are sent as width-modulated pulses of 170, 470 and 770
% ms duration. These are extracted using a 150-Hz lowpass filter and
% 170-ms matched filter.
%
% Construct 150-Hz lowpass filter (b, a)
%
fprintf(fid, '150-Hz lowpass (b, a)\n');
n = 4; f = 150 / fn;
[b a] = ellip(n, .2, 50, f);
for i = 1:length(b)
        fprintf(fid, '%d %13.6e %13.6e\n ', i - 1, b(i), a(i));
end
[h, w] = freqz(b, a, 256);
w = w * fn / pi;
fprintf(1, 'lpf proc gain %5.1f\n', 20 * log10(nb * 256 / sum(abs(h))))
h = 20 * log10(abs(h));
plot(w, h);
set(gca, 'FontSize', 16, 'Xlim', [0 700], 'Ylim', [-60 10])
xlabel('Frequency (Hz)');
ylabel('Response (dB)');
print -deps lpf150;
%
% Construct matched filter for 17 cycles at 100 Hz
%
data(1:170 * 8) = 1;
y3 = xcorr(data, data) * pi / 4;
x = 10 * log10(max(y3));
fprintf(1, 'mod proc gain %.1f\n', x)

fclose(fid);