function res = calculate_spectra(varargin)
% Function for calculating power spectra and cospectra from EC raw data
% INPUT:  x =  variable 1, vector (typically w)
%         y =  variable 2, vector
%         m = length of period to be analyzed, must be a power of two!
%             (m=2^14 for 30min data)
%         fs = sampling freq (Hz)
%         lag = lag time vector given as number of data points
% OUTPUT: freq = frequency vector[Hz]
%         E1 = power spectral energy of variable 1
%         E2 = power spectral energy of variable 2
%         Ec = co-spectral energy of variables 1 and 2
% A.Nordbo 01/2010, originally from U.Rannik,I.Mammarella and S.Launiainen
%
% Used in Peltola et al. (2020): "On the 
% high frequency response correction of eddy covariance fluxes. Part 2: the
% experimental approach and the interdependence with the time lag
% estimation" submitted to Atmospheric Measurement Techniques.
%
% Last update (09/2020) by O. Peltola
% olli.peltola@fmi.fi



    args = varargin;

    x1 = args{1};
    y1 = args{2};
    m = args{3};
    fs = args{4};
    
    n = max(size(x1));	% Number of data points
    
    % fixing the lag between time series
    if nargin > 4
        lag = args{5};
        if lag >= 0
            x = x1(1:n-lag);
            y = y1(1+lag:n);
        else
            lag = -lag;
            x = x1(1+lag:n);
            y = y1(1:n-lag);
        end
    else
        x = x1;
        y = y1;
    end
    
    x = x(:);
    y = y(:);
    
    if size(x,1)<m
        x = [x;x(1:(m-size(x,1)),1)];
        y = [y;y(1:(m-size(y,1)),1)];
    end
    
    index = 1:m;
    
    % detrend
    x = detrend(x(index),'constant');
    y = detrend(y(index),'constant');

    
    % fourier transforms
    Xx = fft(x);
    Yy = fft(y);
    
    Pyy = abs(Yy).^2;           % spectral energy = real(Yy).^2 + imag(Yy).^2 = Yy.*conj(Yy)
    % Pxx = abs(Xx).^2;           % amplitudes squared -> power spectra
    Cxy  = Xx .* conj(Yy); %= Co+iQ = cross spectrum,
    Pxy = real(Cxy); % Cospectrum
    Qxy = imag(Cxy); % quadrature spectrum
    Phxy = angle(Cxy); % phase spectrum
    % Select first half   
    select = [2:m/2];
    
    % frequency interval
    df   = 1/(m/fs);                            
    % frequency [Hz]
    freq = (select - 1)'./m.*fs;     
    
    % power spectral energy of variable 2
    E2  = 2.*Pyy(select)./m^2./df;
    Ec  = 2.*Pxy(select)./m^2./df;
    Eq  = 2.*Qxy(select)./m^2./df;
    Ph = Phxy(select);
    
    
    res = {[freq E2 Ec Eq Ph]};
    
    
end