function [SpectrumOut, Message, StatusFlag] = CepSpec(Input, BandInfo, OutOptions)
% Function to compute the power spectrum from the cepstrum via the method
% described in "Smoothed Nonparametric Spectral Estimation Via Cepstrum
% Thresholding" by Perte Stocia and Niclas Sandgren, IEEE Signal Processing
% Magazine, November 2006 vol 23, number 6
%
% This software copyright by Steve Lutes 11/11/2006. All rights reserved
% Commerical use requires a license. Non-commerical use requires acknowledging
% Steve Lutes as author of this software. Send e-mail to utahraptor@intergate.com for
% more info.
%
% How to use:
%
% [SpectrumOut, Message, StatusFlag] = CepSpec(Input, BandInfo, OutOptions)
%
%---Outputs------------------------------------------------------------------------------
%
% SpectrumOut- an array length(Input) / 2 points long holding the
% power spectrum of Input.
%
% Message - a text string telling of problems, if any. If no trouble, Message = 'Ok' plus
% whatever output option user has specified
%
% StatusFlag- a numeric scalar telling of trouble. = 0 means all ok, = 1 means trouble
%
%---Inputs-------------------------------------------------------------------------------
%
% Input- a numeric array (row or colum vector) of input data
%
% BandInfo- a string of characters. Enter 'narrow' for a narrow-band signal
% or 'wide' for a wide-band signal. Mixed, all upper - or all lower-case
% will work. e.g. 'narrow', 'NARROW', 'Wide'.
%
%
% OutOptions- a string determining how the estimate of the power specta is
% numerically altered. Mixed case, all lower case or all upper case are
% valid. Options are:
%
% [] or 'raw' - nothing happens. Spectral estimate is not altered
%
% 'norm' - All spectral values are normalized to one. This is done
% by dividing all values by the maximum. This
% PowerSpectrumNormalized(i) = PowerSpectrum(i) / max(PowerSpectrum)
%
%
% 'dB10' - PowerSpectrumdB10(i) = 10*log10[PowerSpectrum(i) / max(PowerSpectrum)]
% This is for a power signal
%
% 'dB20' - PowerSpectrumdB20(i) = 20*log10[PowerSpectrum(i) / max(PowerSpectrum)]
% This is for a voltage or current signal
%
% 'psd' - PowerSpectrum_psd(i) = PowerSpectrum(i)/ sum(PowerSpectrum(i))
% This gives the power spectral density function. Dividing out all
% values of PowerSpectrum(i) by sum(PowerSpectrum(i)) makes sure that
% sum(PowerSpectrum_psd(i)) = 1. This is how a probability density function
% behaves.
%
% 'log10' - PowerSpectrum_log10(i) = log10(PowerSpectrum(i)). Log10 is the
% log to the base 10
%
% 'ln' - PowerSpectrum_ln(i) = ln(PowerSpectrum(i)) or log to the base e (natural log)
%
% All numerics ok 11/11/2006
Message = 'Ok'; % Set output message to default value
StatusFlag = 0; % Init trouble flag
OutOptionsList = strvcat('raw','norm','dB10','dB20','psd','log10','ln'); % List of valid output options
% Start code for error trap on Input data
if isempty(Input) == 1
% Trap for empty value of input data. Ok 11/12/2006
disp(['Empty value for input data entered. Please use a row or column vector for the input data.'])
Message = ['Empty value for input data entered. Please use a row or column vector for the input data.'];
StatusFlag = 1;
return
elseif isnumeric(Input) == 0
% Trap for non-numeric value(s) of input data. Ok 11/12/2006
disp(['Non-numeric value(s) for input data entered. Please use a numeric row or column vector for the input data.'])
Message = ['Non-numeric value(s) for input data entered. Please use a numeric row or column vector for the input data.'];
StatusFlag = 1;
return
elseif sum(size(Input)) <= 2
% Trap for scalar value(s) of input data. Ok 11/12/2006
disp(['Scalar for input data entered. Please use a numeric row or column vector for the input data.'])
Message = ['Scalar for input data entered. Please use a numeric row or column vector for the input data.'];
StatusFlag = 1;
return
elseif all(size(Input) > 1) == 1
% Trap for scalar value(s) of input data. Ok 11/12/2006
disp(['Matrix for input data entered. Please use a numeric row or column vector for the input data.'])
Message = ['Matrix for input data entered. Please use a numeric row or column vector for the input data.'];
StatusFlag = 1;
return
end
% End code for error trap on Input data
% Start code to error-trap the bandwidth type specification
if isempty(BandInfo) == 1
% Trap for empty value of BandInfo Ok 11/12/2006
disp(['Empty value for bandwidth spec (wide vs. narrow) entered. Please enter ''wide'' or ''narrow''. Mixed, lower or upper case is ok.'])
Message = ['Empty value for bandwidth spec (wide vs. narrow) entered. Please enter ''wide'' or ''narrow''. Mixed, lower or upper case is ok.'];
StatusFlag = 1;
return
elseif ischar(BandInfo) == 0
% Trap for non-character value of BandInfo Ok 11/12/2006
disp(['Non character value for bandwidth spec (wide vs. narrow) entered. Please enter ''wide'' or ''narrow'' *as a text string*. Mixed, lower or upper case is ok.'])
Message = ['Non character value for bandwidth spec (wide vs. narrow) entered. Please enter ''wide'' or ''narrow'' *as a text string*. Mixed, lower or upper case is ok.'];
StatusFlag = 1;
return
else
BandInfo = lower(BandInfo);
if isempty(strmatch(BandInfo, strvcat('narrow','wide') , 'exact')) == 1
% Trap for non-character value of BandInfo Ok 11/12/2006
disp(['Invalid option for bandwidth spec (wide vs. narrow) entered. Please enter ''wide'' or ''narrow'' *as a text string*. Mixed, lower or upper case is ok.'])
Message = ['Invalid option for bandwidth spec (wide vs. narrow) entered. Please enter ''wide'' or ''narrow'' *as a text string*. Mixed, lower or upper case is ok.'];
StatusFlag = 1;
return
end
end
% End code to error-trap the bandwidth type specification
% Start code for output options error trapping
if isempty(OutOptions) == 1;
OutOptions = ['raw'];
elseif ischar(OutOptions) == 0;
% Trap for non-character value of OutOptions Ok 11/12/2006
disp(['Non-character value entered for OutOptions. Please enter ''raw'', ''norm'', ''dB10'', ''dB20'', ''psd'', ''log10'', or ''ln''. Mixed, lower or upper case is ok. ']);
Message = ['Non-character value entered for OutOptions. Please enter ''raw'', ''norm'', ''dB10'', ''dB20'', ''psd'', ''log10'', or ''ln''. Mixed, lower or upper case is ok. '];
StatusFlag = 1;
return
else
OutOptions = lower(OutOptions);
if isempty(strmatch(OutOptions, OutOptionsList , 'exact')) == 1
% Trap for invalid value of OutOptions
disp(['Invalid output option. Please enter ''raw'', ''norm'', ''dB10'', ''dB20'', ''psd'', ''log10'', or ''ln''. Mixed, lower or upper case is ok. ']);
Message = ['Invalid output option. Please enter ''raw'', ''norm'', ''dB10'', ''dB20'', ''psd'', ''log10'', or ''ln''. Mixed, lower or upper case is ok. '];
StatusFlag = 1;
return
end
end
% End code for output options error trapping
% Make Input a column vector
if size(Input,1) > 1
Input = Input';
end
N = length(Input); % Compute the number of input data points
MidPoint = fix(N/2); % Computer midpoint of input data
% Start code to estimate power spectrum
PwrSpecEst = abs(fft(Input)).^2; % Compute raw power spectrum as Real^2 + Imag^2 i.e. raw periodogram
Cep1 = ifft(log(PwrSpecEst)); % Compute cepstrum
% Start code to threshold cepstrum via equation 36 in Stocia and Sandgren paper
% Start code to determine value of Mu via
% equations 37 and 38 in Stocia and Sandgren paper
if strcmp(BandInfo,'wide') == 1
% Start code for the wideband case
if N < 500
% Case where signal < 500 points
Mu = 4;
else
% Case where signal = or > 500 points
Mu = 5;
end
% End code for the wideband case
elseif strcmp(BandInfo,'narrow') == 1
% Start code for the narrow band case
if N < 500
% Case where signal < 500 points
Mu = 2;
else
% Case where signal = or > 500 points
Mu = 3;
end
% End code for the narrow band case
end
% End code to determine value of Mu
Message = strvcat(Message,['Signal is designated as ' BandInfo '-band.'], ['Mu = ' num2str(Mu)], ['Number of Data Points = ' num2str(N)]);
% Start code for thresholding cepstrum
% Compute thresholds for first, last and (N/2) element of
% cepstrum. Mu is from eqns 37 and 38 in Stocia and Sandgren paper
% Equation for threshold is equation 36 in Stocia and Sandgren paper
FirstAndLastThreshold = (Mu*pi)/sqrt(3*N);
% Compute threshold for the rest of the cepstrum
MiddleThreshold = (Mu*pi)/sqrt(6*N);
% Form vector of threshold values for threshold comparisons
% Enter the proper threshold levels for the first,
% and last cepstrum values into ThresholdVector
ThresholdVector(1) = FirstAndLastThreshold;
ThresholdVector(N) = FirstAndLastThreshold;
% Compute threshold for values of cepstrum from
% 2 to N - 1
ThresholdVector(2:1:N-1) = MiddleThreshold;
% Enter the proper threshold levels for the midpoint
% cepstrum value into ThresholdVector
ThresholdVector(MidPoint) = FirstAndLastThreshold;
% Compute a masking vector. All places where
% abs(real(Cep1)) < or = threshold will be zero
% in the masking vector. If abs(real(Cep1)) > Threshold
% then the corresponding place in masking vector = 1
MaskVector = abs(real(Cep1)) > ThresholdVector;
% Multiply by mask to set the right coeffs = 0
Cep1 = MaskVector.*Cep1;
% End code to threshold cepstrum
% Start code to compute power spectrum estimate from thresholded cepstrum
Index =[1:1:MidPoint]; % Define a index vector [1.....N2]
PwrSpecEstCep = exp(real(fft(Cep1))); % Power spectrum estimate from thresholded cepstrum
% Compute scale factor (equation 41 in Stocia and Sandgren paper)
ScaleConstant = sum(PwrSpecEst(Index).*PwrSpecEstCep(Index))/ sum(PwrSpecEstCep(Index).^2);
SpectrumOut = ScaleConstant*PwrSpecEstCep(Index); % Scale to get final result
% End code to compute power spectrum estimate from thresholded cepstrum
% End code to estimate power spectrum
% Start code to make the output the PSD, DB, normalized, etc
switch OutOptions;
case {'raw'}
% Do nothing to power spectrum
Message = strvcat(Message, 'Raw output option selected.', 'Power spectrum has not been altered');
case {'norm'}
% User wants output spectra normalized by highest value
% Normalize power spectrum by dividing all values by max(power spectrum)
SpectrumOut = SpectrumOut./max(SpectrumOut);
% Give a message
Message = strvcat(Message, 'Normalized output option selected.', 'Power spectrum has normalized by its maximum.');
case {'db10'}
% User wants output spectra in power dB
% i.e. PowerSpectrumdB10(i) = 10*log10[PowerSpectrum(i) / max(PowerSpectrum)
% Compute power spectrum for the dB10 case i.e.
% i.e. PowerSpectrumdB10(i) = 10*log10[PowerSpectrum(i) / max(PowerSpectrum)
SpectrumOut = 10*log10(SpectrumOut./max(SpectrumOut));
% Give a message
Message = strvcat(Message, 'Decibel (power ratio i.e. 10*log10[PowerSpectrum(i) / max(PowerSpectrum) ) output option selected.', 'Power spectrum is in decibels.');
case {'db20'}
% User wants output spectra in voltage or current dB
% i.e. PowerSpectrumdB20(i) = 20*log10[PowerSpectrum(i) / max(PowerSpectrum)
% Compute power spectrum for the dB20 case i.e.
% i.e. PowerSpectrumdB10(i) = 20*log10[PowerSpectrum(i) / max(PowerSpectrum)
SpectrumOut = 20*log10(SpectrumOut./max(SpectrumOut));
% Give a message
Message = strvcat(Message, 'Decibel (voltage or current ratio i.e. 20*log10[PowerSpectrum(i) / max(PowerSpectrum) ) output option selected.', 'Power spectrum is in decibels.');
case {'psd'}
% User wants output as probablity density function i.e the power spectral density function
SpectrumOut = SpectrumOut./sum(SpectrumOut);
% Give a message
Message = strvcat(Message, 'Power spectral density option selected', 'Power spectrum is a probability density function.');
case{'log10'}
% User wants output as PowerSpectrum_log10(i) = log10(PowerSpectrum(i))
SpectrumOut = log10(SpectrumOut);
% Give a message
Message = strvcat(Message, 'Log to base 10 option selected (PowerSpectrumOut(i) = log10(PowerSpectrum(i))', 'Power spectrum is in log to the base 10 units.');
case{'ln'}
% User wants output as PowerSpectrum_n(i) = ln(PowerSpectrum(i))
SpectrumOut = log(SpectrumOut);
% Give a message
Message = strvcat(Message, 'Log to base e (natural log) option selected (PowerSpectrumOut(i) = ln(PowerSpectrum(i))', 'Power spectrum is in log to the base e (natural log) units.');
end
% End code to make the output the PSD, DB, normalized, etc
