% Notes for Matlab spectral analysis exercise

%----------------------------------------------
% (a) Load time series and normalize time to start at 0

  load('rf18L1.txt');
  t = rf18L1(:,1);
  t = t - t(1);
  u = rf18L1(:,2); 
  v = rf18L1(:,3); 
  w = rf18L1(:,4); 
  T = rf18L1(:,5); 
  q = rf18L1(:,6);
  p = rf18L1(:,7); 

%----------------------------------------------
%  b. Spectral analysis of w. Outputs power spectral density Pww as function
%     of frequency f (Hz).

[Pww,f] = psd(w,nfft,ns,window,noverlap,dflag);

PSD Power Spectral Density estimate.

    Pxx = PSD(X,NFFT,Fs,WINDOW) estimates the Power Spectral Density of 
    a discrete-time signal vector X using Welch's averaged, modified
    periodogram method.  
    
%----------------------------------------------
%  c. High pass filtering (passband f > 0.05) and std of w 

  cutoff = 0.05*2/rate;         % Cutoff frequency in Nyquist units
  [b, a] = butter(4,cutoff,'high'); 
  whi = filter(b,a,w);

BUTTER Butterworth digital and analog filter design.

    [B,A] = BUTTER(N,Wn) designs an Nth order lowpass digital
    Butterworth filter and returns the filter coefficients in length 
    N+1 vectors B (numerator) and A (denominator). 
 
    [B,A] = BUTTER(N,Wn,'high') designs a highpass filter.

FILTER One-dimensional digital filter.

    Y = FILTER(B,A,X) filters the data in vector X with the
    filter described by vectors A and B to create the filtered
    data Y. 

%----------------------------------------------
%  d. Integral timescale for w

% I recommend the following use of XCOV.

Tmaxlag = 2;  			% maximum lag (s)
rate = 25;                   	% Data rate (Hz)
maxlag = Tmaxlag * rate		% maximum lag (unitless)

[whi_lag,lags] = xcov(whi,maxlag,'coeff');  % autocorrelation

time_lags = lags / rate;	% lags (s)

XCOV   Cross-covariance function estimates.
 
    XCOV(A), when A is a vector, is the auto-covariance sequence.
 
    XCOV(...,MAXLAG) computes the (auto/cross) covariance over the
    range of lags:  -MAXLAG to MAXLAG, i.e., 2*MAXLAG+1 lags.
    If missing, default is MAXLAG = M-1.
 
    [C,LAGS] = XCOV(...) returns a vector of lag indices (LAGS).
    
    XCOV(...,SCALEOPT), normalizes the covariance according to SCALEOPT:
        coeff    - normalizes the sequence so that the covariances at 
                   zero lag are identically 1.0.

% integral timescale requires integration (summation)

SUM Sum of elements.

    For vectors, SUM(X) is the sum of the elements of X.

%----------------------------------------------
% e. fluxes

% useful functions: std, cov, sum, length

% EXAMPLES 

% flux (these give identical results)
	wu1 = cov(whi,uhi); 
	wu2 = sum( whi .* uhi ) / length(whi)
% variance (these give identical results)
	ww1 = cov(whi)
	ww2 = sum( whi .* whi ) / length(whi)
	ww3 = sum( whi.^2 ) / length(whi)
	wstd = std(whi)
	ww4 = wstd^2
% skewness (these give identical results)
	www =  sum( whi.^3 ) / length(whi)
	skew1 = www / ww1^(3/2)
	skew2 = www / wstd^3