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/2006Message = 'Ok';  % Set output message to default valueStatusFlag = 0;  % Init trouble flagOutOptionsList = strvcat('raw','norm','dB10','dB20','psd','log10','ln'); % List of valid output options% Start code for error trap on Input dataif 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;	returnelseif 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;	returnelseif 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 trappingif 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	endend% End code for output options error trapping% Make Input a column vectorif size(Input,1) > 1 	Input = Input';	endN = length(Input); % Compute the number of input data pointsMidPoint = fix(N/2); % Computer midpoint of input data% Start code to estimate power spectrumPwrSpecEst = abs(fft(Input)).^2; % Compute raw power spectrum as Real^2 + Imag^2 i.e. raw periodogramCep1 = 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 paperif 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 MuMessage = 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 paperFirstAndLastThreshold = (Mu*pi)/sqrt(3*N);% Compute threshold for the rest of the cepstrumMiddleThreshold = (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 ThresholdVectorThresholdVector(1) = FirstAndLastThreshold;ThresholdVector(N) = FirstAndLastThreshold;% Compute threshold for values of cepstrum from% 2 to N - 1ThresholdVector(2:1:N-1) = MiddleThreshold;% Enter the proper threshold levels for the midpoint % cepstrum value into ThresholdVectorThresholdVector(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 = 1MaskVector = abs(real(Cep1)) > ThresholdVector;% Multiply by mask to set the right coeffs = 0Cep1 = MaskVector.*Cep1;% End code to threshold cepstrum% Start code to compute power spectrum estimate from thresholded cepstrumIndex =[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, etcswitch 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