% timeseries4.m
%
% Produce a time series in M chunks of size N, containing 2 sinusoidal
% signals and some noise. Save as a binary file of floats (datafile.bin).
% Takes about 40s to write 1GB (timeseries(1024,1024*256)), 2014-Jan.
%
% CCK 2017-Feb-09 ported to Matlab (now runs in both Matlab & Octave).

function timeseries(M,N, SNR)

% The optional argument SNR specifies the signal-to-noise ratio
% (default = 1.0).

if ( exist("SNR") ~= 1 ) SNR = 1.0; end

% Physical characteristics of signal
e = exp(1); % Matlab lacks this constant, which is standard in Octave.
M           % number of chunks
N           % chunk size
SNR         % signal-to-noise ratio
A = 1.0     % amplitude (V)
T = e       % period (s)
freq = 1/T      % frequency (Hz)
freq2 = e*freq  % second signal
dt = 0.25       % sampling time (s)
% Physical characteristics of noise
sigma = A / sqrt(SNR * 2)   % standard deviation (V)
mean =  0.0;                % mean (V)

fileID = fopen("datafile.bin","w")
phasenoise=0; % initialize
noise_old=0; % initialize
for i=0:(M-1)
   t = dt*( (0:(N-1)) + i*N);
   chunk = A/sqrt(2) * sin(2 * pi * freq * t); % signal 1
   % phasenoise = phasenoise(end) + cumsum( normrnd(0,freq2*dt,1,N) );
   phasenoise = phasenoise(end) + cumsum( freq2*dt*randn(1,N) );
      % phase changes in a random walk freq2*dt radians per time step. After n steps,
      % the expectation value of the phase change is <phi> = sqrt(n) * freq2*dt. 
      % The coherence time, tau, is the time at which <phi> = pi/2,
      % which leads to tau = n*dt = dt * ( pi / (2*freq2*dt) )^2 = 9.9 sec.
   chunk = chunk + A/sqrt(2) * cos(2 * pi * freq2 * t + phasenoise);
   %noise = normrnd(mean, sigma, 1, N);   % noise
   noise = mean + sigma * randn(1, N);   % noise
   noise = noise + 0.9*[noise_old(end), noise(1:end-1)]; % redden the noise
   chunk = chunk + noise;
   noise_old = noise;
   fwrite(fileID, chunk, "float", "ieee-be");
end
fclose(fileID);
