wlp
%
function omegahat = wlp(signal,window, SNR) %WLP Find the Weighted Linear Predictor frequency estimates of % the signals in each column of signal. % % omegahat = wlp( signal [,window[, SNR]] ) % % where window is one of: % 'ckq' the Clarkson-Kootsookos-Quinn estimator [1] % 'kay' Kay's estimator [2] <------- (DEFAULT) % 'lovell' the Lovell-Williamson estimator [3] % 'lpr' Lank-Reed-Pollon estimator [4] % and SNR is the ratio of the sinusoid amplitude squared % to the noise variance (i.e. NOT the usual SNR % definition). The SNR is only required for the 'ckq' % option. % % [1] V. Clarkson, P. J. Kootsookos and B. G. Quinn, % ``Variance Analysis of Kay's Weighted Linear Predictor % Frequency Estimator,'' IEEE Transactions on Signal % Processing, Vol. 42, pp. 2370-2379, 1994. % % [2] S. M. Kay, `` A Fast and Accurate Single Frequency % Estimator,'' IEEE Transactions on Acoustics, Speech % and Signal Processing, Vol. 37(12), pp. 1987-1989, % 1989. % % [3] B.C. Lovell and R.C. Williamson, % "The Statistical Performance of Some Instantaneous % Frequency Estimators," IEEE Trans. on Acoustics, % Speech and Signal Processing, Vol. 40, pp. 1708-1723, % 1992. % % [4] G. W. Lank, I. S. Reed and G. E. Pollon, % `` A Semicoherent Detection Statistic and Doppler % Estimation Statistic,'' IEEE Transactions on % Aerospace and Electronic Systems, Vol. AES-9(2), % pp. 151-165, 1973. % % $Id: wlp.m,v 1.4 1999/01/09 11:10:35 PeterK Exp PeterK $ % File: wlp.m % % Copyright (C) 1999 Peter J. Kootsookos % % This software provided under the GNU General Public Licence, which % is available from the website from which this file originated. For % a copy, write to the Free Software Foundation, Inc., 675 Mass Ave, % Cambridge, MA 02139, USA. %Type cast it to double. signal = double(signal); %-------------------------------------------------------- % Error condition checks [T,N] = size(signal); if any([T N]==0) error('wlp: zero size data not allowed.'); end; if T==1 signal = signal(:); [T,N] = size(signal); end; if ~exist('window') window = 'kay'; end; % Need the analytic signal if isreal(signal) signal = conj(hilbert(signal)); end; T = T-1; t = [1:T-1]'; switch (window) case 'ckq' if ~exist('SNR') SNR = 1; end; Theta = log(1 + SNR + sqrt(SNR^2 + SNR)); WN = ( sinh(T*Theta) - sinh(t*Theta) - sinh((T-t)*Theta) ) ... / ( (T-1)*sinh(T*Theta) - 2*sinh(T*Theta/2)*sinh((T-1)/2*Theta) / ... sinh(Theta/2)); case 'kay' WN = (6*t.*(T-t)); case 'lovell' WN = (6*t.*(T-t)); signal = signal./abs(signal); case 'lrp' WN = ones(T-1,1); end; WN = WN*ones(1,N); omegahat = angle(sum([ signal(1:T-1,:).*conj(signal(3:T+1,:)) ].*WN)) / 2; % Author: Peter J. Kootsookos (kootsoop@ieee.org) % % Based on: P.J. Kootsookos, S.J. Searle and B.G. Quinn, % "Frequency Estimation Algorithms," CRC for Robust and % Adaptive Systems Internal Report, June 1993. %