目录

1.算法描述

2.仿真效果预览

3.MATLAB核心程序

4.完整MATLAB

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核心程序

for iterNum = 1:param.numIteration
    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
05_036_m

4.完整MATLAB

V