%% 12.714 Practical
%% Part 1a: FFT of some common functions and comparison to analtyic forms

% Start with box-car
BW = 0.5 ; AB = 1.0 ; % Width and Amplitude of box car
% Generate sample
N = 129; dt = 2*BW/N; 
BoxCar = AB*ones(1,N);
NF = 1024; % number of points in fft
FFTBoxCar = fft(BoxCar,NF)/N;
% Now rotate the fft to shift the center of the box car time to zero
toff = -floor(N/2)*dt;
f = [0:NF/2-1]/(dt*(NF)); % Freqeuncies in fft

FF = FFTBoxCar(1:NF/2).*exp(-2*pi*i*f*toff);

AnalBoxCar = 2*AB*BW*sin(2*pi*BW*f)./(2*pi*BW*f); AnalBoxCar(1)=2*AB*BW;
%
% With correct time offset, imaginary will be zero.
figure(1);
subplot(2,1,1);
plot(f,AnalBoxCar,'k',f,real(FF(1:NF/2)),'r',f,imag(FF(1:NF/2)),'g'); axis tight
lims = ylim; ylabel('BoxCar FFT'); legend('Analytic','FFT')
subplot(2,1,2);
plot(f,AnalBoxCar-real(FF(1:NF/2)),'r'); axis tight; ylim(lims); 
xlabel('Frequency (cpt)'); ylabel('Difference')

%% Part 1b: FFT of some common functions and comparison to analtyic forms

% Try a triangle function.
BW = 0.5 ; AB = 1.0 ; % Width and Amplitude of box car
% Generate sample
N =129; dt = 2*BW/N; 
BoxCar = AB*ones(1,N);
Triang = conv(BoxCar,BoxCar)/N;
NF = 1024; % number of points in fft
FFTTriang = 2*fft(Triang,NF)/N;
% Now rotate the fft to shift the time to zero
%toff = -(N-1)*dt;
toff = -floor(N-1)*dt;

f = [0:NF/2-1]/(dt*(NF));
FF = FFTTriang(1:NF/2).*exp(-2*pi*i*f*toff);

AnalTriang = 4*AB*BW*(sin(2*pi*BW*f)./(2*pi*BW*f)).^2;

figure(2);
subplot(2,1,1);
plot(f,AnalTriang,'k',f,real(FF(1:NF/2)),'r',f,imag(FF(1:NF/2)),'g'); axis tight; xlim([0 10]);
lims = ylim; ylabel('Triangular FFT'); legend('Analytic','FFT')
subplot(2,1,2);
plot(f,AnalTriang-real(FF(1:NF/2)),'r'); axis tight; ylim([-1 1]); xlim([0 10]);
xlabel('Frequency (cpt)'); ylabel('Difference')

%% Part 2a: Look at repeating effects
% Show effect of repeating of a ramp function.  This shape of spectrum is
% what happens if a slope is left in data.

gt1 = [-4:0.125:3.99]; %  Define a ramp with an even number of points.
% (Some other cases to consider: What happens if -4:0.125:4 used?
lgt1 = length(gt1); 
ff1 = [0:1/lgt1:1-1/lgt1];  ! 
gf1 = fft(gt1);
pgf1 = 10*log10(conj(gf1).*gf1);
figure(3);
subplot(3,1,1); plot(ff1,pgf1); xlim([0 1]); ylim([-20 45]); 
ylabel('Power dB 1 ramp'); xlabel('Frequency');
% Now do two ramps:
% Since there are two we divide amplitude by 2 and since twice as long in
% time we double the frequency resolution.
% Note: Here and below: line becomes dotted due to -Inf values in spectrum
% You can also see what happens when a set number of points is used in the
% FFT.
gt2 = [gt1 gt1];
gf2 = fft(gt2)/2;
ff2 = [0:1/(2*lgt1):1-1/(2*lgt1)];
pgf2 = 10*log10(conj(gf2).*gf2);
subplot(3,1,2); plot(ff2,pgf2); xlim([0 1]); ylim([-20 45]);
ylabel('Power dB 2 ramp'); xlabel('Frequency');
% With 4-ramps; 4 times amplitude and 4 times finer resolution.
gt4 = [gt2 gt2];
gf4 = fft(gt4)/4;
ff4 = [0:1/(4*lgt1):1-1/(4*lgt1)];
pgf4 = 10*log10(conj(gf4).*gf4);
subplot(3,1,3); plot(ff4,pgf4); xlim([0 1]); ylim([-20 45]);
ylabel('Power dB 4 ramp'); xlabel('Frequency');

% 
% Plot together
figure(4); hold off
plot(ff1,pgf1,'k'); hold on;  plot(ff2,pgf2,'b'); plot(ff4,pgf4,'r+'); xlim([0 1]); ylim([-20 45]); hold on
legend('1 ramp','2 ramps','4 ramps');
ylabel('Power dB'); xlabel('Frequency');

%% Part 3: Process noise example
%
% AR(1) Process; First Order Gauss Markov Process
randn('seed',130) ;  % Ensure we get same sequence each time
tau = 100; % Decay time 10 days
dt = 1;   % Data spacing in days
beta = exp(-dt/tau);
sigtot  = 2.0 ; % This is going to be sigma total
sigwn   = 1.0 ; %  White noise sigma (uncorrelated)
sigfog1 = sqrt(sigtot^2-sigwn^2) ; %  Long terms standard deviation for FOGM+WN
sigfog2 = sigtot ; % Long term sigma with just FOGM
sigep1  = sqrt(sigfog1^2*(1-beta^2)); %  Daily noise driving FOGM process
sigep2  = sqrt(sigfog2^2*(1-beta^2)); %  Daily noise driving FOGM process
%
% Generate a random sequence with 10240 samples (28 years of data)

NS = 10240*2; 
WN = randn(1,NS);
fog1 = zeros(1,NS); fog1pwn = zeros(1,NS); % Sample of FOGM model and FOGM model with white noise added
fog2 = zeros(1,NS);       % Pure FOGM
fog1(1) = WN(1)*sigfog1;  % Of this value is sample from long term behavior
fog2(1) = WN(1)*sigfog2; 
for k = 2:NS
    fog1(k) = beta*fog1(k-1) + WN(k)*sigep1;
    fog1pwn(k) = fog1(k) + randn()*sigwn ;    % Term with white noise
    fog2(k) = beta*fog2(k-1) + WN(k)*sigep2;  % Pure FOGM process
end
% Now plot results
figure(5);
subplot(2,1,1)
plot(1:NS, fog1pwn,'k'); axis tight; xlabel('Time days'); ylabel('FOGM+White Noise');
subplot(2,1,2)
plot(1:NS, fog2,'k'); axis tight; xlabel('Time days'); ylabel('FOGM only');
%
% Output the standard devations of both
fprintf('\n12.714: Practical Part 3\n')
fprintf('In FIgure 5: Standard Devation FOGM_WN %6.2f; FOGM only %6.2f mm\n',std(fog1pwn), std(fog2))
%
% Now do ACVS (autocorvaraince squence).
%
% Now take whole sequences and cmompure ACVS and PSD
FT1PWN = fft(fog1pwn)/sqrt(NS); FT2 = fft(fog2)/sqrt(NS);
fall = [0:NS-1]/(NS*dt);  % Total range (Only 1/2 is positive frequencies.
PSD1PWN = FT1PWN.*conj(FT1PWN); PSD2 = FT2.*conj(FT2);
% Compute the autocovarance sequence from the inverse FFT (See
% CDA_SL2L01.m)
ACS1 = real(ifft(PSD1PWN)); ACS2 = real(ifft(PSD2));  % See CDA_SL2L01.m
%
figure(6);  Ntau = 1024 ; % Plot only first 1024 days of ACS
subplot(2,1,1); hold off; 
plot([0:Ntau-1],ACS1([1:Ntau]),'k'); hold on; plot(0,ACS1(1),'r+');  axis tight; xlabel('Tau (days)'); ylabel('ACS FOGM+WN');
subplot(2,1,2); hold off; 
plot([0:Ntau-1],ACS2([1:Ntau]),'k'); hold on; plot(0,ACS2(1),'r+'); axis tight; xlabel('Tau (days)'); ylabel('ACS FOGs Only');
%
% Now plot the PSD's (+ve half only)
figure(7);
subplot(2,1,1); hold off; 
PSD1Th = (sigep2^2./((1-beta)^2 + 4*pi*fall.^2)+sigwn^2)/2;  % Divide by two because we have 
                                                             % folded the
                                                             % FFT spectum.
loglog(fall(1:NS/2-1),PSD1PWN(2:NS/2),'k'); axis tight
hold on; loglog(fall(1:NS/2-1),PSD1Th(2:NS/2),'r'); legend('PSD from FFT','Theory PSD','Location','SouthWest');
xlabel('Frequency (cpd)'); ylabel('PSD FOGM+WN mm^2/cpsd')
subplot(2,1,2); hold off
PSD2Th = sigep2^2./((1-beta)^2 + 4*pi*fall.^2)/2;
loglog(fall(1:NS/2-1),PSD2(2:NS/2),'k'); axis tight
hold on; loglog(fall(1:NS/2-1),PSD2Th(2:NS/2),'r');legend('PSD from FFT','Theory PSD','Location','SouthWest');
xlabel('Frequency (cpd)'); ylabel('PSD FOGM mm^2/cpsd')
%
%% Now see what we can extract by averging spectrum
% Practical Part 3b

M = 10; NG = floor(NS/M);
fals = [0:NG-1]/(NG*dt);  % Total range (Only 1/2 is positive frequencies.
PSD1K = zeros(M,NG); PSD2K = zeros(M,NG);
ACS1K = zeros(M,NG); ACS2K = zeros(M,NG);
% Split data into N groups and average
for k = 1:M
    fog1k = fog1pwn((k-1)*NG+1:k*NG);
    fog2k = fog2((k-1)*NG+1:k*NG);
%   FFT and get ACS
    FT1K = fft(fog1k)/sqrt(NG); FT2K = fft(fog2k)/sqrt(NG);

    PSD1K(k,:) = FT1K.*conj(FT1K); PSD2K(k,:) = FT2K.*conj(FT2K);
%   Compute the autocovarance sequence from the inverse FFT (See
%   CDA_SL2L01.m)
    ACS1K(k,:) = real(ifft(PSD1K(k,:))); ACS2K(k,:) = real(ifft(PSD2K(k,:)));
end
%
% Now average and plot the PSD's (+ve half only)
% Here we do a simple average using mean. Is this the best way?
PSD1ThK = (sigep2^2./((1-beta)^2 + 4*pi*fals.^2)+sigwn^2)/2;
PSD2ThK = sigep2^2./((1-beta)^2 + 4*pi*fals.^2)/2;

figure(8);
subplot(2,1,1); hold off; 
loglog(fals(1:NG/2-1),mean(PSD1K(:,2:NG/2)),'k'); axis tight
hold on; loglog(fals(1:NG/2-1),PSD1ThK(2:NG/2),'r'); legend('Mean PSD from FFT','Theory PSD','Location','SouthWest');
xlabel('Frequency (cpd)'); ylabel('PSD FOGM+WN mm^2/cpsd')
subplot(2,1,2); hold off
loglog(fals(1:NG/2-1),mean(PSD2K(:,2:NG/2)),'k'); axis tight
hold on; loglog(fals(1:NG/2-1),PSD2ThK(2:NG/2),'r');legend('Mean PSD from FFT','Theory PSD','Location','SouthWest');
xlabel('Frequency (cpd)'); ylabel('PSD FOGM mm^2/cpsd')