当前位置:   article > 正文

MCKD_相关峭度_解卷积_源码(部分代码)_maximum correlated kurtosis deconvolution

maximum correlated kurtosis deconvolution
  1. function [y_final f_final ckIter] = mckd(x,filterSize,termIter,T,M,plotMode)
  2. %MAXIMUM CORRELATED KURTOSIS DECONVOLTUION
  3. % code and method by Geoff McDonald (glmcdona@gmail.com), May 2011
  4. % This code file is an external reference for a paper being submitted
  5. % for review.
  6. %
  7. % mckd(x,filterSize,termIter,plotMode,T,M)
  8. %
  9. % Description:
  10. % This method tries to deconvolve a periodic series of impulses from
  11. % a 1d vector. It does this by designing a FIR filter to maximize
  12. % a norm criterion called Correlated Kurtosis. This method is has
  13. % applications in fault detection of rotating machinery (such as
  14. % ball bearing and gear faults).
  15. %
  16. % Algorithm Reference:
  17. % (Paper link coming soon. If you are interested in this, please
  18. % contact me at glmcdona@gmail.com. I will add the link if/when the
  19. % paper is available online)
  20. %
  21. % Inputs:
  22. % x:
  23. % Signal to perform deconvolution on. This should be a 1d vector.
  24. % MCKD will be performed on this vector by designing a FIR
  25. % filter.
  26. %
  27. % filterSize:
  28. % This is the length of the finite impulse filter filter to
  29. % design. Using a value of around 100 is appropriate depending on
  30. % the data. Investigate the performance difference using
  31. % different values.
  32. %
  33. % termIter: (OPTIONAL)
  34. % This is the termination number of iterations. If the
  35. % the number of iterations exceeds this number, the MCKD process
  36. % will complete. Specify [] to use default value of 30.
  37. %
  38. % T:
  39. % This is the period for the deconvolution. The algorithm will
  40. % try to deconvolve periodic impulses separated by this period.
  41. % This period should be specified in number of samples and can be
  42. % fractional (such as 106.29). In the case of a fractional T, the
  43. % method will resample the data to the nearest larger integer T:
  44. % i.e. 106.29 -> 107
  45. % and the y_final output will still be at this resampled factor.
  46. %
  47. % M:
  48. % This is the shift order of the deconvolution algorithm.
  49. % Typically an integer value between 1 and 5 is good. Increasing
  50. % the number increases the number of periodic impulses it tries
  51. % to find in a row. For example M = 5 would try to extract at
  52. % least 5 impulses in a row. When you use a larger M you need a
  53. % better estimate of T. Using too large a M (approx M > 10) will
  54. % result in a loss of numerical precision.
  55. %
  56. % plotMode:
  57. % If this value is > 0, plots will be generated of the iterative
  58. % performance and of the resulting signal.
  59. %
  60. % Outputs:
  61. % y_final:
  62. % The input signal x filtered by the resulting MCKD filter.
  63. % This is obtained simply as: y_final = filter(f_final,1,x);
  64. %
  65. % f_final:
  66. % The final MCKD filter in finite impulse response format.
  67. %
  68. % ckIter:
  69. % Correlated Kurtosis of shift M according to MED iteration.
  70. % ckIter(end) is the final ck.
  71. %
  72. % Example Usage:
  73. % % We want to extract the periodic impulses
  74. % % from the very strong white noise!
  75. % n = 0:999;
  76. % x = 3*(mod(n,100)==0) + randn(size(n));
  77. % [y_final f_final ck_iter] = mckd(x,400,30,100,7,1); % M = 7
  78. % % T = 100
  79. %
  80. %
  81. % Note:
  82. % The solution is not guaranteed to be the optimal solution to the
  83. % correlated kurtosis maximization problem, the solution is just a
  84. % local maximum and therefore a good pick.
  85. % Assign default values for inputs
  86. if( isempty(filterSize) )
  87. filterSize = 100;
  88. end
  89. if( isempty(termIter) )
  90. termIter = 30;
  91. end
  92. if( isempty(plotMode) )
  93. plotMode = 0;
  94. end
  95. if( isempty(T) )
  96. T = 0;
  97. end
  98. if( isempty(M) )
  99. M = 1;
  100. end
  101. % Validate the inputs
  102. if( sum( size(x) > 1 ) > 1 )
  103. error('MCKD:InvalidInput', 'Input signal x must be a 1d vector.')
  104. elseif( sum(size(T) > 1) ~= 0 || T < 0 )
  105. error('MCKD:InvalidInput', 'Input argument T must be a zero or positive scalar.')
  106. elseif( sum(size(termIter) > 1) ~= 0 || mod(termIter, 1) ~= 0 || termIter <= 0 )
  107. error('MCKD:InvalidInput', 'Input argument termIter must be a positive integer scalar.')
  108. elseif( sum(size(plotMode) > 1) ~= 0 )
  109. error('MCKD:InvalidInput', 'Input argument plotMode must be a scalar.')
  110. elseif( sum(size(filterSize) > 1) ~= 0 || filterSize <= 0 || mod(filterSize, 1) ~= 0 )
  111. error('MCKD:InvalidInput', 'Input argument filterSize must be a positive integer scalar.')
  112. elseif( sum(size(M) > 1) ~= 0 || M < 1 || round(M)-M ~= 0 )
  113. error('MCKD:InvalidInput', 'Input argument M must be a positive integer scalar.')
  114. elseif( filterSize > length(x) )
  115. error('MCKD:InvalidInput', 'The length of the filter must be less than or equal to the length of the data.')
  116. end
  117. % Force x into a column vector
  118. x = x(:);
  119. L = filterSize;
  120. OriginalLength = length(x);
  121. % Perform a resampling of x to an integer period if required
  122. if( abs(round(T) - T) > 0.01 )
  123. % We need to resample x to an integer period
  124. T_new = ceil(T);
  125. % The rate transformation factor
  126. Factor = 20;
  127. % Calculate the resample factor
  128. P = round(T_new * Factor);
  129. Q = round(T * Factor);
  130. Common = gcd(P,Q);
  131. P = P / Common;
  132. Q = Q / Common;
  133. % Resample the input
  134. x = resample(x, P, Q);
  135. T = T_new;
  136. else
  137. T = round(T);
  138. end
  139. N = length(x);
  140. % Calculate XmT
  141. XmT = zeros(L,N,M+1);
  142. for( m = 0:M)
  143. for( l = 1:L )
  144. if( l == 1 )
  145. XmT(l,(m*T+1):end,m+1) = x(1:N-m*T);
  146. else
  147. XmT(l, 2:end,m+1) = XmT(l-1, 1:end-1,m+1);
  148. end
  149. end
  150. end
  151. % Calculate the matrix inverse section
  152. Xinv = inv(XmT(:,:,1)*XmT(:,:,1)');
  153. % Assume initial filter as a delayed impulse
  154. f = zeros(L,1);
  155. f(round(L/2)) = 1;
  156. f(round(L/2)+1) = -1;
  157. f_best = f;
  158. ck_best = 0;
  159. iter_best = 0;
  160. % Initialize iteration matrices
  161. y = zeros(N,1);
  162. b = zeros(L,1);
  163. ckIter = [];

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/知新_RL/article/detail/731084
推荐阅读
相关标签
  

闽ICP备14008679号