ssnf.m

编程实现基于小波变换尺度相关性去噪算法

function [s1 a d] = SSNF(s, n, h, g, g1, picDir)
%UNTITLED1 Summary of this function goes here
%  Detailed explanation goes here

% picDir = 'pic';

Ns = length(s);

%小波分解
[a d hj gj] = swt1d_decomp(s,n,h,g);

% figure, title('细节信息矩阵');
% hold on;
% for i=1:n
%     subplot(n,1,i);
%     title(['细节矩阵 ' num2str(i)]);
%     plot(d(:,i));
%     set(gca,'XTick',0:50:Ns);
% end
% saveas(gcf, [picDir filesep 'dwithNoise'], 'png'), close(gcf);

%去除噪声
w = d';
ww = w;
wfmap = zeros(size(w));
%相关系数
cor = w(1:n-1,:).*w(2:n,:);
%用于求规范化的相关系数
pw = sum(w.*w,2);
pcor = sum(cor.*cor,2);
ncor = cor.*repmat(sqrt(pw(1:n-1)./pcor),1,Ns);

% figure, title('相关系数');
% hold on;
% for i=1:n-1
%     subplot(n-1,1,i);
%     title(['相关系数 ' num2str(i)]);
%     plot(ncor(i,:));
%     set(gca,'XTick',0:50:Ns);
% end
% saveas(gcf, [picDir filesep 'ncor'], 'png'), close(gcf);
% 
% figure, title('相关系数 - 小波系数');
% hold on;
% for i=1:n-1
%     subplot(n-1,1,i);
%     title(['相关系数 - 小波系数 ' num2str(i)]);
%     plot(abs(ncor(i,:))-abs(w(i,:)));
%     set(gca,'XTick',0:50:Ns);
% end
% saveas(gcf, [picDir filesep 'ncorminusw'], 'png'), close(gcf);


noiseSigma = -1;
hconv = hj{1};   %hconv = h0*h1*h2...*hj-2
for j=1:n-1
    k = 0;
    while(1)
        k = k+1;
    %提取边缘信息
    edgeIndex = find(abs(ncor(j,:))>=abs(ww(j,:)));
    edgeIndex = edgeIndex(find(wfmap(j,edgeIndex)==0));
    wfmap(j,edgeIndex) = 1;
    Ks = length(edgeIndex);
    %去除已经提取的小波系数
    ww(j,edgeIndex) = 0;
    cor(j,edgeIndex) = 0;
    %更新规范化相关系数
    pw(j) = sum(ww(j,:).*ww(j,:));
    pcor(j) = sum(cor(j,:).*cor(j,:));
    if(pw(j)>0)
        ncor(j,:) = cor(j,:)*sqrt(pw(j)/pcor(j));
    end
    
%     figure;
%     hold on;
%     title(['提出的小波系数 level' num2str(j) '(' num2str(k) ')']);
%     wj = zeros(1,Ns);
%     wj(find(wfmap(j,:)==1)) = w(j,find(wfmap(j,:)==1));
%     plot(wj), set(gca,'XTick',0:50:Ns);
%     saveas(gcf, [picDir filesep 'wj_' num2str(j) '_' num2str(k)],'png');
%     close(gcf);
    
%     figure;
%     wj1 = zeros(1,Ns);
%     wj1(find(wfmap(j,:)==0)) = w(j,find(wfmap(j,:)==0));    
%     subplot(3,1,1)
%     hold on;
%     plot(wj1), set(gca,'XTick',0:50:Ns);
%     title(['剩下的小波系数 level' num2str(j) '(' num2str(k) ')']);
%     subplot(3,1,2)
%     hold on;
%     title(['剩下的相关系数 level' num2str(j) '(' num2str(k) ')']);
%     plot(ncor(j,:)), set(gca,'XTick',0:50:Ns);
%     subplot(3,1,3)
%     hold on;
%     title(['剩下的相关系数 - 小波系数 level' num2str(j) '(' num2str(k) ')']);
%     plot(abs(ncor(j,:))-abs(wj1)), set(gca,'XTick',0:50:Ns);
%     saveas(gcf, [picDir filesep 'wj_left_' num2str(j) '_' num2str(k)],'png');
%     close(gcf);
    
    if(noiseSigma == -1)
    %估计噪声
    noiseSigma = sqrt(pw(j)/(Ns-Ks))/sqrt(sum(g.*g));
    end
    
    %估计当前的噪声
    if(k==1)
        if(j==1)
            nsj = noiseSigma*sqrt(sum(g.*g));
%             g
        else 
%             gj{j};
            hconvj = conv(hconv,gj{j});
            nsj = noiseSigma*sqrt(sum(hconvj.*hconvj));
        end
    end
     
    %判断是否循环
    %符合能量判别条件或者该分辨率下所有小波系数都已提取出来
%     )
%     find(wfmap(j,:))
%     pwj = pw(j)
%     ne = (Ns-Ks)*power(nsj,2)
    if((pw(j)<=(Ns-Ks)*power(nsj,2))||(length(find(wfmap(j,:)==0))==0)||(Ks==0))
        break;
    end
    end
    hconv = conv(hconv,hj{j});
end
%使用wfmap去除噪声系数
w(find(wfmap==0))=0;

d(:,1:n-1) = w(1:n-1,:)';

% figure, title('细节信息矩阵');
% hold on;
% for i=1:n
%     subplot(n,1,i);
% %     title(['去噪后细节矩阵 ' num2str(i)]);
%     plot(d(:,i));
%     set(gca,'XTick',0:50:Ns);
% end
% saveas(gcf, [picDir filesep 'dwithoutNoise'], 'png'), close(gcf);

%小波重构
s1 = swt1d_recon(a,d,n,h,g1);

s1 = s1(:,n+1);