spectral(x,dt,win,Nb,P,Wn,ftype,n)

function varargout = spectral(x,dt,win,Nb,P,Wn,ftype,n)
% Spectrum using periodogram method
%
% USAGE:
%   q = spectral(x,dt,win,Wn,ftype,n)
%   [psdf,f] = spectral(x,dt,win,Wn,ftype,n)
%   [psdf,psdfc,f] = spectral(x,dt,win,Wn,ftype,n)
%
% DESCRIPTION:
%   Calculates the spectrum for x
%   using the periodogram method 
%   If a window other than boxcar is used
%   than the method is refered to as modified
%   periodogram method.
%   The confidence intervals are calculated
%   using the inverse of chi-square CDF.
%   Also includes a filtering option using the 
%   butterworth filter to see the effect of the 
%   filter on the spectrum
%
% INPUT VARIABLES:
%   x - Time series, [vector]
%   dt - Sampling Rate, [scalar]
%   win - Window, one of: 
%       'hanning', 'hamming', 'boxcar'
%   Nb - Band Averaging, number of bands to average
%   P - Probability for confidence intervals
%   Wn - Cut-Off frequencies, used for filtering
%   ftype - Type of filter, 'high', 'low' or 'stop'
%   ncb - Number of coefficients to use in
%       the Butterworth filter
%
% OUTPUT VARIABLES:
%   q - structure with the following fields:
%       xp - detrended x
%       f = Frequencies
%       T - Periods
%       m - Magnitude
%       a - Amplitude
%       s - Power spectrum, Sxx(win), [Power]
%       psdw - Power Spectral Density, Pxx(win), [Power/rad/sample]
%       psdf - Power Spectral Density, Pxx(f), [Power/sample-freq]
%       psdT - Power Spectral Density, Pxx(T), [Power*time-unit]
%       conf - Upper and Lower Confidence Interval multiplication 
%               factors using chi-squared approach 
%
%Copy-Left, Alejandro Sanchez

if nargin<2
    dt = 1;    
end
if nargin<3
    win = 'boxcar';    
end
if nargin<4 || isempty(Nb)
    Nb = 1; 
end
if nargin<5 || isempty(P)
    P = 0.95;
end
if (nargin<7) && (nargin>5)
    ftype = 'low';
end
if (nargin<8) && (nargin>5)
    n = 5;
end

if nargin==0
    spectraldemo
    return
end

x = x(:);
isf = isfinite(x);

if nargout==0
    figure
    t = dt*(0:(length(x)-1));
    subplot(4,1,1)
    plot(t,x,'b')
    hold on
    t = t(isf);
end

x = x(isf);
mx = mean(x);
x = detrend(x); %also removes mean
N = sum(isf);

fn = (1/(2*dt)); %Nyquist frequency
nfft = max(256, 2^(nextpow2(N)));
f = (1:nfft/2)/(nfft/2)*fn; %frequencies, Note: starts from 1
f = f(:); %Frequencies

%Makes a window
try
    w = eval([win,'(',num2str(N),')']);
catch
    error('Invalid Window')
end

%Butterworth Filter
if nargin>5
    [b1,b2] = butter(n,Wn*2*dt,ftype); %Wn/nyquist = Wn*2*dt
    x = filtfilt(b1,b2,x);
    xf = x;
end

if nargout==0
    plot(t,x+mx,'r')
    axis tight
    set(gca,'xlim',[t(1),t(end)])
end

q.xp = x; %save before applying window

x = w.*x; %apply window

% Evaluate the window normalization constant
% A 1/nfft factor has been omitted since it will cancel below
R = w'*w;  % compute sum of squares of window

%Compute Fourier Coefficients
%Pad with zeros to length nfft, makes fft much faster
%and increases frequency resolution, df=1/N/dt;
c = fft(x,nfft); 

%The magnitude is multiplied by one of the following constants
%to compensate for the lost energy when the window is applied
%the values are taken from Emery & Thomson pag. 446-448
%In the case of the boxcar window the constant is 1 and is
%therefore ommitted
if strcmpi(win,'hanning')
    K = 1/2*sqrt(8/3);
elseif strcmpi(win,'hamming')
    K = 1/2*sqrt(5/2);
elseif strcmpi(win,'boxcar')
    K = 1;
else
    error('Invalid Window')
end

%Magnitude
m = K*abs(c(2:nfft/2+1)); %skip first coefficient which is the mean

%Do band averaging now
if Nb>1
    f = bandaverage(f,Nb);
    m = bandaverage(m,Nb);
end
T = 1./f; %Periods

%Amplitude
a = m.*2./R;

%Sxx(win), Angular-Frequency Power Spectrum, [Power]
s = 2*(m.^2)./R; 

%Correction for one-sided spectrum
s(2:end-1) = 2*s(2:end-1);

%Pxx(win), Amgular-Frequency Spectrum
%Power Spectrum Density (PSD), [Power/rad/sample]
psdw = s/2/pi;

%Pxx(f), Frequency Spectrum
%Power Spectrum Density (PSD), [Power/sample-freq]
psdf = s*dt;

%Pxx(T), Period Spectrum
%Power Spectrum Density (PSD), [Power*time-unit]
psdT = 2*pi*psdw./T.^2; 

%Finally Calculate phase
ph = unwrap(angle(c(1:fix(N/2))));

%Confidence Intervals using chi-squared approach
v = 2*Nb;
alpha2 = (1-P)/2;
conf = v./chi2inv([alpha2,1-alpha2],v);

%Collect variables
switch nargout
    case 0
        %do nothing
    case 1
        varargout{1}.f = f;
        varargout{1}.T = T;
        varargout{1}.m = m;
        varargout{1}.a = a;
        varargout{1}.s = s;
        varargout{1}.psdw = psdw;
        varargout{1}.psdf = psdf;
        varargout{1}.psdT = psdT;
        varargout{1}.conf = conf;
    case 2
        varargout{1} = psdf;
        varargout{2} = f;
    case 3
        varargout{1} = psdf;
        varargout{2} = conf;
        varargout{3} = f;
    otherwise
        error('Invalid number of output variables')
end %switch

%Plotting
if nargout==0
    subplot(4,1,2)
    semilogx(T,a)
    set(gca,'xlim',[T(end),T(1)])
    set(gca,'xticklabel',get(gca,'xtick'))
    xlabel('Period')
    ylabel('Amplitude')
    
    subplot(4,1,3)
    loglog(T,psdT)
    hold on
    loglog(T,psdT*conf,'--r')
    axis tight
    set(gca,'xticklabel',get(gca,'xtick'))
    xlabel('Period')
    ylabel('PSD(T)')
   
    %Filter or window
    subplot(4,1,4)
    if exist('Wn','var')
        H = freqz(b1,b2,f*pi*2*dt);
        semilogx(T,abs(H),'b');
        axis tight
        set(gca,'xticklabel',get(gca,'xtick'))
        xlabel('Period')
        ylabel('Filter Response')
    else
        plot(w)
        ylabel('window')
        set(gca,'xticklabel',get(gca,'xtick'))
        axis tight
%         plot(ph)
%         xlabel('Period')
    end
end

return
%----------------------------------------------------
function yy = bandaverage(y,Nb)
%Does the band averaging
Lb = floor(length(y)/Nb);
yy = reshape(y(1:Nb*Lb),Nb,Lb); 
yy = mean(yy,1); 
y = y(:);
return

%-------------------------------------------
function spectraldemo
dt = 1/100;
t=0:dt:200;
sr = 2*randn(size(t)); %white noise
s1 = 4*sin(2*pi*t/1);
s2 = 2*sin(2*pi*t/4);
s3 = 1*sin(t*2*pi/10);
s4 = 3*sin(t*2*pi/20);
s5 = -0.05*t;
s0 = 5;
x = s0 + s1 + s2 + s3 + s4 + s5 + sr;
x(1000:6000) = NaN;
Wn = 1;
ftype = 'low';
win = 'hanning';
Nb = [];
P = 0.95;
n = 5;
spectral(x,dt,win,Nb,P,Wn,ftype,n)

return