1.算法描述
K-SVD算法是一种新型的字典训练法，其基本原理是基于K-SVD算法改进所得到的，其主要过程是字典的训练过程，其具有非常好的自适应性能。该算法的整体流程图如下图所示：
基于KSVD字典学习的图像去噪方法，其可以克服冲击噪声中纹理细节丢失，图像突变干扰等影响因素。该算法的核心内容为设置字典D为DCT字典，然后采用KSVD算法对字典D的原子和相应系数矩阵进行更新，将更新后的D、相应系数矩阵和代入原始图像的估计公式，得到含噪图像的去噪结果。
从线性组合的角度来看，K-SVD训练算法的稀疏模型可表示为：
那么，K-SVD算法的基本原理为在一组基下，获得信号y的一个近似稀疏表示x，且x满足尽可能好的恢复信号x。其中公式的计算过程是一个迭代过程，当字典D固定的时候，通过OMP正交匹配追踪算法找到字典D上Y的稀疏表示的稀疏矩阵D，然后再根据系数矩阵X找到更好的字典D。
2.仿真效果预览
matlab2022a仿真结果如下：
对不同噪声密度的冲击噪声干扰的图像进行滤波测试，这里分别设置脉冲噪声密度nl大小为0.1，0.4，0.7，对测试后图像进行质量评价，评价函数分别采用峰值信噪比（Peak Signal to Noise Ratio，PSNR）和图像相似度（structural similarity，SSIM），其中PSNR计算公式如下：
其中Peak就是指8bits（n=8）表示法的最大值255。MSE指Mean Square Error。PSNR的单位为dB。所以PSNR值越大，就代表失真越少。
SSIM计算公式如下：
当两张影像其中一张为无失真影像，另一张为失真后的影像，二者的结构相似性可以看成是失真影像的影像品质衡量指标。
从上述仿真结果可知，随着噪声密度的增加，测试图像的滤波图像质量逐渐降低，当噪声密度为0.1的时候，滤波后的图像，其PSNR值大概在35左右附近，SSIM值在0.96附近；当噪声密度为0.4的时候，滤波后的图像，其PSNR值大概在29左右附近，SSIM值在0.9附近；当噪声密度为0.7的时候，滤波后的图像，其PSNR值大概在25左右附近，SSIM值在0.78附近。即随着噪声的增加，滤波后的图像质量逐渐降低。
3.MATLAB核心程序
iterNum
% iterNum
% find the coefficients one by one
% ===remove the corrupted pixels and make some change on Dictionary==
for kk = 1:1:size(Data,2)
% kk
tempon = Ionblocks(:,kk);
location = find(tempon==1); %find locations of the good points
tempblock = Data(location,kk);
tempD = Dictionary(location,:);
tempD = tempD./repmat(sqrt(sum(tempD.^2)), size(tempD, 1), 1);
if (Reduce_DC)
vecOfMeans = mean(tempblock);
tempblock = tempblock - vecOfMeans;
end
%OMP压缩感知
if WS == 1
[sols, iters, activationHist] = WOMP(tempD, tempblock, size(Dictionary,2),64/2, 0, 0, 0, 1e-5,errT);
CoefMatrix(:,kk) = sols;
else
CoefMatrix(:,kk) = OMP(tempD,tempblock,L);
end
end
%===================================================================
replacedVectorCounter = 0;
rPerm = randperm(size(Dictionary,2));
for j = rPerm
%优化最佳的字典元素
[betterDictionaryElement,CoefMatrix,addedNewVector] = Modified_I_findBetterDictionaryElement(Data,Ionblocks,[FixedDictionaryElement,Dictionary],j+size(FixedDictionaryElement,2),CoefMatrix ,param.L);
Dictionary(:,j) = betterDictionaryElement;
if (param.preserveDCAtom)
tmpCoef = FixedDictionaryElement\betterDictionaryElement;
Dictionary(:,j) = betterDictionaryElement - FixedDictionaryElement*tmpCoef;
Dictionary(:,j) = Dictionary(:,j)./sqrt(Dictionary(:,j)'*Dictionary(:,j));
end
replacedVectorCounter = replacedVectorCounter+addedNewVector;
end
if (iterNum>=1 & param.displayProgress)
if (param.errorFlag==0)
output.totalerr(iterNum) = sqrt(sum(sum((Data-[FixedDictionaryElement,Dictionary]*CoefMatrix).^2))/prod(size(Data)));
disp(['Iteration ',num2str(iterNum),' Total error is: ',num2str(output.totalerr(iterNum))]);
else
output.numCoef(iterNum) = length(find(CoefMatrix))/size(Data,2);
disp(['Iteration ',num2str(iterNum),' Average number of coefficients: ',num2str(output.numCoef(iterNum))]);
end
end
%清除原来的旧的字典
Dictionary = Modified_I_clearDictionary(Dictionary,CoefMatrix(size(FixedDictionaryElement,2)+1:end,:),Data,Ionblocks);
if (isfield(param,'waitBarHandle'))
waitbar(iterNum/param.counterForWaitBar);
end
end
output.CoefMatrix = CoefMatrix;
Dictionary = [FixedDictionaryElement,Dictionary];SEL = 0;
%冲击密度
nl = 0.5;
%************************************************************************
%读取图像
Is = double(imread('Pic\1.png'));
%由于仿真速度较慢，将图像大小统一变为256*256大小的图像
Is = imresize(Is,[256,256]);
%下面判断图像是否是灰度图
[R,C,K] = size(Is);
if K == 1
I0 = Is;
else
I0 = Is(:,:,1);
end
%加入冲击噪声
[I1,Ion] = impulsenoise(I0,nl,'random-value');
%进行滤波
%字典搜索参数
sigma = 5;
Th = 5;
flagInitialD= 1;
flagOMP = 1;
%去噪声滤波
I2 = medfilt2(I1,[3,3]);
%计算PSNR指标
psnrIdenoise = PSNRs(I2,I0);
%计算SSIM指标
[mssim,ssim_map] = ssim_index(I2,I0);
psnrIdenoise
mssim