% WWV Test Signal Generator
% KD8OXT 18 July 2021
% Generating an audio file for WWV/H scientific modulation test.

% Housekeeping items:
close all; clear all; clc
% If desired, reset RNG so numbers are reproducible:
rng default

% Time specifications:
   Fs = 48000;                   % samples per second
   dt = 1/Fs;                   % seconds per sample
   StopTime = 2;                % seconds
   t = (0:dt:StopTime-dt)';     % seconds


% Generate white noise:
y1 = wgn(Fs*StopTime,1,0);
% normalize to prevent clipping:
y1 = y1/max(abs(y1));
fprintf(strcat("White noise: ",  num2str(length(y1)/Fs), ' seconds.\n'))

% Phase-locked audio signals:
% Tones of 2 kHz, 3 kHz, 4 kHz and 5 kHz (they said that they have 
% transmitted 5 kHz in the past). All tones on at the same time and 
% the same amplitude, then down by 3 dB, then down by 10 dB.

t = (0:dt:1-dt)';     % 1 second
y2a= cos(2*pi*2000*t); % 2 kHz
y2b= cos(2*pi*3000*t); % 3 kHz
y2c= cos(2*pi*4000*t); % 4 kHz
y2d= cos(2*pi*5000*t); % 5 kHz

y2_tones = y2a+y2b+y2c+y2d;
y2_tones = .25*y2_tones; %avoid clipping

% attenuate by 3 dB, down to ~-30 dB:
y2=y2_tones;
y_attenuated = y2; y_prev=y2;
for i = 1:9
    y_attenuated = 1/sqrt(2)*y_prev;
    y_prev = y_attenuated;
    y2 = [y2; y_attenuated];
end
fprintf(strcat("Coherent tones: ",  num2str(length(y2)/Fs), ' seconds.\n'))

% Let's see the final results:
figure, subplot(2, 2, 1)
hold on
plot(y2a), plot(y2b), plot(y2c), plot(y2d), plot(y2, 'LineWidth', 4)
legend('2 kHz', '3 kHz', '4 kHz', '5 kHz', 'Scaled Sum')
ax = gca;
ax.XLim = [0 120]; 
ax.YLim = [-1 1];
title('Coherent Tones')
% Plot our attenuation:
    subplot(2, 2, 2)
    time = (0:dt:size(y2, 1)/Fs-dt)';     % seconds
    semilogy(time, y2)
   xlabel('time (in seconds)');
   title('Attenuation of Coherent Tones');

% Chirp signal:
longchirp = 1; % long chirp: 1 second
t2 = (0:dt:longchirp-dt)'; 
% linear chirps between 0 and 5 kHz
y3 = chirp(t2,0,longchirp,5000); % up
y4 = chirp(t2,5000,longchirp,0); % down

% Short chirp:
shortchirp = .05; % short chirp: .05 second
t3 = (0:dt:shortchirp-dt)'; 
% linear chirps between 0 and 5 kHz
y5 = chirp(t3,0,shortchirp,5000); % up
y6 = chirp(t3,5000,shortchirp,0); % down

% compile all the chirps:
short_up = [y5; zeros(.1*Fs, 1); y5; zeros(.1*Fs, 1); y5; zeros(.1*Fs, 1)];
short_down = [y6; zeros(.1*Fs, 1); y6; zeros(.1*Fs, 1); y6; zeros(.1*Fs, 1)];
long_up = [y3; zeros(.1*Fs, 1); y3; zeros(.1*Fs, 1); y3; zeros(.1*Fs, 1)];
long_down = [y4; zeros(.1*Fs, 1); y4; zeros(.1*Fs, 1); y4; zeros(.1* Fs, 1)];
chirps = [short_up; short_down; zeros(.5*Fs, 1); long_up; long_down];

% End the chirp sequence on an even second:
fprintf(strcat("Adding ", num2str(1-mod(size(chirps, 1)/Fs, 1)), ' seconds of blank time to end on an even second...'))
chirps = [chirps; zeros(Fs - mod(size(chirps, 1), Fs), 1)];
fprintf(strcat(" The chirps take ",  num2str(length(chirps)/Fs), ' seconds total. \n'))

% Plot chirps:
subplot(2,2, 3), pspectrum(chirps, length(chirps), 'spectrogram')
title('Chirp Sequence')

% One cycle of 2.5 kHz:
t4 = (0:dt:1/2500);
spike = sin(2*pi*2500*t4);
subplot(2, 4, 7), plot(t4, spike), title('One cycle of 2.5 kHz')
y7 = zeros(Fs, 1); % preallocate a second to add spikes into
spikecount = 6; % number of spikes to put in, -1
for n = Fs/spikecount:Fs/spikecount:Fs-size(spike, 2)
    y7(n:n+size(spike, 2)-1) = spike;
end

% And another round, this time for 5 kHz:
t5 = (0:dt:1/5000);
spike = sin(2*pi*5000*t5);
subplot(2, 4, 8), plot(t5, spike), title('One cycle of 5 kHz')
y7b = zeros(Fs, 1); % preallocate a second to add spikes into
spikecount = 6; % number of spikes to put in, -1
for n = Fs/spikecount:Fs/spikecount:Fs-size(spike, 2)
    y7b(n:n+size(spike, 2)-1) = spike;
end

% Append the two:
y7 = [y7; y7b];
fprintf(strcat("Pulses: ",  num2str(length(y7)/Fs), ' seconds.\n'))

% Compile into one waveform:
y=[y1; zeros(Fs, 1); y2; zeros(Fs, 1); chirps; zeros(Fs, 1); y7; zeros(Fs, 1); y1; zeros(3*Fs, 1)];

% If desired, add in voice clip:
audiowrite("test-without-voice.wav", y, Fs) %default Fs=8192
voice = zeros(Fs*10, 1); % preallocate 10 seconds
voice_audio = audioread('voiceintro.wav');
voice(1:size(voice_audio, 1)) = voice_audio;
voice = voice/max(abs(voice)); %add amplitude
y = [voice; y];
fprintf(strcat("Voice clip: ",  num2str(length(voice)/Fs), ' seconds.\n'))
fprintf(strcat("Total length: ",  num2str(length(y)/Fs), ' seconds.\n'))

% Plot the signal versus time:
   figure;
   time = (0:dt:size(y, 1)/Fs-dt)';     % seconds
   plot(time, y)
   xlabel('time (in seconds)');
   title('Test Signal');
   zoom xon;
 

% Combine the two into one <45-second clip:
audiowrite("test.wav", y, Fs) %default Fs=8192
writematrix(y, 'test.csv')

% Play clip:
sound(y, Fs);
