www.gusucode.com > signal 工具箱matlab源码程序 > signal/cceps.m
function [xhat,nd,xhat1] = cceps(x,n) %CCEPS Complex cepstrum. % CCEPS(X) returns the complex cepstrum of the real sequence X using % Fourier transform as in its definition [1]. The input is altered to % have no phase discontinuity at +/- pi radians by the application of a % linear phase term (that is, it is circularly shifted (after % zero-padding) by some samples if necessary to have zero phase at pi % radians). The program follows the procedure used in [2] % % [XHAT,ND] = CCEPS(X) returns the number of samples ND of (circular) % delay added to X prior to finding the complex cepstrum. % % [XHAT,ND,XHAT1] = CCEPS(X) returns the complex cepstrum in XHAT1 using % an alternate factorization algorithm developed in [3]. The first % element in the data sequence cannot be zero. % % CCEPS(X,N) zero-pads X to length N, and returns the length N complex % cepstrum of X. % % The comparison between the two algorithms can be seen from the % following table: % % ---------------------------------------------------------------------- % | Algorithm | Pros | Cons | % ---------------------------------------------------------------------- % | Fourier | Can be applied to | Needs phase unwrapping. | % | | any signal. | Output is aliased. | % ---------------------------------------------------------------------- % |Factorization | Does not need phase | The input signal has to have | % | | unwrapping. | an all-zero Z transform. | % | | No aliasing. | Z transform should have no | % | | | zeros on the unit circle. | % ---------------------------------------------------------------------- % % Generally, the results of the two algorithms cannot be used to verify % each other. % % The two results can be used to verify each other only when the first % element of data sequence is positive, the Z transform of the data % sequence has only zeros, all these zeros are inside the unit circle, % and the data sequence is long or padded with enough zeros. % % EXAMPLE: Use CCEPS to show an echo. % Fs = 100; t = 0:1/Fs:1.27; % s1 = sin(2*pi*45*t); % 45Hz sine wave sampled at 100Hz % % % Add an echo of the signal, with half the amplitude, 0.2 seconds after % % the beginning of the signal. % s2 = s1 + 0.5*[zeros(1,20) s1(1:108)]; % % c = cceps(s2); % plot(t,c) % title('The peak at 0.2 shows the echo'); % % See also ICCEPS, RCEPS, HILBERT, and FFT. % Copyright 1988-2013 The MathWorks, Inc. % References: % [1] Oppenheim, A.V. and Schafer, R.W. Discrete-Time Signal % Processing, Prentice-Hall, 1989. % [2] Programs for Digital Signal Processing, IEEE Press, % John Wiley & Sons, 1979, algorithm 7.1. % [3] Steiglitz, K. and Dickinson, B. "Computation of the complex % cepstrum by factorization of the Z-transform," Proc. Int. Conf. % ASSP, 1977, 723-726 narginchk(1,2); % Check for valid input signal chkinput(x); if nargin < 2 h = fft(x); else n = signal.internal.sigcasttofloat(n,'double','cceps','N',... 'allownumeric'); h = fft(x,n); end [ah,nd] = rcunwrap(angle(h)); logh = log(abs(h))+1i*ah; xhat = real(ifft(logh)); if nargout==3 % use alternate rooting algorithm to check original result % Oppenheim & Schafer, p.795 [1] and Steiglitz & Dickinson [3] % - doesn't work with zeros on the unit circle % - only works for finite sequences that can be rooted x = x(:).'; % Force vector to be a row r = roots(x); r = r(r~=0); if ~isempty(r(abs(r) == 1)) error(message('signal:cceps:SignalErr')); end a = r(abs(r)<1); b = r(abs(r)>1); A = x(1)*prod(b); b = 1./b; if nargin == 2 datalen = n; else datalen = length(x); end % class(xhat1) should be same as class(xhat) to enforce precision rules xhat1 = zeros(datalen,1,class(xhat)); %#ok<*ZEROLIKE> xhat1(1) = log(abs(A)); % n > 0 contributions are from zeros inside the unit circle if ~isempty(a) n = 1:datalen-1; [a,n]=meshgrid(a,n'); xhat1 = xhat1 + [0; -((a.^n)./n)*ones(size(a,2),1)]; end % n < 0 contributions are from zeros outside the unit circle if ~isempty(b) n = -(datalen-1):-1; [b,n]=meshgrid(b,n'); xhat1 = xhat1 + [0; ((b.^(-n))./n)*ones(size(b,2),1)]; end if size(x,1)==1 xhat1 = xhat1'; end end %-------------------------------------------------------------------------- function [y,nd] = rcunwrap(x) %RCUNWRAP Phase unwrap utility used by CCEPS. % RCUNWRAP(X) unwraps the phase and removes phase corresponding % to integer lag. See also: UNWRAP, CCEPS. n = length(x); y = unwrap(x); nh = fix((n+1)/2); idx = nh+1; % Special case the index for scalar input. if length(y) == 1 idx = 1; end nd = round(y(idx)/pi); y(:) = y(:)' - pi*nd*(0:(n-1))/nh; % Cast to 'double' to enforce precision rules nd = double(nd); % output nd has no bearing on single precision rules %-------------------------------------------------------------------------- function chkinput(x) % Check for valid input signal signal.internal.sigcheckfloattype(x,'','cceps','X'); if isempty(x) || issparse(x) || ~isreal(x), error(message('signal:cceps:invalidInput', 'X')); end