fastica去噪matlab_MATLAB--数字图像处理频域图像分析

- 日理万妓 2023-01-02 03:29 208阅读 0赞

![2b06628abd864118ab935452e5118bdb.png][]

## 一、实验名称 ##

频域图像分析

## 二、实验目的 ##

1.熟悉MATLAB软件的使用。  
2.掌握频域图像分析的原理及数学运算。

## 三、实验内容 ##

1.自选一幅图像，并对其分别添加一定强度的周期噪声和高斯噪声，然后分别采用高斯模板、中值滤波的时域方法以及傅里叶变换和小波变换的频率滤波方法对该含噪图像进行去噪处理，并基于PSNR值和视觉效果这两个指标来比较这四种滤波方法对两种不同噪声的去噪能力。  
2.编写一个程序，要求实现下列算法:首先将閣像分割为8x8的子图像，对每个予图像进行FFT.对每个了图像中的64个系数。按照每个系数的方差来排序后，舍去小的变換系數，只保留16个系数，实现4: I的图像压缩。  
3.给定一幅行和列都为2的整数次幕图像，用Haar小波基函数对其进行二维小波变换，试着将最低尺度近似分量置零再反变换，结果是什么?如果把垂直方向的细节分量置零，反变换后结果又是什么呢?试解释一下原因。  
4.基于小波变换对图像进行不同压缩比的压缩。在同压缩比情况下，对于基于小波变换和基于傅里叶变换的压缩结果，比较=二者保留原图像能里百分比情况。

## 四、实验仪器与设备 ##

Win10 64位电脑  
MATLAB R2017a

## 五、实验原理 ##

1.傅里叶变换  
从纯粹的数学意义上看，傅里叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数。  
傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数傅里叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，通常用一个二维矩阵表示空间上各点，记为z=f(x,y)。又因空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就必须由梯度来表示，这样我们才能通过观察图像得知物体在三维空间中的对应关系。  
2.小波变换  
小波变换是时间（空间）频率的局部化分析，它通过伸缩平移运算对信号(函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节带噪声信号经过预处理，然后利用小波变换把信号分解到各尺度中，在每一尺度下把属于噪声的小波系数去掉，保留并增强属于信号的小波系数，最后再经过小波逆变换回复检测信号。  
小波变换在去除噪声时可提取并保存对视觉起主要作用的边缘信息，而传统的基于傅里叶变换去除噪声的方法在去除噪声和边沿保持上存在着矛盾，因为傅里叶变换方法在时域不能局部化，难以检测到局域突变信号，在去除噪声的同时，也损失了图像边沿信息。由此可知，与傅里叶变换去除噪声的方法相比较，小波变换法去除噪声具有明显的性能优势。  
3.PSNR算法  
peak的中文意思是顶点。而ratio的意思是比率或比列的。整个意思就是到达噪音比率的顶点信号，psnr一般是用于最大值信号和背景噪音之间的一个工程项目。通常在经过影像压缩之后，输出的影像都会在某种程度与原始影像不同。为了衡量经过处理后的影像品质，我们通常会参考PSNR值来衡量某个处理程序能否令人满意。它是原图像与被处理图像之间的均方误差相对于(2n-1)2的对数值(信号最大值的平方，n是每个采样值的比特数)，它的单位是dB。

## 六、实验过程及代码 ##

t=imread('a1.jpg');
    
    %添加高斯噪声
    t1=imnoise(t,'gaussian',0,0.02);
    [m,n]=size(t);
    t2=t;
    
    %添加周期噪声
    for i=1:m
    for j=1:n
    t2(i,j)=t(i,j)+30*sin(30*i)+30*sin(30*j);
    end
    end
    imshow(t1),title('加入高斯噪声后')
    figure,imshow(t2),title('加入周期噪声后');
    
    %进行高斯滤波、中值滤波
    t3_1=t;
    t3_1(:,:,1)=medfilt2(t1(:,:,1),[3 3]);
    t3_1(:,:,2)=medfilt2(t1(:,:,2),[3 3]);
    t3_1(:,:,3)=medfilt2(t1(:,:,3),[3 3]);
    figure,imshow(t3_1),title('对高斯噪声进行中值滤波','fontsize',16)
    h=fspecial('gaussian',5,3);%确定滤波方式为高斯滤波 5是模板大小  3是方差
    t3_2=imfilter(t2,h);%滤波操作
    figure,imshow(t3_2),title('对周期噪声进行高斯滤波','fontsize',16)
    
    t3_1=t;
    t3_1(:,:,1)=medfilt2(t2(:,:,1),[3 3]);
    t3_1(:,:,2)=medfilt2(t2(:,:,2),[3 3]);
    t3_1(:,:,3)=medfilt2(t2(:,:,3),[3 3]);
    figure,imshow(t3_1),title('对高斯噪声进行中值滤波','fontsize',16)
    h=fspecial('gaussian',5,3);%确定滤波方式为高斯滤波 5是模板大小  3是方差
    t3_2=imfilter(t1,h);%滤波操作
    figure,imshow(t3_2),title('对周期噪声进行高斯滤波','fontsize',16)
    
    
    %对周期噪声傅里叶变换滤波
    t_f=rgb2gray(t2); %将图像灰度化
    t_f=fft2(double(t_f));%利用fft2()函数将图像从时域空间转换到频域空间
    t_f=fftshift(t_f);%将零频平移到中心位置
    [m,n]=size(t_f);
    m_min=round(m/2);
    n_min=round(n/2);
    t_rf=t_f;
    d_0=100;%设置阈值
    for i=1:m
    for j=1:n
    d_1=sqrt((i-m_min)^2+(j-n_min)^2);
    if(d_1>d_0)
    x=0;
    else 
    x=1;
    end
    t_rf(i,j)=x*t_f(i,j);
    end
    end
    t_rf=ifftshift(t_rf);
    t_rf=uint8(real(ifft2(t_rf)));
    figure,imshow(t_rf),title('对周期噪声进行傅里叶变换滤波后')
    
    %对高斯噪声进行傅里叶变换滤波
    t_f=rgb2gray(t1);
    t_f=fft2(double(t_f));
    t_f=fftshift(t_f);
    [m,n]=size(t_f);
    m_min=round(m/2);
    n_min=round(n/2);
    t_rf=t_f;
    d_0=100;
    for i=1:m
    for j=1:n
    d_1=sqrt((i-m_min)^2+(j-n_min)^2);
    if(d_1>d_0)
    x=0;
    else 
    x=1;
    end
    t_rf(i,j)=x*t_f(i,j);
    end
    end
    t_rf=ifftshift(t_rf);
    t_rf=uint8(real(ifft2(t_rf)));
    figure,imshow(t_rf),title('对高斯噪声进行傅里叶变换滤波后')
    
    %对周期噪声进行小波变换滤波
    [c,s]=wavedec2(t2,2,'sym4');  
    a1=wrcoef2('a',c,s,'sym4');    
     figure,imshow(uint8(a1)); title('对周期噪声进行小波变换滤波')
    %对高斯噪声进行小波变换滤波
    [c,s]=wavedec2(t1,2,'sym4');  
    a1=wrcoef2('a',c,s,'sym4');    
    figure, imshow(uint8(a1)); title('对高斯噪声进行小波变换滤波')

SNRP算法

function PSNR = psnr(f1, f2)
    %计算两幅图像的峰值信噪比
    k = 8;
    %k为图像是表示地个像素点所用的二进制位数，即位深。
    fmax = 2.^k - 1;
    a = fmax.^2;
    MSE =(double(im2uint8(f1)) -double( im2uint8(f2))).^2;
    b = mean(mean(MSE));
    PSNR = 10*log10(a/b);
    
    对图像进行4:1的压缩
     t=imread('a6.jpg');
    t=rgb2gray(t);%灰度化
    [k,p]=size(t);
    t=double(t)/255;%归一化 便于计算
    %显示原图
    imshow(t),title('原图','fontsize',16);
    %利用blkproc 进行分块 并对每一块进行fft操作
    t_fft=blkproc(t,[8 8],'fft2(x)');
    %利用im2col进行优化操作 便于计算
    t_block=im2col(t_fft,[8 8],'distinct');
    [t_change,ix]=sort(t_block);%对每一块图像进行排序
    [m,n]=size(t_block);
    nums=48;
    %对后48位系数清零
    for i=1:n
     t_block(ix(1:nums),i)=0;
    end
    t_rchange=col2im(t_block,[8 8],[k p],'distinct');
    t_ifft=blkproc(t_rchange,[8 8],'ifft2(x)');%对每一块进行傅里叶反变换
    figure,
    imshow(t_ifft),title('4:1压缩后','fontsize',16);

haar基函数进行小波变换

%实现一维Haar变换
    function [h]=D1Haar(f,N)
    J=log2(N);
    Y=zeros(J,N);
    f1=zeros(1,N);
    H=zeros(1,N);
    q=f';
    
    %调用倒序函数
    f1=Reverse(q,N);
    for A=1:N
        Y(1,A)=f1(A);
    end
    
    %第一层的迭代
    for C=1:N/2
        Y(2,C)=Y(1,C)+Y(1,C+N/2);
        Y(2,C+N/2)=Y(1,C)-Y(1,C+N/2);
    end
    
    %余下层的迭代
    M=0;
    for B=2:J
        K=2*N/(2^B); 
        F=zeros(1,K);
        for C=1:K/2
            Y(B+1,C)=Y(B,C)+Y(B,C+K/2);
            Y(B+1,C+K/2)=Y(B,C)-Y(B,C+K/2);
        end
        for C=K+1:2*K
            F(1,C-K)=Y(B,C);
        end
        Z=Reverse(F,K);
        for D=1:K
            Y(B+1,D+K)=Z(D);
        end
        for C=2*K+1:N
            Y(B+1,C)=Y(B,C);           
        end
    end
    
    %系数修正
    H(1,1)=(1/N)*Y(J+1,1);
    H(1,2)=(1/N)*Y(J+1,2);
    for a=2:J
        b=2^(a-1);
        c=b^(1/2);
            for d=b+1:2*b
            H(1,d)=(c/N)*Y(J+1,d);
        end
    end
    h=H';
    
    %二维哈尔函数的正变换
    function [y]=D2Haar(x,M,N)
    %先进行按列变换
    for i=1:N
        z=zeros(M,1);
        for c=1:M
           z(c,1)=x(c,i);
        end
        z1=D1Haar(z,M);%调用一维哈尔变换函数D1Haar
        for c=1:M
            x1(c,i)=z1(c,1);
        end
    end
    
    %再进行按行变换
    x2=x1';%将变换好的矩阵进行转置
    for j=1:M
        z=zeros(N,1);
        for c=1:N
            z(c,1)=x2(c,j);
        end
        z1=D1Haar(z,N);%调用一维哈尔变换函数D1Haar
        for c=1:N
            y1(c,j)=z1(c,1);
        end
    end
    y=y1';%再最后将变换好的矩阵作转置
    
    
         %对矩阵作一次哈尔变换
    function [y]=D2Har(x,M,N)
    for i=1:N
        z=zeros(M,1);
        for c=1:M
           z(c,1)=x(c,i);
        end
        z1=D1Haar(z,M);
        for c=1:M
            x1(c,i)=z1(c,1);
        end
    end
    y=x1;
    
    %哈尔函数的逆变换
    function [y]=ID2Haar(X,M,N)
    har1=zeros(M);
    I1=eye(M);%产生M维单位阵
    har1=D2Har(I1,M,M);%调用D2Har产生M维哈尔矩阵
    har2=M*har1';
    I2=eye(N);%产生N维单位阵
    har3=zeros(N);
    har3=D2Har(I2,N,N);%调用D2Har产生M维哈尔矩阵
    har4=N*har3;
    y=har2*X*har4;%哈尔逆变换
    
    %求倒叙函数
    function [F]=Reverse(F,M)
    N=log10(M)/log10(2);
    Y=zeros(1,M);
    for  x=0:M-1
         A=dec2bin(x,N);%十进制转二进制
         B=fliplr(A);%二进制倒叙
         C=bin2dec(B);%二进制转十进制
         Y(x+1)=F(C+1);
    end
    for x=0:M-1
         F(x+1)=Y(x+1);
    end  
    
        %对图片进行哈尔正变化，并对其进行恢复
    clear;
    M=256;
    N=256;
    
    %读取图象
    X=imread('C:UsersLiCongliangDesktop数字图像处理数字图像处理第七次作业tea.jpg');
    subplot(2,2,1);
    imshow(X);
    title('原图像');
    
    %缩放原图象
    x=imresize(X,[M,N]);%将原图象缩放成分辨率为256*256的图象
    subplot(2,2,2);
    imshow(x);
    title('缩放图象');
    
    %对缩放图象进行Haar变换
    y=D2Haar(x,M,N);
    subplot(2,2,3);
    imshow(y);
    title('变换图象');
    
    %对变换后的图象进行Harr逆变换，恢复原图象
    z=ID2Haar(y,M,N);
    %Z=imresize(z,[480,360]);
    subplot(2,2,4);
    imshow(z,[0,255]);
    title('恢复图象');
    
    小波变换进行图像压缩
    X=imread('a5.jpg');
    X=rgb2gray(X);
    subplot(221); imshow(X);
    title('原始图像');
    %对图像用小波进行层小波分解
    [c,s]=wavedec2(X,2,'haar');
    %提取小波分解结构中的一层的低频系数和高频系数
    cal=appcoef2(c,s,'haar',1);
    ch1=detcoef2('h',c,s,1);      %水平方向
    cv1=detcoef2('v',c,s,1);      %垂直方向
    cd1=detcoef2('d',c,s,1);      %斜线方向
    %各频率成份重构
    a1=wrcoef2('a',c,s,'haar',1);
    h1=wrcoef2('h',c,s,'haar',1);
    v1=wrcoef2('v',c,s,'haar',1);
    d1=wrcoef2('d',c,s,'haar',1);
    c1=[a1,h1;v1,d1];
    subplot(222),imshow(c1,[]);
    title ('分解后低频和高频信息');
    %进行图像压缩
    %保留小波分解第一层低频信息
    %首先对第一层信息进行量化编码
    ca1=appcoef2(c,s,'haar',1);
    ca1=wcodemat(ca1,440,'mat',0);
    %改变图像高度并显示
    ca1=0.5*ca1;
    subplot(223);imshow(cal,[]);
    title('第一次压缩图像');
    %保留小波分解第二层低频信息进行压缩
    ca2=appcoef2(c,s,'haar',2);
    %首先对第二层信息进行量化编码
    ca2=wcodemat(ca2,440,'mat',0);
    %改变图像高度并显示
    ca2=0.25*ca2;
    subplot(224);imshow(ca2,[]);
    title('第二次压缩图像');

七、实验结果与分析

图 1原图

![d43754d827feed1b18c119bffc292b4d.png][]

1.加入周期噪声、高斯噪声

![3f7652daf330df271717d6f421588307.png][]

![c72654e698cc36e105006b14db02f049.png][]

2.对添加了高斯噪声和周期噪声的图像进行高斯滤波

![f6152afbaaa3e0eca75ca1155cd57160.png][]

![92df4e655343090ad686d9aa270cebf4.png][]

PSNR值  
1.对高斯噪声进行高斯滤波后 23.0287  
2.对周期噪声进行高斯滤波后 23.4837

2.中值滤波

![38869b9549fd46dfee2b7edf3280a3dc.png][]

![0ebc8dd9e949ea4238892b2fb3927b1a.png][]

PSNR值：  
1.对高斯噪声进行中值滤波 23.9931  
2.对周期噪声进行中值滤波 24.3134

3.傅里叶变换滤波

![395a834911aa1e2bea81c3ce18df28c4.png][]

![5d7d73051d34b92ad6b602e9a9ae66de.png][]

PSNR值：  
1.对添加了高斯噪声的图像进行傅里叶变换滤波 20.4922  
2.对添加了周期噪声的图像进行傅里叶变换滤波 18.9736

4.小波变换滤波

![0e450b0ae1999fbee891d6e7db90bf2c.png][]

![1612764a48f6ed4e46a62e02809aab89.png][]

PSNR值：  
1.对添加了高斯噪声的图像进行小波变换滤波 23.4712  
2.对添加了周期噪声的图像进行小波变换滤波 24.4525

分析：  
对于高斯噪声，高斯滤波和傅里叶变换滤波声的除噪效果较好，中值滤波效果较差，小波变换滤波的处理效果也比较好  
对于周期噪声，中值滤波和高斯滤波效果不是很好，傅里叶变换变换滤波对噪声的去处效果比较好，对于原图像损坏不大，小波变换对原图的损坏较大，但是图片可以看出噪声也去除的比较好。

5.图像压缩(4:1压缩) 原图-左 压缩后-右

![3f95fa98e14862bf9d4a8a60c37c8d3a.png][]

![791a3d8b7e5875e3dac3a8b38101aa1a.png][]

![fa84465a1bb912eb04eed0546e1284e0.png][]

![749c7ceb94c549bf5d514659f0fa6224.png][]

分析：  
图像压缩算法就是先将一副图像分成很多小块，然后分别对这些小块进行变换，这里采用的是傅里叶变换，然后过滤掉冗余的像素点，然后再利用反变换得到压缩后的图像即可。  
小波变换  
1.定义  
小波变换是时间（空间）频率的局部化分析，它通过伸缩平移运算对信号(函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节带噪声信号经过预处理，然后利用小波变换把信号分解到各尺度中，在每一尺度下把属于噪声的小波系数去掉，保留并增强属于信号的小波系数，最后再经过小波逆变换回复检测信号。  
2.优点  
小波变换在去除噪声时可提取并保存对视觉起主要作用的边缘信息，而传统的基于傅里叶变换去除噪声的方法在去除噪声和边沿保持上存在着矛盾，因为傅里叶变换方法在时域不能局部化，难以检测到局域突变信号，在去除噪声的同时，也损失了图像边沿信息。由此可知，与傅里叶变换去除噪声的方法相比较，小波变换法去除噪声具有明显的性能优势。  
Haar基函数进行小波变换  
图 2原图

![214c08b577ef7b013a65cf097bf3ce8b.png][]

图 3 haar变换

![6f27dfc7c32f3297c040dd3979349223.png][]

图 4 haar反变换后

![e6307dfa24e1c8b7d113db4faca644cf.png][]

图 5 最低分量近似置零

![edfead315c2f2072ba5051030d1272cb.png][]

图 6 垂直分量置零

![e678fd6aa15e5823d58d887a74efa3d9.png][]

小波变换进行图像压缩与傅里叶变换压缩对比

1.压缩比 1:2(左-小波压缩 右-傅里叶压缩)

![3b9e3b278d5d659cabc4b46e4599fb68.png][]

![3fd1b447b02513e13d2a56bfcde892ed.png][]

2.压缩比 1:4(左-小波压缩 右-傅里叶压缩)

![bd3bb151db86af8eefb6ec37fa6730a8.png][]

![6e44d92ba57c47681041d645698ed202.png][]

## 八、实验总结及心得体会 ##

通过这次实验，学到了很多。特别是在傅里叶变换和小波变换等方面，开始的时候连傅里叶变换的基础基础也不懂，后来在csdn上看了一篇讲解傅里叶变换的文章，豁然开朗，傅里叶变换居然可以将一个时域信号转化到频域，而且自己还对与i有了更加深刻的理解。虽然傅里叶变换可以把信号从时域转换到频域，但是频域与时域的对应关系却无法一一对应，所以诞生了小波变换。小波变换的特别之处就是可以把一个时域上的信息转换为时域-频域一一对应，这对应特殊信号的提取是有很好的效果，在一定程度上比傅里叶变换更厉害。但是在傅里叶、小波等基础概念知识方面，自己还是涉猎的比较少，原理的论证公式太复杂了。

![3a38fa719d697bd0012f21c0b852afc5.png][]

[2b06628abd864118ab935452e5118bdb.png]: /images/20221120/3d8a5c47eca7499796379d61336f1213.png
[d43754d827feed1b18c119bffc292b4d.png]: /images/20221120/1345901b509d4b63b434ec13eb1de048.png
[3f7652daf330df271717d6f421588307.png]: /images/20221120/ed5d7897b3fd4bf7bd2b11054a381c6a.png
[c72654e698cc36e105006b14db02f049.png]: /images/20221120/580bb0c1c9c64cd1ab4f26eb9be13b9b.png
[f6152afbaaa3e0eca75ca1155cd57160.png]: /images/20221120/14e68bbd37e1472fbb9e0b2d5c7fd6a6.png
[92df4e655343090ad686d9aa270cebf4.png]: /images/20221120/f1af906848f041289b23b39285017607.png
[38869b9549fd46dfee2b7edf3280a3dc.png]: /images/20221120/67ec87b39c4d42af92a6586d9a33a2b1.png
[0ebc8dd9e949ea4238892b2fb3927b1a.png]: /images/20221120/b74028eb0ee04576be664fb9b8932589.png
[395a834911aa1e2bea81c3ce18df28c4.png]: /images/20221120/f6ada5e163ac4d6a83354b3038b356f5.png
[5d7d73051d34b92ad6b602e9a9ae66de.png]: /images/20221120/146908763e324baab49ba98c6c372a37.png
[0e450b0ae1999fbee891d6e7db90bf2c.png]: /images/20221120/63412cdbdbcf445cb172e561e35fa81e.png
[1612764a48f6ed4e46a62e02809aab89.png]: /images/20221120/2b6dd86643534d7ca7cec2d695cce941.png
[3f95fa98e14862bf9d4a8a60c37c8d3a.png]: /images/20221120/10bd050b7c734692807de1cbbfe9fdef.png
[791a3d8b7e5875e3dac3a8b38101aa1a.png]: /images/20221120/8449b40f304d4003923bb04e7b2e900e.png
[fa84465a1bb912eb04eed0546e1284e0.png]: /images/20221120/184f080b9af641ad9210e7d53a5693ff.png
[749c7ceb94c549bf5d514659f0fa6224.png]: /images/20221120/27d5fbb23c354551ae921e1993793be5.png
[214c08b577ef7b013a65cf097bf3ce8b.png]: /images/20221120/19bd7a8b9da14ef49dee087c8a73c66d.png
[6f27dfc7c32f3297c040dd3979349223.png]: /images/20221120/665c35d4da1a4b4f92bbbea421d8f894.png
[e6307dfa24e1c8b7d113db4faca644cf.png]: /images/20221120/aaa03b6cbf0f4c3780b02eacbac0bf0a.png
[edfead315c2f2072ba5051030d1272cb.png]: /images/20221120/9960464e2bb1489290402663e79425f2.png
[e678fd6aa15e5823d58d887a74efa3d9.png]: /images/20221120/d8cc546de0f148e1b10f9a98ff3637e7.png
[3b9e3b278d5d659cabc4b46e4599fb68.png]: /images/20221120/1b03580988d74c7eb40b5fb07e7120de.png
[3fd1b447b02513e13d2a56bfcde892ed.png]: /images/20221120/c96dcedc10db473ea8635ecab392cfdc.png
[bd3bb151db86af8eefb6ec37fa6730a8.png]: /images/20221120/8259ad5b8cef4d8e9770c3aa59e1ca32.png
[6e44d92ba57c47681041d645698ed202.png]: /images/20221120/35e365b99a854700862bca4d2602adcc.png
[3a38fa719d697bd0012f21c0b852afc5.png]: /images/20221120/4096ec131729480696c9260873f0234f.png