最小熵反褶积（MED）matlab

function [y_final f_final kurtIter] = med(x,filterSize,termIter,termDelta,plotMode) if( sum(size(x)>1) == 1 ) x = x(:); end L = filterSize; autoCorr = zeros(1,L); for column = 1:size(x,2); for k = 0:L-1 x2 = zeros(size(x,1),1); x2(k+1:end) = x(1:end-k,column); % 计算自相关 autoCorr(k+1) = autoCorr(k+1) + sum(x(:,column).*x2); end end autoCorr = autoCorr / size(x,2); % 平均值归一化 A = toeplitz(autoCorr); A_inv = inv(A); % 矩阵初始化 f = zeros(L,1); y = zeros(size(x,1),size(x,2)); b = zeros(L,1); kurtIter = []; f(2) = 1; % 循环调整最小熵滤波器系数 n = 1; while n == 1 || ( n < termIter && ( (kurt(filter(f,1,x)) - kurtIter(n-1)) > termDelta ) ) % 输出信号 y = filter(f,1,x); % 计算峭度值 kurtIter(n) = kurt(y); % 计算矩阵 g = weighted av{ crosscorr( y.^3, x) } yc = y.^3; weightedCrossCorr = zeros(L,1); for column = 1:size(x,2); for k = 0:L-1 x2 = zeros(size(x,1),1); x2(k+1:end) = x(1:end-k,column); % 计算互相关 weightedCrossCorr(k+1) = weightedCrossCorr(k+1) + sum((y(:,column).^3).*x2); end end weightedCrossCorr = weightedCrossCorr / size(x,2); % 新的滤波器系数 f = A_inv*weightedCrossCorr; f = f/sqrt(sum(f.^2)); % 滤波结果标准化 n = n + 1; end %更新结果 f_final = f; y_final = filter(f_final,1,x); kurtIter(n) = kurt(y_final); % 绘图 if(plotMode>0) figure; subplot(2,1,1) plot(x) title('Input Signal(s)') ylabel('Value') xlabel('Sample Number') axis tight subplot(2,1,2) plot(y_final) title('Signal(s) Filtered by MED') ylabel('Value') xlabel('Sample Number') axis tight end end function [result] = kurt(x) %计算输入信号的峭度值 result = mean( (sum((x-ones(size(x,1),1)*mean(x)).^4)/(size(x,1)-2))./(std(x).^4) ); end