% Create a chirped waveform meant to roughly resemble that from GW150914. I started
% by playing around with Graphing Calculator (chirp.gcf).
%
% CCK 2015-Feb-18
%
function [time, data, waveform, frequency, noise] = chirp(stdev, ar1)
%
% Arguments: 
%  stdev = standard deviation of the noise.
%  ar1 = lag-1 autocorrelation of the noise.
%
% Outputs:
%  time = time axis, in seconds.
%  data = chirped waveform plus (red) noise.
%  waveform = just the waveform, without noise.
%  frequency = time-varying frequency of the waveform.
%  noise = just the noise, equivalent to data - waveform.
%
% Example:
%  [time, data, waveform, frequency, noise] = chirp(0.15,0.9);
% Above result will loosely resemble plots in figure 1 of the GW150914 discovery paper.

% The "time" scale of my analytic chirp envelope is of order 2, while the GW data lasts ~0.1 s.
% So I'm going to begin with a time axis in seconds and then convert to x, a time in units of 0.05s
% that can be plugged into my expressions.
%
% CCK 2018-Feb-21 updated to run in Matlab
time = -1:0.0001:1.5;  % time axis
conversion = 20; % conversion factor from time in seconds to pseudo-time x
x = conversion*(time-0.43);
envelope = 1.5 ./ (exp(-x) + exp(10*x));
freq = x + 20*(1+2*exp(2*x)); % in units to match x
frequency = conversion*freq/(2*pi); % frequency in Hz
chirp = sin(freq .* x);
waveform = chirp .* envelope;

figure(1)
hold off
plot([time,fliplr(time)], [envelope,fliplr(-envelope)], '--b', 'linewidth',2.0)
xlabel("time (s)")
hold on
plot(time, chirp, '-c', 'linewidth',2.0)
plot(time, log2(frequency/128),'-g', 'linewidth',2.0)
plot(time, waveform, '-r', 'linewidth',4.0)
axis([0.2,0.5,-1.5,2.5])
legend("envelope","chirp","log_2(f)-7 (Hz)","waveform")
legend("boxoff")
legend("location","northwest")
print("chirp1.pdf","-dpdfwrite")

% Now create red noise specified by stdev and sigma. Add to the waveform.
noise(1)=stdev*randn(1);
for i=2:numel(time)
   noise(i) = ar1*noise(i-1) + stdev*sqrt(1-ar1^2)*randn(1);
end
data = waveform + noise;

figure(2)
hold off
% plot(time, data, 'r', 'linewidth', 2.0, time, noise, 'b', 'linewidth',1.0) % does not work in Matlab
plot(time, data, 'r', 'linewidth', 2.0)
hold on
plot(time, noise, 'b', 'linewidth',1.0)
xlabel("time (s)")
legend("data","noise")
legend("boxoff")
legend("location","northwest")
axis([min(time),max(time),-2,2])
print("chirp2long.pdf","-dpdfwrite")
axis([0.25,0.47,-2,2])
print("chirp2.pdf","-dpdfwrite")
