多点最优最小熵解卷积算法[附源码]

% Simple vibration fault model    close all    n = 0:4999;    h = [-0.05 0.1 -0.4 -0.8 1 -0.8 -0.4 0.1 -0.05];    faultn = 0.05*(mod(n,50)==0);    fault = filter(h,1,faultn);    noise = wgn(1,length(n),-25);    x = sin(2*pi*n/30) + 0.2*sin(2*pi*n/60) + 0.1*sin(2*pi*n/15) + fault;    xn = x + noise;
    % No window. A 5-sample recangular window would be ones(5,1).    window = ones(1,1);
    % 1000-sample FIR filters will be designed    L = 1000;
    % Recover the fault signal of period 50    [MKurt f y] = momeda(xn,L,window,50,1);

function [MKurt f y] = momeda(x,filterSize,window,period,plotMode)    % MULTIPOINT OPTIMAL MINUMUM ENTROPY DECONVOLUTION ADJUSTED    %       code by Geoff McDonald (glmcdona@gmail.com), 2015    %       Used in my research with reference to paper:    %    % momeda(x,filterSize,window,period,plotMode)    %  Multipoint Optimal Minimum Entropy Deconvolution (MOMEDA) computation algorithm. This proposed    %  method solves the optmial solution for deconvolving a periodic train of impulses from a signal.    %  It is best-suited in application to rotating machine faults from vibration signals, to deconvolve    %  the impulse-like vibration associated with many gear and bear faults.    %    %  This method is derived in the Algorithm Reference section.    %    %    % Algorithm Reference:    %    Preparing to publish.    %    G.L. McDonald, <others>, Multipoint Optimal Minimum Entropy Deconvolution and Convolution    %    Fix: Application to Vibration Fault Detection, unpublished    %    % Inputs:    %    x:     %       Signal to generate apply MOMEDA on. Generally this should be around the range    %       of 1000 to 10,000 samples covering at least 5 rotations of the elements in the machine.    %     %    filterSize:    %       This is the length of the finite inpulse filter filter to     %       design. Generally a number on the order of 500 or 1000 is good, but may    %       depend on the dataset length.    %     %    window:    %       This is the window that be convolved with the impulse train target. Generally, a    %       rectangular window works well, eg [1 1 1 1 1].    %     %    period:    %       This is the periods to test as the spectrum x-axis. It should be a decimal range, like:    %           range = 5:0.1:300;    %     %    plotMode:    %       If this value is > 0, plots will be generated of the iterative    %       performance and of the resulting signal.    %    % Outputs:    %    MKurt:    %       The Multipoint Kurtosis of the Deconvolution result.    %    %    f:    %       Optimal FIR filter designed.    %    %    y:    %       Filtered output signal.    %           % Assign default values for inputs    if( isempty(filterSize) )        filterSize = 300;    end    if( isempty(plotMode) )        plotMode = 0;    end    if( isempty(window) )        window = ones(1,1);    end        if( sum( size(x) > 1 ) > 1 )        error('MOMEDA:InvalidInput', 'Input signal x must be 1d.')    elseif(  sum(size(plotMode) > 1) ~= 0 )        error('MOMEDA:InvalidInput', 'Input argument plotMode must be a scalar.')    elseif( sum(size(filterSize) > 1) ~= 0 || filterSize <= 0 || mod(filterSize, 1) ~= 0 )        error('MOMEDA:InvalidInput', 'Input argument filterSize must be a positive integer scalar.')    elseif( sum(size(window) > 1) > 1 )        error('MOMEDA:InvalidInput', 'Input argument window must be 1d.')    elseif( period <= length(window) )        error('MOMEDA:InvalidInput', 'Period should be larger than the length of the window.')    elseif( filterSize >= length(x) )        error('MOMEDA:InvalidInput', 'Input argument filterSize must be smaller than the length of input signal x.')    end        L = filterSize;    x = x(:); % A column vector        %%% Calculte X0 matrix    N = length(x);    X0 = zeros(L,N);        for( l =1:L )        if( l == 1 )            X0(l,1:N) = x(1:N);        else            X0(l,2:end) = X0(l-1, 1:end-1);        end    end                            % "valid" region only    X0 = X0(:,L:N-1);   % y = f*x where only valid x is used                        % y = Xm0'*x to get valid output signal        autocorr = X0*X0';    autocorr_inv = pinv(autocorr);        % Built the impulse train vector separated the by periods    t = zeros(N-L,1);    points{1} = 1:period:(size(X0,2)-1);    points{1} = round(points{1});    t(points{1},1) = 1;        % Apply the windowing function to the target vectors    t = filter(window, 1, t);        % Calculate the spectrum of optimal filters    f = autocorr_inv * X0 * t;
    % Calculate the spectrum of outputs    y = X0'*f;        % Calculate the spectrum of PKurt values for each output    MKurt = mkurt(y,t);        % Plot the result    if( plotMode > 0 )        figure;        subplot(3,1,1)        plot(x)        title('Input signal');        xlabel('Sample number');                subplot(3,1,2)        plot(y)        title('Output signal');        xlabel('Sample number');                subplot(3,1,3)        stem(f)        title('Designed optimal FIR filter');        xlabel('Sample number');    endend
function [result] = mkurt(x,target)    % This function simply calculates the summed kurtosis of the input    % signal, x, according to the target vector positions.    result = zeros(size(x,2),1);    for i = 1:size(x,2)        result(i) = ( (target(:,i).^4)'*(x(:,i).^4) )/(sum(x(:,i).^2)^2) * sum(abs(target(:,i)));    endend