function [f,k,p_S,p_A,p_T] = wk99(data,T);

%
%  [f,k,p_S,p_A,p_T,p_bg] = wk99(data,T)
%  
%  Frequency-wavenumber diagram as described in Wheeler and Kiladis (1999).
%  Note that this script DOES NOT calculate the background spectrum and
%  calculate the ratio of the symmetric and anti-symmetric parts to the
%  background. This can be done using code something like:
%   > sum(p_S,3) ./ two_dim_smoothing(sum(p_S,3))
%  to give a symmetric spectrum "anomaly" as presented in WK99
%
%  The input quantities are as follows:
%
%   data - data from which to generate the WK99 diagram. Should have
%          dimensions of [ time x longitude x latitude ]. The latitudes
%          must be symmetric about the equator
%   T    - length of segments from which to generate independent spectra
%          (e.g., 96 in the WK99 paper). The final spectra are then generated
%          by averaging over all the segments.
%
%  The output quantities are as follows:
%  
%   f    - frequencies [ rad / "sample rate" ]
%   k    - wavenumbers [ rad / "sample rate" ]
%   p_S  - 2D spectral density of symmetric variability
%   p_A  - 2D spectral density of anti-symmetric variability
%   p_T  - 2D spectral density of total variability
%  

% some sizes
X = size(data,2); Y = size(data,3);

% sumber of segments to average over
nseg = floor(size(data,1)/T); if ( (nseg*T)+T/2 > size(data,1) ) nseg = nseg-1; end;

% generate tapering window
taperwin = repmat(tukeywin(T,0.25),1,2*nseg);

% allocate matrices
p_T = zeros(X,T,Y); p_S = zeros(X,T,Y); p_A = zeros(X,T,Y);

% loop over latitudes
for j = 1:Y
  % get [time x latitude] matrix of var (var_eq)
  var_eq = data(:,:,j); % var at lat
  var_eq2 = data(:,:,Y-j+1); % var at -lat
  var_eq_S = 0.5*( var_eq + var_eq2 ); % symmetric component
  var_eq_A = 0.5*( var_eq - var_eq2 ); % anti-symmetric component
  var_eq_S = detrend(var_eq_S); var_eq_A = detrend(var_eq_A); % detrend temporally
  var_eq_S = detrend(var_eq_S')'; var_eq_A = detrend(var_eq_A')'; % detrend spatially
  % -----
  % TOTAL
  % segments
  for i = 1:size(var_eq,2)
    var_seg_tmp = reshape( var_eq(1:(nseg*T),i) , T , nseg ); % 96 day segments
    var_seg_tmp = [var_seg_tmp reshape( var_eq((T/2)+[1:(nseg*T)],i) , T , nseg )]; % overlapping segments (96/2 day overlap)
    var_seg_tmp = detrend(var_seg_tmp).*taperwin; % taper t.s.
    var_seg(:,i,:) = reshape( var_seg_tmp , [size(var_seg_tmp,1) 1 size(var_seg_tmp,2)]);
  end
  % calculate average w-k spectra over all segments
  p_seg = zeros(X,T);%length(floor(T/2)+2:T),X);
  for seg = 1:size(var_seg,3)
    p0 = abs(fftshift(circshift(fft2(var_seg(:,:,seg)'),[0 -1])));% fliplr required to make (T/2+1,X/2+1) element at (w,k)=(0,0)
    p_seg = p_seg + fliplr(p0);%(floor(T/2)+2:end,:));
  end
  p_T(:,:,j) = p_seg./size(var_seg,3);
  % ---------
  % SYMMETRIC
  % segments
  for i = 1:size(var_eq_S,2)
    var_seg_tmp = reshape( var_eq_S(1:(nseg*T),i) , T , nseg ); % T-day segments
    var_seg_tmp = [var_seg_tmp reshape( var_eq_S((T/2)+[1:(nseg*T)],i) , T , nseg )]; % overlapping segments (T/2 day overlap)
    var_seg_tmp = detrend(var_seg_tmp).*taperwin; % taper t.s.
    var_seg(:,i,:) = reshape( var_seg_tmp , [size(var_seg_tmp,1) 1 size(var_seg_tmp,2)]);
  end
  % calculate average w-k spectra over all segments
  p_seg = zeros(X,T);%length(floor(T/2)+2:T),X);
  for seg = 1:size(var_seg,3)
    p0 = abs(fftshift(circshift(fft2(var_seg(:,:,seg)'),[0 -1])));% fliplr required to make (T/2+1,X/2+1) element at (w,k)=(0,0)
    p_seg = p_seg + fliplr(p0);%(floor(T/2)+2:end,:));
  end
  p_S(:,:,j) = p_seg./size(var_seg,3);
  % --------------
  % ANTI-SYMMETRIC
  % segments
  for i = 1:size(var_eq_A,2)
    var_seg_tmp = reshape( var_eq_A(1:(nseg*T),i) , T , nseg ); % 96 day segments
    var_seg_tmp = [var_seg_tmp reshape( var_eq_A((T/2)+[1:(nseg*T)],i) , T , nseg )]; % overlapping segments (96/2 day overlap)
    var_seg_tmp = detrend(var_seg_tmp).*taperwin; % taper t.s.
    var_seg(:,i,:) = reshape( var_seg_tmp , [size(var_seg_tmp,1) 1 size(var_seg_tmp,2)]);
  end
  % calculate average w-k spectra over all segments
  p_seg = zeros(X,T);%length(floor(T/2)+2:T),X);
  for seg = 1:size(var_seg,3)
    p0 = abs(fftshift(circshift(fft2(var_seg(:,:,seg)'),[0 -1])));% fliplr required to make (T/2+1,X/2+1) element at (w,k)=(0,0)
    p_seg = p_seg + fliplr(p0);%(floor(T/2)+2:end,:));
  end
  p_A(:,:,j) = p_seg./size(var_seg,3);
end

% create f and k vectors (units: rad/"cell no")
f = 2*pi*[0:T-1]/T; f(find(f>=pi)) = f(find(f>=pi))-2*pi; f = fftshift(f);
k = 2*pi*[0:X-1]/X; k(find(k>=pi)) = k(find(k>=pi))-2*pi; k = fftshift(k);
