cepspec.m

无参数功率谱估计

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 [e-mail] 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