第一篇：MATLAB函数处理图像实现膨胀腐蚀
MATLAB函数处理图像实现膨胀腐蚀
一、实验目的
1、了解二值形态学的基本运算
2、掌握二值图像膨胀、腐蚀的基本方法
3、编程实现膨胀、腐蚀
二、实验要求
1、使用imdilate函数进行图像膨胀，并观察膨胀后图像的变化。
2、使用imerode函数进行图像腐蚀，观察腐蚀后的图像变化情况。
三、实验原理
膨胀：将与物体接触的所有背景点合并到该物体中，使边界向外部扩张的过程。利用它可以填补物体中的空洞。B对X膨胀所产生的二值图像D是满足以下条件的点（x,y）的集合:如果B的原点平移到点（x,y），那么它与X的交集非空。数学表达式：CAB
腐蚀：一种消除边界点，使边界向内部收缩的过程。利用它可以消除小而且无意义的物体。B对X腐蚀所产生的二值图像E是满足以下条件的点（x,y)的集合：如果B的原点平移到点（x,y），那么B将完全包含于X中。数学表达式：CAB
膨胀处理：一种消除边界点，使边界点向内部收缩的过程。
腐蚀处理：将与物体接触的所有背景点合并到该物体中，使边界向外部扩张的过程。
四、实验步骤
1.图像膨胀的Matlab实现：
可以使用imdilate函数进行图像膨胀,imdilate函数需要两个基本输入参数，即待处理的输入图像和结构元素对象。结构元素对象可以是strel函数返回的对象，也可以是一个自己定义的表示结构元素邻域的二进制矩阵。此外，imdilate还可以接受两个可选参数：PADOPT(padopt)——影响输出图片的大小、PACKOPT(packopt).——说明输入图像是否为打包的二值图像(二进制图像)。步骤1，首先创建一个包含矩形对象的二值图像矩阵。>> BW=zeros(9,10);>> BW(4:6,4:7)=1 BW = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 步骤2，使用一个3×3的正方形结构元素对象对创建的图像进行膨胀。>> SE=strel('square',3)SE = Flat STREL object containing 9 neighbors.Neighborhood: 1 1 1 1 1 1 1 1 1 步骤3，将图像BW和结构元素SE传递给imdilate函数。>> BW2=imdilate(BW,SE)BW2 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0
0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 步骤4，显示结果。>> imshow(BW,'notruesize')>> imshow(BW2,'notruesize')2.图像腐蚀的Matlab实现：
可以使用imerode函数进行图像腐蚀。imerode函数需要两个基本输入参数：待处理的输入图像以及结构元素对象。此外，imerode函数还可以接受3个可选参数：PADOPT(padopt)——影响输出图片的大小、PACKOPT(packopt).——说明输入图像是否为打包的二值图像(二进制图像)。M——指定原始图像的行数。以下程序示例说明了如何对某一副具体图像进行腐蚀操作，腐蚀前后的效果对比如图末。
步骤1，读取图像cameraman.tif（该图像是Matlab当前目录下自带的图片）>> BW1=imread('cameraman.tif');步骤2，创建一个任意形状的结构元素对象 >> SE=strel('arbitrary',eye(5));步骤3，以图像BW1和结构元素SE为参数调用imerode函数进行腐蚀操作。>> BW2=imerode(BW1,SE);步骤4，显示操作结果 >> imshow(BW1)>> figure,imshow(BW2)
五、实验代码及结果
代码：
imerode函数，该函数能够实现二值图像的腐蚀操作； imdilate函数，该函数能够实现二值图像的膨胀操作； bw=imread(‘d:image1’)bw=rgb2gray(bw)se1=strel(‘disk’,11);
se2=strel(‘line’,11,90);bw2= imdilate(bw,se2);bw1=imerode(bw,se1);imshow(bw),title(‘原图’)figure,imshow(bw2), title(‘膨胀后的图像’)figure,imshow(bw1), title(‘腐蚀后的图像’)
结果：
原
图
膨胀后的图像
腐蚀后的图像
六、实验心得体会
通过本次的实验，我了解了二值形态学的基本运算，掌握了二值图像膨胀、腐蚀的基本方法，并且会运用编程实现膨胀、腐蚀，本次的实验目的已经完成，意识到在以后的生活中要了解做事情的目的，注重每一个与细节，认真思考遇到的所有问题，提高自己各方面的能力。感谢尹强老师教会我们理论与实践知识，也让我明白了什么是学习，怎么样学习，为以后的生活奠定的基础与指引了方向。
第二篇：matlab数字图像处理 膨胀和腐蚀
基于Matlab的腐蚀和膨胀的边缘检测
一、实验目的: 掌握运用Matlab软件对灰度与二值图像的膨胀与腐蚀的处理方法。
二、实验环境（软件条件）：
Windws2000/XP
MATLAB 7.x
三、实验内容：
1、图像膨胀的Matlab实现 ① 实验原理：
膨胀：给图像中的对象边界添加像素。
在操作中，输出图像中所有给定像素的状态都是通过对输入图像的相应像素及邻域使用一定的规则进行确定。在膨胀操作时，输出像素值是输入图像相应像素邻域内所有像素的最大值。在二进制图像中，如果任何像素值为1，那么对应的输出像素值为1。
可以使用imdilate函数进行图像膨胀,imdilate函数需要两个基本输入参数，即待处理的输入图像和结构元素对象。结构元素对象可以是strel函数返回的对象，也可以是一个自己定义的表示结构元素邻域的二进制矩阵。此外，imdilate还可以接受两个可选参数：PADOPT(padopt)——影响输出图片的大小、PACKOPT(packopt).——说明输入图像是否为打包的二值图像(二进制图像)。② 实验步骤：
A、首先创建一个包含矩形对象的二值图像矩阵。R=zeros(9,10);R(4:6,4:7)=1 R = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 B、使用一个3×3的正方形结构元素对象对创建的图像进行膨胀。
C=strel('square',3)C = Flat STREL object containing 9 neighbors.Neighborhood: 1 1 1 1 1 1 1 1 1 C、将图像R和结构元素C传递给imdilate函数。R1=imdilate(R,C)R1 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 D、显示结果： imshow(R,'notruesize')imshow(R1,'notruesize')③ 实验结果： 膨胀前后效果图：
2、图像腐蚀的Matlab实现 ① 实验原理：
腐蚀：删除对象边界某些像素。
在操作中，输出图像中所有给定像素的状态都是通过对输入图像的相应像素及邻域使用一定的规则进行确定。在腐蚀操作中，输出像素值是输入图像相应像素邻域内所有像素的最小值。在二进制图像中，如果任何一个像素值为0，那么对应的输出像素值为0。
可以使用imerode函数进行图像腐蚀。imerode函数需要两个基本输入参数：待处理的输入图像以及结构元素对象。此外，imerode函数还可以接受3个可选参数：PADOPT(padopt)——影响输出图片的大小、PACKOPT(packopt).——说明输入图像是否为打包的二值图像(二进制图像)。M——指定原始图像的行数。② 实验内容及步骤：
A、读取图像cameraman.tif（该图像是Matlab当前目录下自带的图片）R=imread('cameraman.tif');B、创建一个任意形状的结构元素对象 C=strel('arbitrary',eye(5));C、以图像R和结构元素C为参数调用imerode函数进行腐蚀操作。R1=imerode(R,C);D、显示操作结果 imshow(R)； figure,imshow(C)； ③ 实验结果：
图像cameraman.tif 腐蚀前后的效果对比：
3、膨胀和腐蚀联合操作 ① 实验内容及步骤： A、创建结构元素： clear;close all C= strel('rectangle',[40 30]);
%结构元素必须具有适当的大小,既可以删电流线又可以删除矩形.B、使用结构元素腐蚀图像:
%将会删除所有直线,但也会缩减矩形 R=imread('circbw.tif');R1=imerode(R,C);imshow(R1);figure,imshow(R);C、恢复矩形为原有大小,使用相同的结构元素对腐蚀过的图像进行膨胀.R2=imdilate(R1,C);figure,imshow(R2);② 实验结果： 最终效果如下图: 原始图像 腐蚀后的图像 膨胀后的图像
4、边缘检测 ① 实验原理：
对于一副灰度二进制图像,如果图像像素值为1,则该像素的状态为ON,如果其像素值为0,则该像素的状态为OFF。在一副图像中，如果图像某个像素满足以下两个条件：
1.该像素状态为ON； 2.该像素邻域中有一个或多个像素状态为OFF。则认为该像素为边缘像素。
Matlab中提供了专门的函数bwperim,可以用于判断一副二进制图像中的哪些像素为边缘像素。② 实验内容及步骤：
以下程序代码示例就是利用bwperim函数，对图像circbw.tif进行边缘检测，其边缘像素检测效果如尾图。clear;close all R=imread('circbw.tif');R1=bwperim(R);
imshow(R)； figure,imshow(R1)； ② 实验结果：
心得体会：
本次实验利用形态学运算—腐蚀和膨胀对图像进行了处理。数学形态学的基本思想是：用具有一定形态的结构元素去量度和提取图像中的对应形状，以达到图像分析和识别的目的。它的优点有：有效滤除噪声、保留图像中原有信息、算法易于用并行处理方法有效实现等等。
总的来说，在此次实验中遇到了较多问题，经过对相关知识的复习以及对相关资料的查询，才逐个地解决了。在今后的实验课中要有充分的准备与预习，不懂的问题要提前查资料，这样在实验中才会胸有成竹，才不会手忙脚乱。
第三篇：matlab图像处理小结
1.function [center, r] = solve_circle(pt1, pt2, pt3)
2.%Effect: solve the circle which across points 'pt1', 'pt2' and 'pt3' 3.%Inputs:
4.%pt1, pt2, pt3: [x, y]
5.%center: the circle center [x0;y0] 6.%r: the radius of the circle 7.%Author: Su dongcai at 2012/1/2 8.A = zeros(2, 2);B = zeros(2, 1);9.[A(1, :), B(1)] = circle2line(pt1, pt2);10.[A(2, :), B(2)] = circle2line(pt2, pt3);11.center = AB;
12.r = norm(pt1'(y2^2 + y2^2)18.%(a-x2)^2 +(b-y2)^2 = r^2 | 19.%Inputs:
20.%pt1, pt2: [x1, y1], [x2, y2] 21.%Outputs:
22.%A: 2[x1-x2, y1-y2]
23.%B:(x1^2 + y1^2)pt2);
26.B = norm(pt1)^2-norm(pt2)^2;
close all;clear;clc;>> i=imread('rice.png');%>> imshow(i);>> background=imopen(i,strel('disk',15));>> i2=imsubtract(i,background);%>> figure,imshow(i2);>> i3=imadjust(i2,stretchlim(i2),[0 1]);%>> figure,imshow(i3);>> level=graythresh(i3);>> bw=im2bw(i3,level);%>> figure,imshow(bw);>> [labeled,numobjects]=bwlabel(bw,4);graindata=regionprops(labeled,'all');
close all;clear;clc;i=imread('rice.png');background=imopen(i,strel('disk',15));i2=imsubtract(i,background);i3=imadjust(i2,stretchlim(i2),[0 1]);level=graythresh(i3);bw=im2bw(i3,level);[labeled,numobjects]=bwlabel(bw,4);data=regionprops(labeled,'all');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.2 close all;clear;clc;>> i=imread('r.jpg');%>> figure,imshow(i);>> imgray=rgb2gray(i);>> figure,imshow(imgray)>> background=imopen(imgray,strel('disk',15));>> i2=imsubtract(imgray,background);%>> figure,imshow(i2);>> i3=imadjust(i2,stretchlim(i2),[0 1]);%>> figure,imshow(i3);>> level=graythresh(i3);>> bw=im2bw(i3,level);%>> figure,imshow(bw);>> imnobord=imclearborder(bw,4);%>> figure,imshow(imnobord);>> [labeled,numobjects]=bwlabel(bw,4);>> rgb_label=label2rgb(labeled,@spring,'c','shuffle');>> figure,imshow(rgb_label);>> graindata=regionprops(labeled,'all');hold on;for k=1:numobjects lab=sprintf('%d',k);text(graindata(k).Centroid(1),graindata(k).Centroid(2),lab,'Color','k');end hold off;%剔除碎米粒
>> idxdown=find([graindata.Area]<150);%剔除碎米粒 little=ismember(labeled,idxdown);figure,imshow(little);
[lab_little,num_little]=bwlabel(little,4);rgb_little=label2rgb(lab_little,@spring,'c','shuffle');figure,imshow(rgb_little);
little_data=regionprops(lab_little,'all');hold on;for k=1:num_little lab=sprintf('%d',k);text(little_data(k).Centroid(1),little_data(k).Centroid(2),lab,'Color','k');end hold off;%>> graindata(idxdown,:)=[];%剔除碎米粒 %剔除连接米粒
>> idxup=find([graindata.Area]>250);%剔除连接米粒 big=ismember(labeled,idxup);figure,imshow(big);
[lab_big,num_big]=bwlabel(big,4);rgb_big=label2rgb(lab_big,@spring,'c','shuffle');figure,imshow(rgb_big);
big_data=regionprops(lab_big,'all');hold on;for k=1:num_big lab=sprintf('%d',k);text(big_data(k).Centroid(1),big_data(k).Centroid(2),lab,'Color','k');end hold off;%>> graindata(numup,:)=[];%剔除连接米粒 %获取完整米粒
idxsuit=find([graindata.Area]>=150&[graindata.Area]<=250);suit=ismember(labeled,idxsuit);figure,imshow(suit);%获取完整米粒 [lab_suit,num_suit]=bwlabel(suit,4);suit_data=regionprops(lab_suit,'all');hold on;for k=1:num_suit signature=sprintf('%d',k);text(suit_data(k).Centroid(1),suit_data(k).Centroid(2),signature,'Color','r');end hold off;%获取完整米粒 whos graindata whos little_data whos big_data whos suit_data
>> graindata >> mean([graindata.Area])>> mean([graindata.Eccentricity])>> mean([graindata.MajorAxisLength])>> mean([graindata.MinorAxisLength])>> mean([graindata.EquivDiameter])>> figure,hist([graindata.Area],20);>> figure,hist([graindata.Eccentricity],20);>> figure,hist([graindata.MajorAxisLength],20);>> figure,hist([graindata.MinorAxisLength],20);>> figure,hist([graindata.EquivDiameter],20);
data=[graindata.Area] data=[graindata.Centroid] data=[graindata.BoundingBox] data=[graindata.SubarrayIdx] data=[graindata.MajorAxisLength] data=[graindata.MinorAxisLength] data=[graindata.Eccentricity] data=[graindata.Orientation] data=[graindata.ConvexHull] data=[graindata.ConvexImage] data=[graindata.ConvexArea] data=[graindata.Image] data=[graindata.FilledImage] data=[graindata.FilledArea] data=[graindata.EulerNumber] data=[graindata.Extrema] data=[graindata.EquivDiameter] data=[graindata.Solidity] data=[graindata.Extent] data=[graindata.PixelIdxList] data=[graindata.PixelList]
Area 计算各个连通区域中的象素总数 BoundingBox 包含相应区域的最小矩形 Centroid 给出每个区域的质心
MajorAxisLength 与区域具有相同标准二阶中心矩（又叫标准差）的椭圆的长轴长度 MinorAxisLength 与区域具有相同标准二阶中心矩的椭圆的短轴长度 Eccentricity 与区域具有相同标准二阶中心矩的椭圆的离心率
Orientation 与区域具有相同标准二阶中心矩的椭圆的长轴与x轴的交角 Image 二值图像，与某区域具有相同大小的逻辑矩阵。
FilledImage 与上相同，唯一区别是这是个做了填充的逻辑矩阵！本例中和上面的没有区别，只有 区域有空洞时才有明显差别。
FilledArea 是标量，填充区域图像中的 on 像素个数
ConvexHull 是p行2列的矩阵，包含某区域的最小凸多边形 ConvexImage 二值图像，用来画出上述的区域最小凸多边形 ConvexArea 是标量，填充区域凸多边形图像中的 on 像素个数 EulerNumber 等于图像中目标个数减去这些目标中空洞的个数 Extrema 8行2列矩阵，八方向区域极值点
EquivDiameter 是标量，等价直径：与区域具有相同面积的圆的直径.计算公式为：sqrt(4*Area/pi)
Solidity 是标量，同时在区域和其最小凸多边形中的像素比例。计算公式为： Area/ConvexArea，这也是个仿射特征，实际上反映出区域的固靠性程度。
Extent 是标量，同时在区域和其最小边界矩形中的像素比例。计算公式为：Area除以边界矩 形面积，这也是个仿射特征，实际上反映出区域的扩展范围程度。
PixelIdxList p元向量，存储区域像素的索引下标
PixelList p行ndims(L)列矩阵，存储上述索引对应的像素坐标 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 基于特定原则的区域选择
当你要基于特定准则条件选择某个区域时，将函数 ismember 和 regionprops 联合使用是很有用处的。例如：创建一个只包含面积大于80的二值图像，用以下命令
idx = find([stats.Area] > 80);BW2 = ismember(L,idx);regionprops函数的扩展思路
在regionprops函数的基础上，你可以使用它提供的基本数据来扩展它的功能，比如我就将区域的曲率数据和骨架数据作为它的另外属性值来开发，从而希望它能用来做更细致的特征提取。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.2 P221图像粒度测定（雪花）>> i=imread('snowflakes.png');>> figure,imshow(i);>> %(2)>> clahei=adapthisteq(i,'numtiles',[10 10]);>> clahei=imadjust(clahei);>> imshow(clahei);>> gi=imadjust(im2double(i),[],[0 1]);>> figure,imshow(gi),title('adjusted grayscale image');>> %(3)>> se=strel('disk',10);>> topi=imtophat(gi,se);>> figure,imshow(topi),title('top-hat image');>> %(4)>> for counter=0:22 remain=imopen(clahei,strel('disk',counter));intensity_area(counter+1)=sum(remain(:));end >> figure,plot(intensity_area,'m-*'),grid on;>> title('sum of opening(pixels)');>> title('sum of opening values in opened image as a function of radius');>> xlabel('radius of opening(pixels)');>> ylabel('pixel value sum of opened objects(intensity)');>> >> >> >> for counter=0:20 remain=imopen(topi,strel('disk',counter));surfarea(counter+1)=sum(remain(:));end >> figure,plot(surfarea,'m-*'),grid on;>> set(gca,'xtick',[0 2 4 6 8 10 12 14 16 18 20]);>> title('surface area of opened objects as a function of radius');>> xlabel('radius of opening(pixels)');>> ylabel('surface area of opened objects(pixels)');>> %(5)>> intensity_area_prime=diff(intensity_area);>> figure,plot(intensity_area_prime,'m-*'),grid on;>> title('Granulometry(size distrubution)of snowflakes');>> set(gca,'xtick',[0 2 4 6 8 10 12 14 16 18 20 22]);>> xlabel('radius of snowflakes(pixels)');>> ylabel('sum of pixel values in snowflakes as a function of radius');>> derivsurfarea=diff(surfarea);>> figure,plot(derivsurfarea,'m-*'),grid on;>> title('granulometry(size distribution)of stars');>> xlabel('radius of stars(pixels)');>> ylabel('loss of pixels between two successive openings');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.2 花椒检测 clc;clear;close all;i=imread('gj.jpg');imshow(i);icanny=edge(i,'canny');imshow(icanny);se90=strel('line',2,90);se0=strel('line',2,0);bwsdil=imdilate(icanny,[se90 se0]);figure,imshow(bwsdil),title('dilated');ifill=imfill(bwsdil,'holes');figure,imshow(ifill);
%bwero=imerode(bwsdil,[se90 se0]);%figure,imshow(bwero);%i2fill=imfill(bwero,'holes');%figure,imshow(bwero);%imshow(i2fill);
%bwnobord=imclearborder(bwsdil,4);%figure,imshow(bwnobord);bwnobord=imclearborder(ifill,4);figure,imshow(bwnobord);se=strel('disk',5);bwc=imclose(bwnobord,se);bwco=imopen(bwnobord,se);figure,imshow(bwc);figure,imshow(bwco);%mask=bwsdil&bwco;%figure,imshow(mask);clc [labeled,numobjects]=bwlabel(bwco);numobjects
jdata=regionprops(labeled,'all');%jdata
jarea=[jdata.Area];mean(jarea)max(jarea)min(jarea)hist(jarea,255)jdata.Eccentricity %std([jdata.Eccentricity]）/(Mean([jdata.Eccentricity])jstd=std([jdata.Eccentricity])jmean=Mean([jdata.Eccentricity])jcv=jstd/jmean
>> std([jdata.Area])/ mean([jdata.Area])%面积的变异系数
>> std([jdata.Eccentricity])/ mean([jdata.Eccentricity])%椭圆的变异系数 >> std([jdata.MajorAxisLength])/ mean([jdata.MajorAxisLength])>> std([jdata.MinorAxisLength])/ mean([jdata.MinorAxisLength])>> std([jdata.EquivDiameter])/ mean([jdata.EquivDiameter])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.06.06 rice.png close all;clear;clc >> i=imread('rice.png');imshow(i);background=imopen(i,strel('disk',15));figure,imshow(background);i2=imsubtract(i,background);figure,imshow(i2);i3=imadjust(i2,stretchlim(i2),[0 1]);figure,imshow(i3);level=graythresh(i3);bw=im2bw(i3,level);figure,imshow(bw);imnobord=imclearborder(bw);[label,numobjects]=bwlabel(imnobord,4);numobjects rgb_label=label2rgb(label,@spring,'c','shuffle');figure,imshow(rgb_label);graindata=regionprops(label,'all');graindata
>> numdown=find([graindata.Area]<150);>> graindata(numdown,:)=[];>> numup=find([graindata.Area]>250);>> graindata(numup,:)=[];>> graindata
>> std([graindata.Area])/ mean([graindata.Area])%面积的变异系数
>> std([graindata.Eccentricity])/ mean([graindata.Eccentricity])%椭圆的变异系数
>> std([graindata.MajorAxisLength])/ mean([graindata.MajorAxisLength])>> std([graindata.MinorAxisLength])/ mean([graindata.MinorAxisLength])>> std([graindata.EquivDiameter])/ mean([graindata.EquivDiameter])%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.06.06 rice的垩白度检测 >> clear;close all;clc;>> rgb=imread('r.jpg');>> close all;>> imshow(rgb);>> i=rgb2gray(rgb);>> j=medfilt2(i,[5 5]);>> figure,imshow(i);>> figure,imshow(j);>> imhist(j,256);>> t=0.3;>> v=imadjust(j,[t 1],[],1);>> imhist(v,256);>> t_c=0.6;>> bw_v=im2bw(v,0.01);>> chalk=imadjust(v,[t_c 1],[],1);>> bw_chalk=im2bw(chalk,0.01);>> figure,imshow(v);>> figure,imshow(bw_v);>> figure,imshow(chalk);>> figure,imshow(bw_chalk);>> degree_chalkness=bwarea(bw_chalk)/bwarea(bw_v)*100 >> bw=im2bw(j,t);>> figure,imshow(bw);>> se=(ones(3,3));>> bw1=imerode(bw,se);%两次腐蚀 >> figure,imshow(bw1);>> bw2=imerode(bw1,se);>> figure,imshow(bw2);
>> [l,num]=bwlabel(bw2);%标记腐蚀后的大米图像 >> t_chalk=100;%设置垩白面积的下限 >> compare=(l)&(chalk>t_chalk);%>> compare=(bw2)&(bw_chalk>t_chalk);>> [r,c]=find(compare);%标记垩白米粒的位置 >> result=bwselect(l,c,r);%显示只含有垩白米粒的图像 >> figure,imshow(result);
>> [l_chalk,num_chalk]=bwlabel(result);%标记垩白米粒图像，便于计数 >> rate_chalky_grains=num_chalk/num*100;>> rate_chalky_grains
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.17 bwmorph函数 >> help bwmorph BWMORPH Perform morphological operations on binary image.BW2 = BWMORPH(BW1,OPERATION)applies a specific morphological operation to the binary image BW1.BW2 = BWMORPH(BW1,OPERATION,N)applies the operation N times.N can be Inf, in which case the operation is repeated until the image no longer changes.OPERATION is a string that can have one of these values: 'bothat' Subtract the input image from its closing 'bridge' Bridge previously unconnected pixels 'clean' Remove isolated pixels(1's surrounded by 0's)'close' Perform binary closure(dilation followed by erosion)'diag' Diagonal fill to eliminate 8-connectivity of background 'dilate' Perform dilation using the structuring element ones(3)'erode' Perform erosion using the structuring element ones(3)'fill' Fill isolated interior pixels(0's surrounded by 1's)'hbreak' Remove H-connected pixels 'majority' Set a pixel to 1 if five or more pixels in its 3-by-3 neighborhood are 1's 'open' Perform binary opening(erosion followed by dilation)'remove' Set a pixel to 0 if its 4-connected neighbors are all 1's, thus leaving only boundary pixels 'shrink' With N = Inf, shrink objects to points;shrink objects with holes to connected rings 'skel' With N = Inf, remove pixels on the boundaries of objects without allowing objects to break apart 'spur' Remove end points of lines without removing small objects completely.'thicken' With N = Inf, thicken objects by adding pixels to the exterior of objects without connected previously unconnected objects 'thin' With N = Inf, remove pixels so that an object without holes shrinks to a minimally connected stroke, and an object with holes shrinks to a ring halfway between the hold and outer boundary 'tophat' Subtract the opening from the input image
Class Support-------------The input image BW1 can be numeric or logical.It must be 2-D, real and nonsparse.The output image BW2 is logical.Examples--------BW1 = imread('circles.png');imview(BW1)BW2 = bwmorph(BW1,'remove');BW3 = bwmorph(BW1,'skel',Inf);imview(BW2)imview(BW3)
See also erode, dilate, bweuler, bwperim.Reference page in Help browser doc bwmorph
BW1 = imread('circles.png');figure,imshow(BW1)BW2 = bwmorph(BW1,'erode');figure,imshow(BW2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %边界提取 b=bwmorph(bw,'remove');b=bwperim(bw,8);%又叫边界象素测定 b=edge(bw,'canny');%又叫边界提取 %去除孤立象素点
nosinglepixel=bwmorph(bw,'clean');%去除小面积物体
nosmall=bwareaopen(bw,CNN);%阈值处理再取反
bw=~im2bw(i,graythresh(i));
%开运算（消除小物体）与闭运算（填充物体内细小空洞）se=strel('disk',6);iopen=imopen(bw,se);iclose=imclose(bw,se);%腐蚀与膨胀联合操作 %（1）创建结构元素 se=strel('rectangle',[40 30]);%（2）使用结构元素腐蚀图像 bw1=imread('circbw.tif');bw2=imerode(bw1,se);imshow(bw2);%（3）逆操作，回复矩形原来大小 bw3=imdilate(bw2,se);figure,imshow(bw3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.18花椒子
%直接对灰度图进行canny运算 >> i=imread('nut.bmp');>> figure,imshow(i);>> ig=rgb2gray(i);>> figure,imshow(ig);%igcanny=edge(ig,'canny');%igcfill=imfill(igcanny,'hole');igcanny_thresh=edge(ig,'canny',(graythresh(ig)*.1));igcfill=imfill(igcanny_thresh,'hole');>> figure,imshow(igcfill);
%先对灰度图滤波，再进行canny运算
>> imed=medfilt2(ig);%中值滤波后对图像边界有一定的损伤！！>> imedcanny=edge(imed,'canny');>> imedfill=imfill(imedcanny,'hole');>> figure,imshow(imedfill);>> nosmall=bwareaopen(imedfill,150);>> figure,imshow(nosmall);
%注意：若对灰度图像先拉氏锐化，在canny提取边界，效果不大好！！%结论：无需拉氏锐化，也不必中值滤波，可直接canny提取边界！！>> ifill=igcfill|imedfill;>> figure,imshow(ifill);>> nosmall=bwareaopen(ifill,150);>> figure,imshow(nosmall);
%当t=0.55时,阈值处理再canny运算的效果 >> imhist(ig);>> t=0.55;>> v=imadjust(ig,[0 t],[],1);>> vcanny=edge(v,'canny');>> vfill=imfill(vcanny,'hole');>> figure,imshow(vfill);>> ifill=igcfill|vfill;>> figure,imshow(ifill);>> nosmall=bwareaopen(ifill,150);>> figure,imshow(nosmall);
%当t=0.6时,阈值处理再canny运算的效果的效果 >> t=0.6;>> v=imadjust(ig,[0 t],[],1);>> vcanny=edge(v,'canny');>> vfill=imfill(vcanny,'hole');>> figure,imshow(vfill);>> ifill=igcfill|vfill;>> figure,imshow(ifill);>> nosmall=bwareaopen(ifill,150);>> figure,imshow(nosmall);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %处理花椒子
>> i=imread('nut.bmp');%figure,imshow(i);ig=rgb2gray(i);figure,imshow(ig);>> imed=medfilt2(ig);imedcanny=edge(imed,'canny');imedfill=imfill(imedcanny,'hole');%figure,imshow(imedfill);nosmall=bwareaopen(imedfill,150);>> figure,imshow(nosmall);>> [labeled,numobjects]=bwlabel(nosmall,4);>> rgb_label=label2rgb(labeled,@spring,'c','shuffle');%>> figure,imshow(rgb_label);>> nutdata=regionprops(labeled,'all');>> min([nutdata.Solidity])
>> rectangle('Position', [253.5000 207.5000 26.0000 28.0000])%画矩形
>> rectangle('Position', [250.5000 50.5000 27.0000 26.0000])>> figure,imshow(nutdata(1).Image)%只显示1号物体的图像
>> figure,imshow(nutdata(1).ConvexImage)%画出1号物体的凸多边形 >> std([nutdata.Eccentricity])/ mean([nutdata.Eccentricity])std([nutdata.Area])/ mean([nutdata.Area])std([nutdata.Solidity])/ mean([nutdata.Solidity])>> std([nutdata.Centroid])/ mean([nutdata.Centroid])std([nutdata.MajorAxisLength])/ mean([nutdata.MajorAxisLength])std([nutdata.MinorAxisLength])/ mean([nutdata.MinorAxisLength])std([nutdata.Orientation])/ mean([nutdata.Orientation])std([nutdata.EquivDiameter])/ mean([nutdata.EquivDiameter])std([nutdata.Extent])/ mean([nutdata.Extent])std([nutdata.Extrema])/ mean([nutdata.Extrema])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %处理花椒皮 close all;clc;clear;>> i=imread('p.bmp');imshow(i);ig=rgb2gray(i);figure,imshow(ig);imed=medfilt2(ig);imedcanny=edge(imed,'canny');figure,imshow(imedcanny);>> se90=strel('line',2,90);se0=strel('line',2,0);bwsdil=imdilate(imedcanny,[se90 se0]);figure,imshow(bwsdil),title('dilated');ifill=imfill(bwsdil,'holes');figure,imshow(ifill);>> bwero=imerode(ifill,[se90 se0]);>> figure,imshow(bwero);>> nosmall=bwareaopen(bwero,150,4);>> figure,imshow(nosmall);>> nobord=imclearborder(nosmall,4);>> figure,imshow(nobord);>> [labeled,numobjects]=bwlabel(nobord,4);>> numobjects >> pdata=regionprops(labeled,'all');>> max([pdata.Solidity])>> std([pdata.Solidity])/mean([pdata.Solidity])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %处理混合图像 >> clear;clc;close all;>> i=imread('m.bmp');%>> figure,imshow(i);>> ig=rgb2gray(i);figure,imshow(ig);imed=medfilt2(ig);%>> figure,imshow(imed);imedcanny=edge(imed,'canny');%>> figure,imshow(imedcanny);>> se90=strel('line',2,90);se0=strel('line',2,0);bwsdil=imdilate(imedcanny,[se90 se0]);%figure,imshow(bwsdil),title('dilated');ifill=imfill(bwsdil,'holes');%figure,imshow(ifill);bwero=imerode(ifill,[se90 se0]);%figure,imshow(bwero);>> nosmall=bwareaopen(bwero,150,4);%figure,imshow(nosmall);nobord=imclearborder(nosmall,4);figure,imshow(nobord);>> [labeled,numobjects]=bwlabel(nobord,4);>> numobjects >> rgb_label=label2rgb(labeled,@spring,'c','shuffle');figure,imshow(rgb_label);>> mexdata=regionprops(labeled,'all');hold on;%以下内容画在同一figure中 centr=[mexdata.Centroid];%寻找重心位置 nums=1:numobjects;for k = 1:numobjects soli=mexdata(k).Solidity;soli_string=sprintf('%2.2f',soli);%等价于转字符串 % signal=num2str(nums(k));signal=sprintf('%d',k);%直接使用打印语句打印序号 text(centr(2*k-1),centr(2*k),signal)%按序标记物体
text(centr(2*k-1)-30,centr(2*k)-30,soli_string)%标注每个Solidity值 end
for k=1:numobjects plot(mexdata(k).ConvexHull(:,1),mexdata(k).ConvexHull(:,2),...'b','Linewidth',2)end
%画出1和2号物体的外接矩形
%>> rectangle('position',[9.5000 224.5000 62.0000 63.0000])%>> rectangle('position',[65.5000 141.5000 34.0000 39.0000])%画出每个物体的外接矩形 bb=[mexdata.BoundingBox];for k=1:numobjects rectangle('position',[bb(4*k-3)bb(4*k-2)bb(4*k-1)bb(4*k)])end
%>> figure,imshow(mexdata(1).Image)%只显示1号物体的图像
%>> figure,imshow(mexdata(1).ConvexImage)%画出1号物体的凸多边形 %>> figure,imshow(mexdata(2).Image)%只显示2号物体的图像
%>> figure,imshow(mexdata(2).ConvexImage)%画出2号物体的凸多边形 %画出单个物体的凸多边形的填充图形 for k=1:numobjects figure,imshow(mexdata(k).ConvexImage)end
%只显示Solidity>0.92的物体的图像 >> idx = find([mexdata.Solidity] > 0.92);>> BW2 = ismember(labeled,idx);>> figure,imshow(BW2)
>> mexdata=regionprops(labeled,'all');>> %只显示Solidity<0.92的物体的图像 idx = find([mexdata.Solidity] < 0.92);bw2 = ismember(labeled,idx);figure,imshow(bw2)%mexdata.Solidity;
>> numdown=find([mexdata.Solidity]<0.92);mexdata(numdown,:)=[];>> mexdata
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.19 %roipoly函数的用法 I = imread('eight.tif');c = [222 272 300 270 221 194];r = [21 21 75 121 121 75];BW = roipoly(I,c,r);imview(I), imview(BW)
%可以使用下面的方法创建相应的向量： regionprops(L,'Area');allArea = [stats.Area];
%创建一个只包含面积大于80的二值图像 idx = find([stats.Area] > 80);BW2 = ismember(L,idx);
%只显示某个下标所对应的物体图像 bw2=ismember(L,N);figure,imshow(bw2);
%在调用regionprops之前必须将二值图像转变为标注矩阵 L = bwlabel(BW);%或者
L = double(BW);
%将matlab数据写到excel中 a=ones(3);success = xlswrite('c:/matlab/work/myworkbook.xls',a,'A2:C4')%将行矩阵转换为列矩阵 a=[1 2 3 4 5 6];b=transpose(a);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.22球形物体的检测和标识（循环检测和标识算法）clc;clear;close all;%Step 1: Read image %Step 2: Threshold the image %Step 3: Remove the noise %Step 4: Find the boundaries %Step 5: Determine which objects are round >> RGB = imread('pillsetc.png');imshow(RGB)>> I = rgb2gray(RGB);threshold = graythresh(I);bw = im2bw(I,threshold);imshow(bw)>> % remove all object containing fewer than 30 pixels bw = bwareaopen(bw,30);>> figure,imshow(bw)>> % fill a gap in the pen's cap se = strel('disk',2);bw = imclose(bw,se);>> figure,imshow(bw)>> % fill any holes, so that regionprops can be used to estimate % the area enclosed by each of the boundaries bw = imfill(bw,'holes');>> figure,imshow(bw)>> [B,L] = bwboundaries(bw,'noholes');>> % Display the label matrix and draw each boundary figure,imshow(label2rgb(L, @jet, [.5.5.5]))>> hold on for k = 1:length(B)boundary = B{k};plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)end >> stats = regionprops(L,'Area','Centroid');>> stats = regionprops(L,'Area','Centroid');threshold = 0.94;% loop over the boundaries for k = 1:length(B)% obtain(X,Y)boundary coordinates corresponding to label 'k' boundary = B{k};% compute a simple estimate of the object's perimeter delta_sq = diff(boundary).^2;perimeter = sum(sqrt(sum(delta_sq,2)));
% obtain the area calculation corresponding to label 'k' area = stats(k).Area;
% compute the roundness metric metric = 4*pi*area/perimeter^2;
% display the results metric_string = sprintf('%2.2f',metric);% mark objects above the threshold with a black circle if metric > threshold centroid = stats(k).Centroid;plot(centroid(1),centroid(2),'ko');end
text(boundary(1,2)-35,boundary(1,1)+13,metric_string,'Color','y',...'FontSize',14,'FontWeight','bold');end >> title(['Metrics closer to 1 indicate that ',...'the object is approximately round']);
第四篇：基于MATLAB图像处理报告
基于MATLAB图像处理报告
一、设计题目
图片叠加。
二、设计要求
将一幅礼花图片和一幅夜景图片做叠加运算，使达到烟花夜景的美图效果。
三、设计方案
3.1、设计思路
利用matlab强大的图像处理功能，通过编写程序，实现对两幅图片的像素进行线性运算，利用灰度变换的算法使图片达到预期的效果。
3.2、软件介绍
MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。可以直接调用，用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户直接进行下载就可以用。
3.3、常见简单程序语句及算法分析
（1）CLC；清零。（2）CLEAR ； 清内存。
（3）r=imread（‘路径图片名.jpg’）；读入一幅图片。（4）imshow（r）；显示图片r。（5）g=rgb2gray(r)；灰度转变。
（6）%imhist（g）；g1=histeq（g）；figure，imsho（g1）；直方图均衡化。（7）Imwrite（‘g1，路径图片名.jpg’）;保存图片。
3.4、图片叠加及灰度变换分析
程序1 [m,n,l]=size(C);for i=1:m for j=1:n for k=1:l D(i,j,k)=C(i,j,k)+ B(i,j,k);end end end 此程序的主要功能是对两幅图片通过算法来实现叠加的效果，程序中的几个变量都是像素点的值，通过三个循环使得两幅图片对应的值线性相加，最大值应该是以255输出，超过255也是以255输出。
程序2 J = imadjust(I,[low_in high_in], [low_out high_out])此程序变换的原理是：如果原图像f(x, y)的灰度范围是[m, M]，我们希望对图像的灰度范围进行线性调整，调整后的图像g(x, y)的灰度范围是[n, N]，那么下述变换：g(x,y)
N
n
f(x,y)
m
n就可以实现这一要求。MATLAB图像处理工具箱中提供的imadjust函数，可以实现上述的线性变换对比度调整。
四、设计步骤（1）处理之前，我们先来看看两幅原图，一幅是带有礼花的图片，另一幅是一幅东方明珠的夜景图。
图 1
图 2 通过图片我们发现，礼花的图片非常的亮，而夜景图则显得有些灰暗。我们推测，如果让礼花和夜景的图片叠加到一起，会不会由于礼花的亮度太大而掩盖了城市的夜光，达不到我们想要的效果。
（2）打开matlab图像处理软件，7.0及以上版本就可以，打开之后，显示界面如下（我的版本是2012b）。
图3（3）新建script文件，点击’New script‘，或点击’New‘，选择script。
图4（4）在打开的界面之中输入程序。图片的位置是你所要用的图片的准确位置，尽量写详细点，减少计算机的读取时间，加快读取速度。如图所示。
图5（5）点击运行按钮，即界面上的绿色按钮。
图6（6）查看效果，如下图。
图7 通过处理后的图片，我们看到由于礼花太亮，完全将城市的夜景掩盖住了，效果不是太理想。我们设想加入灰度变换函数，将礼花的图片变的暗一点，将夜景的亮度提高一点，看看效果怎么样，加入灰度变换程序，如下图。
图8（7）软件调试后运行程序，如下图。
图9 结果显示加入灰度变换的函数之后，图片效果较之前好多了，达到了预期的目的。
五、软件代码
通过matlab进行图像处理，代码如下： clear;clc;A=imread('C:UsersAdministratorDesktop作业礼花.jpg');C=imresize(A,[300,400]);B=imread('C:UsersAdministratorDesktop作业夜景提亮天空中加入礼花.jpg');[m,n,l]=size(C);for i=1:m for j=1:n for k=1:l D(i,j,k)=C(i,j,k)+ B(i,j,k);end end end figure imshow(D);imwrite('C:UsersAdministratorDesktop作业礼花效果图1.jpg')调试之后的程序代码如下： clear;clc;A=imread('C:UsersAdministratorDesktop作业礼花.jpg');B=imread('C:UsersAdministratorDesktop作业夜景提亮天空中加入礼花.jpg');A1=imadjust(A,[0,1],[0,0.9]);B1=imadjust(B,[0.3,0.6],[0,1]);C=imresize(A1,[300,400]);[m,n,l]=size(C);for i=1:m for j=1:n for k=1:l D(i,j,k)=C(i,j,k)+ B1(i,j,k);end end end figure imshow(D);imwrite('C:UsersAdministratorDesktop作业礼花效果图2.jpg')处理后的图片如下：
图 10
六、结果分析
通过两幅图的对比，发现第二幅图片较第一幅，效果明显增强。是由于加入灰度变换函数，使原图的灰度值发生变化，以达到实际的效果。
效果对比图11
七、心得体会
通过这次任务，以前在课堂上没太听明白怎么回事的东西，通过上网，查资料，以及用软件处理，通通实践了一遍，加深了对这门课程的认识和理解。Matlab是一款功能很强大的应用软件，它不仅可以对图像进行处理，而且可以进行各种数字计算和符号计算功能，具有绘图功能，语言体系等等。这次的任务我们乐在其中，喜悦的是看到了成果，内心充满了满足感和成就感，就如同看到了图片中的烟花，有种过年的感觉。不过这些都只是皮毛而已，要想真正地掌握它，还得更进一步地学习理论知识。
参考文献：
【冈萨雷斯 数字图像处理（MATLAB版）】[美] RafaelC.Gonzalez RichardE.Woods StevenL.Eddins 著 电子工业出版社。
【数字图像处理及MATLAB实现】杨杰 主编 电子工业出版社。
第五篇：MATLAB实现数字信号处理
数字信号处理
说 明 书
目录
一．摘要…………………………………3 二．课程设计目的………………………3 三．设计内容……………………………3 四．设计原理……………………………4 4.1．语音信号的采集…………………………….4 4.2．滤波器……………………………………….4 4．21．IIR滤波器原理…………………………………….4 4.22．FIR滤波器原理………………………………………5 五．设计步骤……………………………6 5.1录制女音………………………………………6 5.2采样语音信号并画出时域波形和频谱图……7 5.3采用双线性变换法设计IIR滤波器…………10 5.4窗函数法设计FFR滤波器………………......12 5.5用IIR滤波器对信号进行滤波………………14 5.6用FIR滤波器对信号进行滤波………………16 5.7男女声语音信号频谱特点分析………………19 5.8有背景噪声的信号分析………………………20 六．心得体会…………………………….22 七．参考文献…………………………….23
一．摘要：
这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。
MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。
在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。
二．课程设计目的：
综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。
三．设计内容：
内容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数 法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。
四．设计原理：
4.1．语音信号的采集
熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有：
a：y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。Dtype为采样数据的存储格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据；
b：wavplay(y,fs)；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y；
c：wavwrite((y,fs,wavfile);创建音频文件； d：y=wavread(file);读取音频文件；
关于声音的函数还有sound()；soundsc()；等。4.2滤波器： 4．21.IIR滤波器原理
冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。
双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也 就消除了频谱混叠现象。
双线性法设计IIR数字滤波器的步骤：
1）将数字滤波器的频率指标{ k}由Wk＝（2/T）*tan(wk),转换为模拟滤波器的频率指标{k}.2)由模拟滤波器的指标设计H(s).3)由H(s)转换为H(z)21z1H(z)H(s)sT1z1
4.22.FIR滤波器原理
FIR滤波器与IIR滤波器特点不同，IIR滤波器的相位是非线性的，若需线性相位则要采用全通网络进行相位校正。而有限长单位冲激响应（FIR）数字滤波器就可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。
由于FIR系统的冲激响应就是其系统函数各次项的系数，所以设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲激响应作为H（z）的系数，冲激响应长度N就是系统函数H（z）的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。这种时域设计、频域检验的方法一般要反复几个回合，不像IIR DF设计靠解析公式一次计算成功。给出的理想滤波器频率响应是，它是w的周期函数，周期
由傅立叶反变换导出，即
hd(n)1Hd(ejw)ejwndw2，再将hd(n)与窗函数，因此可展开成傅氏级数w(n)相乘就可以得到h(n)。、的计算可采用傅氏变换的现成公式和程序，窗函数也是现成的。但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算，整个设计可用计算机编程来做。
窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度，常用的几种窗函数有：
矩形窗
w(n)=RN(n)；
Hanning窗
；
Hamming窗
；
Blackmen窗
；
Kaiser窗。
式中Io(x)为零阶贝塞尔函数。
五．设计步骤：
5.1录制女音：
利用MATLAB中的函数录制声音。function nvyin()fs=11025;
%采样频率
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('
开始录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double');
%录制声音3秒
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('
录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('
播放录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs);
%播放录音
str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('
播放录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音');
%声音的存储
5.2采样语音信号并画出时域波形和频谱图
读取语音信号，画出其时域波形和频谱图，与截取后的语音信号的时域波形和频谱图比较，观察其变化。程序如下：
[x,fs,bits]=wavread('女音.wav');
%读取声音
N=length(x);
%计数读取信号的点数 t=(1:N)/fs;
%信号的时域采样点 f0=fs/N;
%采样间隔 n=1:N/2;
%取信号的一半 figure(1);subplot(2,2,1);
%把画图区域划分为2行2列，指定第一个图 plot(t, x);
%画出声音采样后的时域波形 title('原女音信号的时域波形');
%给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid;
%添加网格
y=fft(x);
%对信号做N点FFT变换 k=(n-1)*f0;
%频域采样点
subplot(2,2,3);
%把画图区域划分为2行2列，指定第三个图 plot(k,abs(y(n)));
%绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性');
%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;
%添加网格
subplot(2,2,4);
%把画图区域划分为2行2列，指定第四个图 if y~=0
%判断指数是否为0
plot(k,20*log10(abs(y(n))));
%画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid;
%添加网格
%实际发出声音落后录制动作半拍的现象的解决 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]);
%截取语音信号 N=length(x1);
%计数读取信号的点数 t=(1:N)/fs;
%信号的时域采样点 f0=fs/N;
%采样间隔 n=1:N/2;
%取信号的一半
figure(2);subplot(2,2,1);
%把画图区域划分为2行2列，指定第一个图 plot(t,x1);
%画出声音采样后的时域波形 title('截取后女音信号的时域波形');
%给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid;
%添加网格
y1=fft(x1);
%对信号做N点FFT变换
subplot(2,2,3);
%把画图区域划分为2行2列，指定第三个图 k=(n-1)*f0;
%频域采样点
plot(k,abs(y(n)));
%绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性');
%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;
%添加网格
subplot(2,2,4);
%把画图区域划分为1行2列，指定第二个图 if y1~=0
%判断指数是否为0
plot(k,20*log10(abs(y1(n))));
%画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid;
%添加网格
原女音信号的时域波形10.5A/幅0振-0.5-10123时间/tFFT变换后声音的频谱特性FFT变换后声音的频谱特性30050A200贝/值分/幅0幅100振00200040006000-***频率/HzHz 截取后女音信号的时域波形10.5振幅/A0-0.5-10123FFT变换后声音的频谱特性50时间/tFFT变换后声音的频谱特性300200振幅/分贝幅值/A01000020004000频率/Hz6000-5002000Hz40006000
结果分析：
由原女音信号的时域波形可知录取开始时实际发出声音大概落后3500个采样点，我们把前3500点去除即可解决实际发出声音落后录制动作半拍的现象。由原女音的的频谱图和截取后声音的频谱图可看出，对声音的截取并不会影响它们频谱分布。
5.3采用双线性变换法设计IIR滤波器：
人的声音频率一般在（1~~4）kHZ之间，则我们只需要设计一个带通滤波器即可滤去声音频带以外的无用噪声，得到比较清晰的声音。根据声音的频谱图分析，设计一个带通滤波器性能指标如下：
fp1＝1000 Hz，fp2＝3000 Hz，fsc1＝500 Hz，fsc2＝3500Hz，As＝100dB，Ap＝1dB，fs=10000 程序如下：
%iir带通的代码： %w=2*pi*f/fs Ap=1;
%通带波纹系数
Az=100;
%最小阻带衰减
wp=[0.2 0.6];
%归一化通带数字截止频率 wz=[0.1 0.7];
%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);
%估计契比雪夫I型滤波器阶数 [b,a]=cheby1(N,Ap,wn);
%N指定滤波器阶数，wn归一化
截 %止频率，Ap通带波动
[h,w]=freqz(b,a);
%求数字滤波器的复频率响应 figure(1);subplot(2,1,1);plot(w/pi,abs(h));
%绘制数字滤波器的频谱图 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型带通滤波器的幅频响应');subplot(2,1,2);if abs(h)~=0
%判断指数是否为0
plot(w/pi,20*log10(abs(h)));
%绘制数字滤波器频谱的分贝图 end grid;xlabel('omega/pi');ylabel('振幅(分贝)');title('契比雪夫Ⅰ型带通滤波器的幅频响应');契比雪夫Ⅰ型带通滤波器的幅频响应1振幅（幅值）0.5000.10.20.50.60.70.8/契比雪夫Ⅰ型带通滤波器的幅频响应0.30.40.910振幅(分贝)-200-400-60000.10.20.30.40.5/0.60.70.80.91
5.4窗函数法设计FFR滤波器
线性相位FIR滤波器通常采用窗函数法设计。窗函数法设计FIR滤波器的基本思想是：根据给定的滤波器技术指标，选择滤波器长度N和窗函数ω（n），使其具有最窄宽度的主瓣和最小的旁瓣。其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h(n)。工程中常用的窗函数共有6种，即矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、汉明（Hamming）窗、布莱克曼（Blackman）窗和凯泽（Kaiser）窗。
这次设计我采用的是布莱克曼来设计给定数字带通滤波器的参数如下： wp1=0.3pi, wp2=0.6pi, wz1=0.2pi, wz2=0.7pi, Ap=1dB, Az=70dB 程序如下：
Ap=1;
%通带波纹系数 Az=100;
%最小阻带衰减 fs=10000;
%采样频率 wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));
%---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW);
%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);
%---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%理想带通滤波器 wdWindow=blackman(N);
%布拉克曼窗 hr=wdWindow.*hdWindow';
%点乘
n=0:N-1;
%阶数 subplot(2,2,1);stem(n,wdWindow);
%绘制布拉克曼窗时域波形 xlabel('时间');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1);
%求滤波器频率响应 subplot(2,2,3);plot(W/pi,abs(H))
%绘制滤波器频域波形 xlabel('omega/pi');ylabel('振幅');title('FIR带通滤波器幅频特性');subplot(2,2,4);
if abs(H)~=0
%判断指数是否为0
plot(W/pi,20*log10(abs(H)));
%画滤波器频谱的分贝图 end xlabel('omega/');ylabel('振幅/分贝');title('FIR带通滤波器幅频特性');grid;
%添加网格 %---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录ideallp.m function hd = ideallp(wc,N);% 理想低通滤波器的脉冲响应子程序 % hd = 点0 到 N-1之间的理想脉冲响应 % wc = 截止频率(弧度)% N = 理想滤波器的长度
tao =(N-1)/2;
% 理想脉冲响应的对称中心位置 n = [0:(N-1)];
% 设定脉冲响应长度 m = n-tao + eps;
% 加一个小数以避免零作除数
hd = sin(wc*m)./(pi*m);
% 理想脉冲响应
布拉克曼窗1振幅0.500406080时间FIR带通滤波器幅频特性500振幅/分贝20FIR带通滤波器幅频特性1.51振幅-50-100-15000.5/10.5000.5/1
5.5用IIR滤波器对信号进行滤波
用自己设计的IIR滤波器分别对采集的信号进行滤波，在Matlab中，IIR滤波器利用函数filter对信号进行滤波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x);
%计数读取信号的点数 t=(1:N)/fs;
%信号的时域采样点 f0=fs/N;
%采样间隔 n=1:N/2;
%取信号的一半 y=fft(x);
%对信号做N点FFT变换 k=(n-1)*f0;
%频域采样点
subplot(2,1,1);
%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));
%绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性');
%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;
%iir带通的代码：
Ap=1;
%通带波纹系数
Az=100;
%最小阻带衰减
wp=[0.2 0.6];
%归一化通带数字截止频率 wz=[0.1 0.7];
%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);
%估计契比雪夫I型滤波器阶数
[b,a]=cheby1(N,Ap,wn);
%N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x);
%对声音滤波 wavplay(x1)wavwrite(x1,'IIR滤波后女音.wav');N=length(x1);
%计数读取信号的点数 t=(1:N)/fs;
%信号的时域采样点 f0=fs/N;
%采样间隔 n=1:N/2;
%取信号的一半
y=fft(x1);
%对信号做N点FFT变换 k=(n-1)*f0;
%频域采样点
subplot(2,1,2);
%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));
%绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性');
%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;
滤波前女音的频谱特性300幅值/A***030004000频率/Hz滤波后女音的频谱特性500060006040幅值/A***0频率/Hz400050006000
结果分析：
由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3400Hz的频谱成分。回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。
5.6用FIR滤波器对信号进行滤波
用自己设计的FIR滤波器分别对采集的信号进行滤波，在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波 程序如下：
[x,fs,bits]=wavread('女音.wav');N=length(x);
%计数读取信号的点数
t=(1:N)/fs;
%信号的时域采样点 f0=fs/N;
%采样间隔 n=1:N/2;
%取信号的一半
y=fft(x);
%对信号做N点FFT变换 k=(n-1)*f0;
%频域采样点
subplot(2,1,1);
%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));
%绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性');
%给图形加注标签说明 xlabel('频率/omega');ylabel('幅值/A');grid;
%FIR带通滤波器代码 fs=10000;wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));
%---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW);
%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);
%---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);wdWindow=blackman(N);hr=wdWindow.*hdWindow';x1=fftfilt(hr,x);
%对声音滤波 wavplay(x1)wavwrite(x1,'FIR滤波后女音.wav');N=length(x1);
%计数读取信号的点数 t=(1:N)/fs;
%信号的时域采样点 f0=fs/N;
%采样间隔 n=1:N/2;
%取信号的一半
y=fft(x1);
%对信号做N点FFT变换 k=(n-1)*f0;
%频域采样点
subplot(2,1,2);
%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));
%绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性');
%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;
滤波前女音的频谱特性300200幅值/A***004000频率/l滤波后女音的频谱特性500060006040幅值/A20005001000***03000频率/Hz***0
结果分析：
由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3500Hz的频谱成分。和用IIR滤波器滤波一样，回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。5.7男女声语音信号频谱特点分析
换一个男音录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点。程序如下：
[x,fs,bits]=wavread('女音.wav');
%读取声音
N=length(x);
%计数读取信号的点数 t=(1:N)/fs;
f0=fs/N;
n=1:N/2;
y=fft(x);
k=(n-1)*f0;
subplot(2,1,1);
plot(k,abs(y(n)));
title('FFT变换后女音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid;
[x,fs,bits]=wavread('明明.wav');
N=length(x);
t=(1:N)/fs;
f0=fs/N;
n=1:N/2;
y=fft(x);
k=(n-1)*f0;
subplot(2,1,2);
plot(k,abs(y(n)));
title('FFT变换后男音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid;
%信号的时域采样点
%采样间隔
%取信号的一半
%对信号做N点FFT变换
%频域采样点
%把画图区域划分为2行1列，指定第一个图%绘制原始语音信号的幅频响应图
%给图形加注标签说明
%添加网格
%读取声音
%计数读取信号的点数
%信号的时域采样点
%采样间隔
%取信号的一半
%对信号做N点FFT变换
%频域采样点
%把画图区域划分为2行1列，指定第二个图%绘制原始语音信号的幅频响应图
%给图形加注标签说明
%添加网格
axis([0 6000 0 300])；
%改变横纵坐标便于比较频谱图
FFT变换后女音的频谱特性300200幅值/A***00频率/FFT变换后男音的频谱特性***200幅值/A***00频率/400050006000
结果分析：
通过比较上面女音频谱图和男音频谱图可知，男音的频谱集中在低频部分，高频成分底，谱线较平滑，声音听起来低沉。5.8有背景噪声的信号分析
从硬盘中把一段噪声（频谱能量集中在某个小范围内）叠加到语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除； 程序如下：
z=wavread('女音.wav',[1 24000]);
%读取声音在1-24000之间 f=wavread('noise.wav',[1 24000]);x=z+f;wavplay(x);fs=11025;N=length(x);f0=fs/N;
%采样间隔
n=1:N;
%取信号的一半 y=fft(x,N);%对信号做N点FFT变换
k=(n-1)*f0;
%频域采样点
subplot(2,1,1);
%把画图区域划分为1行2列，指定第二个图 plot(k,abs(y(n)));
%绘制原始语音信号的幅频响应图 title('加噪声后声音的频谱特性');
%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格
%iir带通滤波器的代码：
Ap=1;
%通带波纹系数
Az=70;
%最小阻带衰减
wp=[0.2 0.7];
%归一化通带数字截止频率 wz=[0.1 0.8];
%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);
%估计契比雪夫I型滤波器阶数
[b,a]=cheby1(N,Ap,wn);
%N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x);
%对声音滤波 wavplay(x1)；
wavwrite(x1,'滤除噪音后女音.wav');N=length(x1);f0=fs/N;
%采样间隔 n=1:N;
%取信号的一半
y1=fft(x1,N);
%对信号做fs点FFT变换
subplot(2,1,2);
%把画图区域划分为1行2列，指定第二个图 k=(n-1)*f0;
%频域采样点
plot(k,abs(y1(n)));
%绘制原始语音信号的幅频响应图 title('滤除噪声后声音的频谱特性');
%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格
加噪声后声音的频谱特性3000幅值/A***0008000频率/Hz滤除噪声后声音的频谱特性***030幅值/A***000频率/Hz80001000012000
结果分析
观察加噪声后声音的频谱图可知，噪音频率主要在4000Hz处，只要我们设计一个，滤波器滤除大概在4000Hz的频谱即可，回放滤波后的语音信号，可证噪音基本滤除。
六．心得体会：
通过这次课程设计，让我对MATLAB的基本应用有了更深的了解，还有数字信号处理在MATLAB中的一些函数的用法。通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。
在这次实验中，也遇到了很多问题，比如画信号频谱的分贝图时（20*log10(abs(y))）指数为零时的处理。滤波器的设计也花了好大的功夫，刚开始不会设计参数，一头雾水，通过同学的指导和讨论，得知通过观察信号的频谱图，看噪音频率集中在那一部分，设计滤波器把其滤除即可。可反复设置参数直到滤波后语音信号的效果好为止。
七．参考文献：
（1）《MATLAB LabVIEW SystemView》翁剑枫 叶志前 编著, 机械工业出版社；
（2）《MATLAB及在电子信息课程中的应用》陈怀琛 吴大正 高西全编著，电子工业出版社；
（3）《MATLAB在数字信号处理中的应用》（弟2版）薛年喜 编著，清华大学出版社；
（4）《MATLAB扩展编程》何强 何英
编著，清华大学出版社；（5）《MATLAB7简明教程》吴清 曹辉林 编著，清华大学出版社；（6）MATLAB5.3精要.编程及高级应用》程卫国 冯峰 王雪梅 刘艺 编著，机械工程出版社。
点击咨询 