function [y_final f_final ckIter] = mckd(x,filterSize,termIter,T,M,plotMode) %MAXIMUM CORRELATED KURTOSIS DECONVOLTUION % code and method by Geoff McDonald (glmcdona@gmail.com), May 2011 % This code file is an external reference for a paper being submitted % for review. % % mckd(x,filterSize,termIter,plotMode,T,M) % % Description: % This method tries to deconvolve a periodic series of impulses from % a 1d vector. It does this by designing a FIR filter to maximize % a norm criterion called Correlated Kurtosis. This method is has % applications in fault detection of rotating machinery (such as % ball bearing and gear faults). % % Algorithm Reference: % (Paper link coming soon. If you are interested in this, please % contact me at glmcdona@gmail.com. I will add the link if/when the % paper is available online) % % Inputs: % x: % Signal to perform deconvolution on. This should be a 1d vector. % MCKD will be performed on this vector by designing a FIR % filter. % % filterSize: % This is the length of the finite impulse filter filter to % design. Using a value of around 100 is appropriate depending on % the data. Investigate the performance difference using % different values. % % termIter: (OPTIONAL) % This is the termination number of iterations. If the % the number of iterations exceeds this number, the MCKD process % will complete. Specify [] to use default value of 30. % % T: % This is the period for the deconvolution. The algorithm will % try to deconvolve periodic impulses separated by this period. % This period should be specified in number of samples and can be % fractional (such as 106.29). In the case of a fractional T, the % method will resample the data to the nearest larger integer T: % i.e. 106.29 -> 107 % and the y_final output will still be at this resampled factor. % % M: % This is the shift order of the deconvolution algorithm. % Typically an integer value between 1 and 5 is good. Increasing % the number increases the number of periodic impulses it tries % to find in a row. For example M = 5 would try to extract at % least 5 impulses in a row. When you use a larger M you need a % better estimate of T. Using too large a M (approx M > 10) will % result in a loss of numerical precision. % % plotMode: % If this value is > 0, plots will be generated of the iterative % performance and of the resulting signal. % % Outputs: % y_final: % The input signal x filtered by the resulting MCKD filter. % This is obtained simply as: y_final = filter(f_final,1,x); % % f_final: % The final MCKD filter in finite impulse response format. % % ckIter: % Correlated Kurtosis of shift M according to MED iteration. % ckIter(end) is the final ck. % % Example Usage: % % We want to extract the periodic impulses % % from the very strong white noise! % n = 0:999; % x = 3*(mod(n,100)==0) + randn(size(n)); % [y_final f_final ck_iter] = mckd(x,400,30,100,7,1); % M = 7 % % T = 100 % % % Note: % The solution is not guaranteed to be the optimal solution to the % correlated kurtosis maximization problem, the solution is just a % local maximum and therefore a good pick. % Assign default values for inputs if( isempty(filterSize) ) filterSize = 100; end if( isempty(termIter) ) termIter = 30; end if( isempty(plotMode) ) plotMode = 0; end if( isempty(T) ) T = 0; end if( isempty(M) ) M = 1; end % Validate the inputs if( sum( size(x) > 1 ) > 1 ) error('MCKD:InvalidInput', 'Input signal x must be a 1d vector.') elseif( sum(size(T) > 1) ~= 0 || T < 0 ) error('MCKD:InvalidInput', 'Input argument T must be a zero or positive scalar.') elseif( sum(size(termIter) > 1) ~= 0 || mod(termIter, 1) ~= 0 || termIter <= 0 ) error('MCKD:InvalidInput', 'Input argument termIter must be a positive integer scalar.') elseif( sum(size(plotMode) > 1) ~= 0 ) error('MCKD:InvalidInput', 'Input argument plotMode must be a scalar.') elseif( sum(size(filterSize) > 1) ~= 0 || filterSize <= 0 || mod(filterSize, 1) ~= 0 ) error('MCKD:InvalidInput', 'Input argument filterSize must be a positive integer scalar.') elseif( sum(size(M) > 1) ~= 0 || M < 1 || round(M)-M ~= 0 ) error('MCKD:InvalidInput', 'Input argument M must be a positive integer scalar.') elseif( filterSize > length(x) ) error('MCKD:InvalidInput', 'The length of the filter must be less than or equal to the length of the data.') end % Force x into a column vector x = x(:); L = filterSize; OriginalLength = length(x); % Perform a resampling of x to an integer period if required if( abs(round(T) - T) > 0.01 ) % We need to resample x to an integer period T_new = ceil(T); % The rate transformation factor Factor = 20; % Calculate the resample factor P = round(T_new * Factor); Q = round(T * Factor); Common = gcd(P,Q); P = P / Common; Q = Q / Common; % Resample the input x = resample(x, P, Q); T = T_new; else T = round(T); end N = length(x); % Calculate XmT XmT = zeros(L,N,M+1); for( m = 0:M) for( l = 1:L ) if( l == 1 ) XmT(l,(m*T+1):end,m+1) = x(1:N-m*T); else XmT(l, 2:end,m+1) = XmT(l-1, 1:end-1,m+1); end end end % Calculate the matrix inverse section Xinv = inv(XmT(:,:,1)*XmT(:,:,1)'); % Assume initial filter as a delayed impulse f = zeros(L,1); f(round(L/2)) = 1; f(round(L/2)+1) = -1; f_best = f; ck_best = 0; iter_best = 0; % Initialize iteration matrices y = zeros(N,1); b = zeros(L,1); ckIter = []; % Iteratively adjust the filter to minimize entropy n = 1; delta = 0; while n == 1 || ( n <= termIter ) % Compute output signal y = (f'*XmT(:,:,1))'; % Generate yt yt = zeros(N,M); for( m = 0:M ) if( m == 0 ) yt(:,m+1) = y; else yt(T+1:end,m+1) = yt(1:end-T,m); end end % Calculate alpha alpha = zeros(N,M+1); for( m = 0:M ) alpha(:,m+1) = (prod(yt(:,[1:m (m+2):size(yt,2)]),2).^2).*yt(:,m+1); end % Calculate beta beta = prod(yt,2); % Calculate the sum Xalpha term Xalpha = zeros(L,1); for( m = 0:M ) Xalpha = Xalpha + XmT(:,:,m+1)*alpha(:,m+1); end % Calculate the new filter coefficients f = sum(y.^2) / (2*sum(beta.^2)) * Xinv * Xalpha; f = f/sqrt(sum(f.^2)); % Calculate the ck value fo this iteration ckIter(n) = sum(prod(yt,2).^2)/( sum(y.^2)^(M+1) ); % Update the best match filter if required if( ckIter(n) > ck_best ) ck_best = ckIter(n); f_best = f; iter_best = n; end n = n + 1; end % Update the final result f_final = f_best; y_final = filter(f_final,1,x); % Resample the final result if( length(y_final) ~= OriginalLength ) y_final = y_final(1:OriginalLength); end % Plot the results if( plotMode > 0 ) figure; subplot(4,1,1) plot(x) title('Input Signal') ylabel('Value') xlabel('Sample Number') axis tight subplot(4,1,2) plot(y_final) title('Signal Filtered by MCKD') ylabel('Value') xlabel('Sample Number') axis tight subplot(4,1,3) stem(f_final) xlabel('Sample Number') ylabel('Value') title('Final Filter, Finite Impulse Response') axis tight subplot(4,1,4) plot(ckIter); title('CK_m versus algorithm iteration') xlabel('Iteration') ylabel('CK_m(T)') axis tight hold on plot([iter_best], ckIter(iter_best),'-oblack','LineWidth',2) legend('CK_m versus iteration','Best match location') endend
