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
