specwelch(x,dt,win,Nsg,pnv,Nb,P,Wn,ftype,n)

function varargout = specwelch(x,dt,win,Nsg,pnv,Nb,P,Wn,ftype,n)
% Spectrum using Welch's method
%
% USAGE:
%   q = specwelch(x,dt,w,Nsg,pnv,Wn,ftype,n)
%   [psdf,f] = specwelch(x,dt,w,Nsg,pnv,Wn,ftype,n)
%   [psdf,conf,f] = specwelch(x,dt,w,Nsg,pnv,Wn,ftype,n)
%
% DESCRIPTION:
%   Calculates the spectrum for x
%   using Welch's 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'
%   Nsg - Number of Segments (>=1)
%   pnv - Percentage Noverlap of Segments (0-100)
%   Nb - Band Averaging, number of bands to average
%   Wn - Cut-Off frequencies, used for filtering
%   ftype - Type of filter, 'high', 'low' or 'stop'
%   n - 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 || isempty(win)
    win = 'boxcar';    
end
if nargin<4 || isempty(Nsg)
    Nsg = 6;    
end
if nargin<5 || isempty(pnv)
    pnv = 0; %Percent
end
if nargin<6 || isempty(Nb)
   Nb = 1; 
end
if nargin<7 || isempty(P)
    P = 0.95;
end
if nargin<8 && nargin>6
    ftype = 'low';
end
if nargin<9 && nargin>6
    n = 5;
end

%******** IMPORTANT *********
pnv = pnv/100;

if nargin==0
    figure
    specwelchdemo
    return
end

isf = isfinite(x);

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

x = x(isf); %Remove NaN's
x = x(:);
N = sum(isf);

%Calculate number of Segments
% Nsg = (N-nv)./(L-nv); nv = pnv*L;
L = ceil(N/(Nsg - Nsg*pnv + pnv));
nv = floor(pnv*L); %Check- L*Nsg - nv*(Nsg-1)

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

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

q.xp = x; %save before applying filter

%Butterworth Filter
if exist('Wn','var')
    [b1,b2] = butter(n,Wn*2*dt,ftype); %fc/nyquist = fc*2*dt
    x = filtfilt(b1,b2,x);
    hold on
    plot(t,x,'b')
end

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

% 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

%The magnitude is multiplied by 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
index = 1:L;
ind4c = 2:nfft/2+1;
m = 0*f;
s = 0*f;
x = datawrap(x,L*Nsg);
for j=1:Nsg
    %Remove Trend, Apply Window 
    xx = w.*detrend(x(index));
    %Compute fourier coefficients padding with zeros
    %to length nfft, makes fft much faster
    %and increases frequency resolution, df=1/nfft/dt;
    c = fft(xx,nfft); 
    mj = K*abs(c(ind4c)); %skip first coefficient which is mean
    m = m + mj;
    %Sxx(win), Angular-Frequency Power Spectrum, [Power]
    s = s + mj.^2;
    index = index + L - nv; %new index
end

% Average the sum of the periodograms
m = m/Nsg; %Magnitude
s = 4*s/Nsg/R; %Includes correction for one-sided spectrum

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

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

%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; 

%Confidence Intervals using chi-squared approach
v = 2*Nsg*Nb; %degrees of freedom
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

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(f,psdf)
    hold on
    loglog(f,psdf*conf,'--r')
    axis tight
    set(gca,'xticklabel',get(gca,'xtick'))
    xlabel('Frequency')
    ylabel('PSD(f)')

    %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
        w = eval([win,'(',num2str(N),')']);
        plot(w)
        set(gca,'xticklabel',get(gca,'xtick'))
        xlabel('Period')
        ylabel('window')
        axis tight
    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 specwelchdemo
dt = 1/100;
t=0:dt:200;
sr = 2*randn(size(t)); %white noise
s1 = 4*sin(2*pi*t/1 + pi/2);
s2 = 2*sin(2*pi*t/4);
s3 = 1*sin(t*2*pi/10);
s4 = 3*sin(t*2*pi/0.5);
s5 = -0.05*t;
s0 = 5;
x = s0 + s1 + s2 + s3 + s4 + s5 + sr;
x(1000:6000) = NaN;
w = 'hanning';
Wn = 2;
ftype = 'low';
n = 3;
Nsg = 10;
pnv = 30;
Nb = 1;
specwelch(x,dt,w,Nsg,pnv,Nb);

return