第一篇：MATLAB函數(shù)處理圖像實現(xiàn)膨脹腐蝕

MATLAB函數(shù)處理圖像實現(xiàn)膨脹腐蝕

一、實驗?zāi)康?/p>

1、了解二值形態(tài)學的基本運算

2、掌握二值圖像膨脹、腐蝕的基本方法

3、編程實現(xiàn)膨脹、腐蝕

二、實驗要求

1、使用imdilate函數(shù)進行圖像膨脹，并觀察膨脹后圖像的變化。

2、使用imerode函數(shù)進行圖像腐蝕，觀察腐蝕后的圖像變化情況。

三、實驗原理

膨脹：將與物體接觸的所有背景點合并到該物體中，使邊界向外部擴張的過程。利用它可以填補物體中的空洞。B對X膨脹所產(chǎn)生的二值圖像D是滿足以下條件的點（x,y）的集合:如果B的原點平移到點（x,y），那么它與X的交集非空。數(shù)學表達式：C?A?B

腐蝕：一種消除邊界點，使邊界向內(nèi)部收縮的過程。利用它可以消除小而且無意義的物體。B對X腐蝕所產(chǎn)生的二值圖像E是滿足以下條件的點（x,y)的集合：如果B的原點平移到點（x,y），那么B將完全包含于X中。數(shù)學表達式：C?A?B

膨脹處理：一種消除邊界點，使邊界點向內(nèi)部收縮的過程。

腐蝕處理：將與物體接觸的所有背景點合并到該物體中，使邊界向外部擴張的過程。

四、實驗步驟

1.圖像膨脹的Matlab實現(xiàn)：

可以使用imdilate函數(shù)進行圖像膨脹,imdilate函數(shù)需要兩個基本輸入?yún)?shù)，即待處理的輸入圖像和結(jié)構(gòu)元素對象。結(jié)構(gòu)元素對象可以是strel函數(shù)返回的對象，也可以是一個自己定義的表示結(jié)構(gòu)元素鄰域的二進制矩陣。此外，imdilate還可以接受兩個可選參數(shù)：PADOPT(padopt)——影響輸出圖片的大小、PACKOPT(packopt).——說明輸入圖像是否為打包的二值圖像(二進制圖像)。步驟1，首先創(chuàng)建一個包含矩形對象的二值圖像矩陣。>> 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的正方形結(jié)構(gòu)元素對象對創(chuàng)建的圖像進行膨脹。>> SE=strel('square',3)SE = Flat STREL object containing 9 neighbors.Neighborhood: 1 1 1 1 1 1 1 1 1 步驟3，將圖像BW和結(jié)構(gòu)元素SE傳遞給imdilate函數(shù)。>> 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，顯示結(jié)果。>> imshow(BW,'notruesize')>> imshow(BW2,'notruesize')2.圖像腐蝕的Matlab實現(xiàn)：

可以使用imerode函數(shù)進行圖像腐蝕。imerode函數(shù)需要兩個基本輸入?yún)?shù)：待處理的輸入圖像以及結(jié)構(gòu)元素對象。此外，imerode函數(shù)還可以接受3個可選參數(shù)：PADOPT(padopt)——影響輸出圖片的大小、PACKOPT(packopt).——說明輸入圖像是否為打包的二值圖像(二進制圖像)。M——指定原始圖像的行數(shù)。以下程序示例說明了如何對某一副具體圖像進行腐蝕操作，腐蝕前后的效果對比如圖末。

步驟1，讀取圖像cameraman.tif（該圖像是Matlab當前目錄下自帶的圖片）>> BW1=imread('cameraman.tif');步驟2，創(chuàng)建一個任意形狀的結(jié)構(gòu)元素對象 >> SE=strel('arbitrary',eye(5));步驟3，以圖像BW1和結(jié)構(gòu)元素SE為參數(shù)調(diào)用imerode函數(shù)進行腐蝕操作。>> BW2=imerode(BW1,SE);步驟4，顯示操作結(jié)果 >> imshow(BW1)>> figure,imshow(BW2)

五、實驗代碼及結(jié)果

代碼：

imerode函數(shù)，該函數(shù)能夠?qū)崿F(xiàn)二值圖像的腐蝕操作； imdilate函數(shù)，該函數(shù)能夠?qū)崿F(xiàn)二值圖像的膨脹操作； 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(‘腐蝕后的圖像’)

結(jié)果：

原

圖

膨脹后的圖像

腐蝕后的圖像

六、實驗心得體會

通過本次的實驗，我了解了二值形態(tài)學的基本運算，掌握了二值圖像膨脹、腐蝕的基本方法，并且會運用編程實現(xiàn)膨脹、腐蝕，本次的實驗?zāi)康囊呀?jīng)完成，意識到在以后的生活中要了解做事情的目的，注重每一個與細節(jié)，認真思考遇到的所有問題，提高自己各方面的能力。感謝尹強老師教會我們理論與實踐知識，也讓我明白了什么是學習，怎么樣學習，為以后的生活奠定的基礎(chǔ)與指引了方向。

第二篇：matlab數(shù)字圖像處理 膨脹和腐蝕

基于Matlab的腐蝕和膨脹的邊緣檢測

一、實驗?zāi)康? 掌握運用Matlab軟件對灰度與二值圖像的膨脹與腐蝕的處理方法。

二、實驗環(huán)境（軟件條件）：

Windws2000/XP

MATLAB 7.x

三、實驗內(nèi)容：

1、圖像膨脹的Matlab實現(xiàn) ① 實驗原理：

膨脹：給圖像中的對象邊界添加像素。

在操作中，輸出圖像中所有給定像素的狀態(tài)都是通過對輸入圖像的相應(yīng)像素及鄰域使用一定的規(guī)則進行確定。在膨脹操作時，輸出像素值是輸入圖像相應(yīng)像素鄰域內(nèi)所有像素的最大值。在二進制圖像中，如果任何像素值為1，那么對應(yīng)的輸出像素值為1。

可以使用imdilate函數(shù)進行圖像膨脹,imdilate函數(shù)需要兩個基本輸入?yún)?shù)，即待處理的輸入圖像和結(jié)構(gòu)元素對象。結(jié)構(gòu)元素對象可以是strel函數(shù)返回的對象，也可以是一個自己定義的表示結(jié)構(gòu)元素鄰域的二進制矩陣。此外，imdilate還可以接受兩個可選參數(shù)：PADOPT(padopt)——影響輸出圖片的大小、PACKOPT(packopt).——說明輸入圖像是否為打包的二值圖像(二進制圖像)。② 實驗步驟：

A、首先創(chuàng)建一個包含矩形對象的二值圖像矩陣。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的正方形結(jié)構(gòu)元素對象對創(chuàng)建的圖像進行膨脹。

C=strel('square',3)C = Flat STREL object containing 9 neighbors.Neighborhood: 1 1 1 1 1 1 1 1 1 C、將圖像R和結(jié)構(gòu)元素C傳遞給imdilate函數(shù)。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、顯示結(jié)果： imshow(R,'notruesize')imshow(R1,'notruesize')③ 實驗結(jié)果： 膨脹前后效果圖：

2、圖像腐蝕的Matlab實現(xiàn) ① 實驗原理：

腐蝕：刪除對象邊界某些像素。

在操作中，輸出圖像中所有給定像素的狀態(tài)都是通過對輸入圖像的相應(yīng)像素及鄰域使用一定的規(guī)則進行確定。在腐蝕操作中，輸出像素值是輸入圖像相應(yīng)像素鄰域內(nèi)所有像素的最小值。在二進制圖像中，如果任何一個像素值為0，那么對應(yīng)的輸出像素值為0。

可以使用imerode函數(shù)進行圖像腐蝕。imerode函數(shù)需要兩個基本輸入?yún)?shù)：待處理的輸入圖像以及結(jié)構(gòu)元素對象。此外，imerode函數(shù)還可以接受3個可選參數(shù)：PADOPT(padopt)——影響輸出圖片的大小、PACKOPT(packopt).——說明輸入圖像是否為打包的二值圖像(二進制圖像)。M——指定原始圖像的行數(shù)。② 實驗內(nèi)容及步驟：

A、讀取圖像cameraman.tif（該圖像是Matlab當前目錄下自帶的圖片）R=imread('cameraman.tif');B、創(chuàng)建一個任意形狀的結(jié)構(gòu)元素對象 C=strel('arbitrary',eye(5));C、以圖像R和結(jié)構(gòu)元素C為參數(shù)調(diào)用imerode函數(shù)進行腐蝕操作。R1=imerode(R,C);D、顯示操作結(jié)果 imshow(R)； figure,imshow(C)； ③ 實驗結(jié)果：

圖像cameraman.tif 腐蝕前后的效果對比：

3、膨脹和腐蝕聯(lián)合操作 ① 實驗內(nèi)容及步驟： A、創(chuàng)建結(jié)構(gòu)元素： clear;close all C= strel('rectangle',[40 30]);

%結(jié)構(gòu)元素必須具有適當?shù)拇笮?既可以刪電流線又可以刪除矩形.B、使用結(jié)構(gòu)元素腐蝕圖像:

%將會刪除所有直線,但也會縮減矩形 R=imread('circbw.tif');R1=imerode(R,C);imshow(R1);figure,imshow(R);C、恢復矩形為原有大小,使用相同的結(jié)構(gòu)元素對腐蝕過的圖像進行膨脹.R2=imdilate(R1,C);figure,imshow(R2);② 實驗結(jié)果： 最終效果如下圖: 原始圖像 腐蝕后的圖像 膨脹后的圖像

4、邊緣檢測 ① 實驗原理：

對于一副灰度二進制圖像,如果圖像像素值為1,則該像素的狀態(tài)為ON,如果其像素值為0,則該像素的狀態(tài)為OFF。在一副圖像中，如果圖像某個像素滿足以下兩個條件：

1.該像素狀態(tài)為ON； 2.該像素鄰域中有一個或多個像素狀態(tài)為OFF。則認為該像素為邊緣像素。

Matlab中提供了專門的函數(shù)bwperim,可以用于判斷一副二進制圖像中的哪些像素為邊緣像素。② 實驗內(nèi)容及步驟：

以下程序代碼示例就是利用bwperim函數(shù)，對圖像circbw.tif進行邊緣檢測，其邊緣像素檢測效果如尾圖。clear;close all R=imread('circbw.tif');R1=bwperim(R);

imshow(R)； figure,imshow(R1)； ② 實驗結(jié)果：

心得體會：

本次實驗利用形態(tài)學運算—腐蝕和膨脹對圖像進行了處理。數(shù)學形態(tài)學的基本思想是：用具有一定形態(tài)的結(jié)構(gòu)元素去量度和提取圖像中的對應(yīng)形狀，以達到圖像分析和識別的目的。它的優(yōu)點有：有效濾除噪聲、保留圖像中原有信息、算法易于用并行處理方法有效實現(xiàn)等等。

總的來說，在此次實驗中遇到了較多問題，經(jīng)過對相關(guān)知識的復習以及對相關(guān)資料的查詢，才逐個地解決了。在今后的實驗課中要有充分的準備與預習，不懂的問題要提前查資料，這樣在實驗中才會胸有成竹，才不會手忙腳亂。

第三篇：matlab圖像處理小結(jié)

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 計算各個連通區(qū)域中的象素總數(shù) BoundingBox 包含相應(yīng)區(qū)域的最小矩形 Centroid 給出每個區(qū)域的質(zhì)心

MajorAxisLength 與區(qū)域具有相同標準二階中心矩（又叫標準差）的橢圓的長軸長度 MinorAxisLength 與區(qū)域具有相同標準二階中心矩的橢圓的短軸長度 Eccentricity 與區(qū)域具有相同標準二階中心矩的橢圓的離心率

Orientation 與區(qū)域具有相同標準二階中心矩的橢圓的長軸與x軸的交角 Image 二值圖像，與某區(qū)域具有相同大小的邏輯矩陣。

FilledImage 與上相同，唯一區(qū)別是這是個做了填充的邏輯矩陣！本例中和上面的沒有區(qū)別，只有 區(qū)域有空洞時才有明顯差別。

FilledArea 是標量，填充區(qū)域圖像中的 on 像素個數(shù)

ConvexHull 是p行2列的矩陣，包含某區(qū)域的最小凸多邊形 ConvexImage 二值圖像，用來畫出上述的區(qū)域最小凸多邊形 ConvexArea 是標量，填充區(qū)域凸多邊形圖像中的 on 像素個數(shù) EulerNumber 等于圖像中目標個數(shù)減去這些目標中空洞的個數(shù) Extrema 8行2列矩陣，八方向區(qū)域極值點

EquivDiameter 是標量，等價直徑：與區(qū)域具有相同面積的圓的直徑.計算公式為：sqrt(4*Area/pi)

Solidity 是標量，同時在區(qū)域和其最小凸多邊形中的像素比例。計算公式為： Area/ConvexArea，這也是個仿射特征，實際上反映出區(qū)域的固靠性程度。

Extent 是標量，同時在區(qū)域和其最小邊界矩形中的像素比例。計算公式為：Area除以邊界矩 形面積，這也是個仿射特征，實際上反映出區(qū)域的擴展范圍程度。

PixelIdxList p元向量，存儲區(qū)域像素的索引下標

PixelList p行ndims(L)列矩陣，存儲上述索引對應(yīng)的像素坐標 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 基于特定原則的區(qū)域選擇

當你要基于特定準則條件選擇某個區(qū)域時，將函數(shù) ismember 和 regionprops 聯(lián)合使用是很有用處的。例如：創(chuàng)建一個只包含面積大于80的二值圖像，用以下命令

idx = find([stats.Area] > 80);BW2 = ismember(L,idx);regionprops函數(shù)的擴展思路

在regionprops函數(shù)的基礎(chǔ)上，你可以使用它提供的基本數(shù)據(jù)來擴展它的功能，比如我就將區(qū)域的曲率數(shù)據(jù)和骨架數(shù)據(jù)作為它的另外屬性值來開發(fā)，從而希望它能用來做更細致的特征提取。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %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])%面積的變異系數(shù)

>> std([jdata.Eccentricity])/ mean([jdata.Eccentricity])%橢圓的變異系數(shù) >> 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])%面積的變異系數(shù)

>> std([graindata.Eccentricity])/ mean([graindata.Eccentricity])%橢圓的變異系數(shù)

>> 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;%設(shè)置堊白面積的下限 >> 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);%標記堊白米粒圖像，便于計數(shù) >> rate_chalky_grains=num_chalk/num*100;>> rate_chalky_grains

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.17 bwmorph函數(shù) >> 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));

%開運算（消除小物體）與閉運算（填充物體內(nèi)細小空洞）se=strel('disk',6);iopen=imopen(bw,se);iclose=imclose(bw,se);%腐蝕與膨脹聯(lián)合操作 %（1）創(chuàng)建結(jié)構(gòu)元素 se=strel('rectangle',[40 30]);%（2）使用結(jié)構(gòu)元素腐蝕圖像 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提取邊界，效果不大好?。?結(jié)論：無需拉氏銳化，也不必中值濾波，可直接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;%以下內(nèi)容畫在同一figure中 centr=[mexdata.Centroid];%尋找重心位置 nums=1:numobjects;for k = 1:numobjects soli=mexdata(k).Solidity;soli_string=sprintf('%2.2f',soli);%等價于轉(zhuǎn)字符串 % 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函數(shù)的用法 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)

%可以使用下面的方法創(chuàng)建相應(yīng)的向量： regionprops(L,'Area');allArea = [stats.Area];

%創(chuàng)建一個只包含面積大于80的二值圖像 idx = find([stats.Area] > 80);BW2 = ismember(L,idx);

%只顯示某個下標所對應(yīng)的物體圖像 bw2=ismember(L,N);figure,imshow(bw2);

%在調(diào)用regionprops之前必須將二值圖像轉(zhuǎn)變?yōu)闃俗⒕仃?L = bwlabel(BW);%或者

L = double(BW);

%將matlab數(shù)據(jù)寫到excel中 a=ones(3);success = xlswrite('c:/matlab/work/myworkbook.xls',a,'A2:C4')%將行矩陣轉(zhuǎn)換為列矩陣 a=[1 2 3 4 5 6];b=transpose(a);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2006.6.22球形物體的檢測和標識（循環(huán)檢測和標識算法）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圖像處理報告

一、設(shè)計題目

圖片疊加。

二、設(shè)計要求

將一幅禮花圖片和一幅夜景圖片做疊加運算，使達到煙花夜景的美圖效果。

三、設(shè)計方案

3.1、設(shè)計思路

利用matlab強大的圖像處理功能，通過編寫程序，實現(xiàn)對兩幅圖片的像素進行線性運算，利用灰度變換的算法使圖片達到預期的效果。

3.2、軟件介紹

MATLAB是matrix&laboratory兩個詞的組合，意為矩陣工廠（矩陣實驗室）。是由美國mathworks公司發(fā)布的主要面對科學計算、可視化以及交互式程序設(shè)計的高科技計算環(huán)境。它將數(shù)值分析、矩陣計算、科學數(shù)據(jù)可視化以及非線性動態(tài)系統(tǒng)的建模和仿真等諸多強大功能集成在一個易于使用的視窗環(huán)境中，為科學研究、工程設(shè)計以及必須進行有效數(shù)值計算的眾多科學領(lǐng)域提供了一種全面的解決方案，并在很大程度上擺脫了傳統(tǒng)非交互式程序設(shè)計語言（如C、Fortran）的編輯模式，代表了當今國際科學計算軟件的先進水平。

MATLAB和Mathematica、Maple并稱為三大數(shù)學軟件。它在數(shù)學類科技應(yīng)用軟件中在數(shù)值計算方面首屈一指。MATLAB可以進行矩陣運算、繪制函數(shù)和數(shù)據(jù)、實現(xiàn)算法、創(chuàng)建用戶界面、連接其他編程語言的程序等，主要應(yīng)用于工程計算、控制設(shè)計、信號處理與通訊、圖像處理、信號檢測、金融建模設(shè)計與分析等領(lǐng)域。MATLAB的基本數(shù)據(jù)單位是矩陣，它的指令表達式與數(shù)學、工程中常用的形式十分相似，故用MATLAB來解算問題要比用C，F(xiàn)ORTRAN等語言完成相同的事情簡捷得多，并且MATLAB也吸收了像Maple等軟件的優(yōu)點，使MATLAB成為一個強大的數(shù)學軟件。在新的版本中也加入了對C，F(xiàn)ORTRAN，C++，JAVA的支持。可以直接調(diào)用，用戶也可以將自己編寫的實用程序?qū)氲組ATLAB函數(shù)庫中方便自己以后調(diào)用，此外許多的MATLAB愛好者都編寫了一些經(jīng)典的程序，用戶直接進行下載就可以用。

3.3、常見簡單程序語句及算法分析

（1）CLC；清零。（2）CLEAR ； 清內(nèi)存。

（3）r=imread（‘路徑圖片名.jpg’）；讀入一幅圖片。（4）imshow（r）；顯示圖片r。（5）g=rgb2gray(r)；灰度轉(zhuǎn)變。

（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 此程序的主要功能是對兩幅圖片通過算法來實現(xiàn)疊加的效果，程序中的幾個變量都是像素點的值，通過三個循環(huán)使得兩幅圖片對應(yīng)的值線性相加，最大值應(yīng)該是以255輸出，超過255也是以255輸出。

程序2 J = imadjust(I,[low_in high_in], [low_out high_out])此程序變換的原理是：如果原圖像f(x, y)的灰度范圍是[m, M]，我們希望對圖像的灰度范圍進行線性調(diào)整，調(diào)整后的圖像g(x, y)的灰度范圍是[n, N]，那么下述變換：g(x,y)

N

n

f(x,y)

m

n就可以實現(xiàn)這一要求。MATLAB圖像處理工具箱中提供的imadjust函數(shù)，可以實現(xiàn)上述的線性變換對比度調(diào)整。

四、設(shè)計步驟（1）處理之前，我們先來看看兩幅原圖，一幅是帶有禮花的圖片，另一幅是一幅東方明珠的夜景圖。

圖 1

圖 2 通過圖片我們發(fā)現(xiàn)，禮花的圖片非常的亮，而夜景圖則顯得有些灰暗。我們推測，如果讓禮花和夜景的圖片疊加到一起，會不會由于禮花的亮度太大而掩蓋了城市的夜光，達不到我們想要的效果。

（2）打開matlab圖像處理軟件，7.0及以上版本就可以，打開之后，顯示界面如下（我的版本是2012b）。

圖3（3）新建script文件，點擊’New script‘，或點擊’New‘，選擇script。

圖4（4）在打開的界面之中輸入程序。圖片的位置是你所要用的圖片的準確位置，盡量寫詳細點，減少計算機的讀取時間，加快讀取速度。如圖所示。

圖5（5）點擊運行按鈕，即界面上的綠色按鈕。

圖6（6）查看效果，如下圖。

圖7 通過處理后的圖片，我們看到由于禮花太亮，完全將城市的夜景掩蓋住了，效果不是太理想。我們設(shè)想加入灰度變換函數(shù)，將禮花的圖片變的暗一點，將夜景的亮度提高一點，看看效果怎么樣，加入灰度變換程序，如下圖。

圖8（7）軟件調(diào)試后運行程序，如下圖。

圖9 結(jié)果顯示加入灰度變換的函數(shù)之后，圖片效果較之前好多了，達到了預期的目的。

五、軟件代碼

通過matlab進行圖像處理，代碼如下： clear;clc;A=imread('C:UsersAdministratorDesktop作業(yè)禮花.jpg');C=imresize(A,[300,400]);B=imread('C:UsersAdministratorDesktop作業(yè)夜景提亮天空中加入禮花.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作業(yè)禮花效果圖1.jpg')調(diào)試之后的程序代碼如下： clear;clc;A=imread('C:UsersAdministratorDesktop作業(yè)禮花.jpg');B=imread('C:UsersAdministratorDesktop作業(yè)夜景提亮天空中加入禮花.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作業(yè)禮花效果圖2.jpg')處理后的圖片如下：

圖 10

六、結(jié)果分析

通過兩幅圖的對比，發(fā)現(xiàn)第二幅圖片較第一幅，效果明顯增強。是由于加入灰度變換函數(shù)，使原圖的灰度值發(fā)生變化，以達到實際的效果。

效果對比圖11

七、心得體會

通過這次任務(wù)，以前在課堂上沒太聽明白怎么回事的東西，通過上網(wǎng)，查資料，以及用軟件處理，通通實踐了一遍，加深了對這門課程的認識和理解。Matlab是一款功能很強大的應(yīng)用軟件，它不僅可以對圖像進行處理，而且可以進行各種數(shù)字計算和符號計算功能，具有繪圖功能，語言體系等等。這次的任務(wù)我們樂在其中，喜悅的是看到了成果，內(nèi)心充滿了滿足感和成就感，就如同看到了圖片中的煙花，有種過年的感覺。不過這些都只是皮毛而已，要想真正地掌握它，還得更進一步地學習理論知識。

參考文獻：

【岡薩雷斯 數(shù)字圖像處理（MATLAB版）】[美] RafaelC.Gonzalez RichardE.Woods StevenL.Eddins 著 電子工業(yè)出版社。

【數(shù)字圖像處理及MATLAB實現(xiàn)】楊杰 主編 電子工業(yè)出版社。

第五篇：MATLAB實現(xiàn)數(shù)字信號處理

數(shù)字信號處理

說 明 書

目錄

一．摘要…………………………………3 二．課程設(shè)計目的………………………3 三．設(shè)計內(nèi)容……………………………3 四．設(shè)計原理……………………………4 4.1．語音信號的采集…………………………….4 4.2．濾波器……………………………………….4 4．21．IIR濾波器原理…………………………………….4 4.22．FIR濾波器原理………………………………………5 五．設(shè)計步驟……………………………6 5.1錄制女音………………………………………6 5.2采樣語音信號并畫出時域波形和頻譜圖……7 5.3采用雙線性變換法設(shè)計IIR濾波器…………10 5.4窗函數(shù)法設(shè)計FFR濾波器………………......12 5.5用IIR濾波器對信號進行濾波………………14 5.6用FIR濾波器對信號進行濾波………………16 5.7男女聲語音信號頻譜特點分析………………19 5.8有背景噪聲的信號分析………………………20 六．心得體會…………………………….22 七．參考文獻…………………………….23

一．摘要：

這次課程設(shè)計的主要目的是綜合運用本課程的理論知識進行頻譜分析以及濾波器設(shè)計，通過理論推導得出相應(yīng)結(jié)論，并利用MATLAB或者ＤＳＰ開發(fā)系統(tǒng)作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應(yīng)用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。通過對聲音的采樣，將聲音采樣后的頻譜與濾波。

MATLAB全稱是Matrix Laboratory，是一種功能強大、效率高、交互性好的數(shù)值和可視化計算機高級語言，它將數(shù)值分析、矩陣運算、信號處理和圖形顯示有機地融合為一體，形成了一個極其方便、用戶界面友好的操作環(huán)境。經(jīng)過多年的發(fā)展，已經(jīng)發(fā)展成為一種功能全面的軟件，幾乎可以解決科學計算中所有問題。MATLAB軟件還提供了非常廣泛和靈活的用于處理數(shù)據(jù)集的數(shù)組運算功能。

在本次課程設(shè)計中，主要通過MATLAB來編程對語音信號處理與濾波，設(shè)計濾波器來處理數(shù)字信號并對其進行分析。

二．課程設(shè)計目的：

綜合運用本課程的理論知識進行頻譜分析以及濾波器設(shè)計，通過理論推導得出相應(yīng)結(jié)論，并利用MATLAB作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應(yīng)用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。

三．設(shè)計內(nèi)容：

內(nèi)容：錄制一段個人自己的語音信號，并對錄制的信號進行采樣；畫出采樣后語音信號的時域波形和頻譜圖；給定濾波器的性能指標，采用窗函數(shù) 法和雙線性變換法設(shè)計濾波器，并畫出濾波器的頻率響應(yīng)；然后用自己設(shè)計的濾波器對采集的信號進行濾波，畫出濾波后信號的時域波形和頻譜，并對濾波前后的信號進行對比，分析信號的變化；回放語音信號；換一個與你性別相異的人錄制同樣一段語音內(nèi)容，分析兩段內(nèi)容相同的語音信號頻譜之間有什么特點；再錄制一段同樣長時間的背景噪聲疊加到你的語音信號中，分析疊加前后信號頻譜的變化，設(shè)計一個合適的濾波器，能夠把該噪聲濾除。

四．設(shè)計原理：

4.1．語音信號的采集

熟悉并掌握MATLAB中有關(guān)聲音（wave）錄制、播放、存儲和讀取的函數(shù)，在MATLAB環(huán)境中，有關(guān)聲音的函數(shù)有：

a：y=wavrecord(N,fs,Dtype);利用系統(tǒng)音頻輸入設(shè)備錄音，以fs為采樣頻率，默認值為11025，即以11025HZ進行采樣。Dtype為采樣數(shù)據(jù)的存儲格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度進行采樣，其它三種都是16位采樣結(jié)果轉(zhuǎn)換為指定的MATLAB數(shù)據(jù)；

b：wavplay(y,fs)；利用系統(tǒng)音頻輸出設(shè)備播放，以fs為播放頻率，播放語音信號y；

c：wavwrite((y,fs,wavfile);創(chuàng)建音頻文件； d：y=wavread(file);讀取音頻文件；

關(guān)于聲音的函數(shù)還有sound()；soundsc()；等。4.2濾波器： 4．21.IIR濾波器原理

沖激響應(yīng)不變法是使數(shù)字濾波器在時域上模擬濾波器，但是它們的缺點是產(chǎn)生頻率響應(yīng)的混疊失真，這是由于從s平面到z平面是多值的映射關(guān)系所造成的。

雙線性變換法是使數(shù)字濾波器的頻率響應(yīng)與模擬濾波器的頻率響應(yīng)相似的一種變換方法。為了克服多值映射這一缺點，我們首先把整個s平面壓縮變換到某一中介的s1平面的一條橫帶里，再通過變換關(guān)系將此橫帶變換到整個z平面上去，這樣就使得s平面與z平面是一一對應(yīng)的關(guān)系，消除了多值變換性，也 就消除了頻譜混疊現(xiàn)象。

雙線性法設(shè)計IIR數(shù)字濾波器的步驟：

1）將數(shù)字濾波器的頻率指標{ ?k}由Wk＝（2/T）*tan(wk),轉(zhuǎn)換為模擬濾波器的頻率指標{?k}.2)由模擬濾波器的指標設(shè)計H(s).3)由H(s)轉(zhuǎn)換為H(z)21?z?1H(z)?H(s)s?T1?z?1

4.22.FIR濾波器原理

FIR濾波器與IIR濾波器特點不同，IIR濾波器的相位是非線性的，若需線性相位則要采用全通網(wǎng)絡(luò)進行相位校正。而有限長單位沖激響應(yīng)（FIR）數(shù)字濾波器就可以做成具有嚴格的線性相位，同時又可以具有任意的幅度特性。

由于FIR系統(tǒng)的沖激響應(yīng)就是其系統(tǒng)函數(shù)各次項的系數(shù)，所以設(shè)計FIR濾波器的方法之一可以從時域出發(fā)，截取有限長的一段沖激響應(yīng)作為H（z）的系數(shù)，沖激響應(yīng)長度N就是系統(tǒng)函數(shù)H（z）的階數(shù)。只要N足夠長，截取的方法合理，總能滿足頻域的要求。這種時域設(shè)計、頻域檢驗的方法一般要反復幾個回合，不像IIR DF設(shè)計靠解析公式一次計算成功。給出的理想濾波器頻率響應(yīng)是，它是w的周期函數(shù)，周期

由傅立葉反變換導出，即

hd(n)?1Hd(ejw)ejwndw?2?，再將hd(n)與窗函數(shù)，因此可展開成傅氏級數(shù)w(n)相乘就可以得到h(n)。、的計算可采用傅氏變換的現(xiàn)成公式和程序，窗函數(shù)也是現(xiàn)成的。但整個設(shè)計過程不能一次完成，因為窗口類型和大小的選擇沒有解析公式可一次算，整個設(shè)計可用計算機編程來做。

窗函數(shù)的傅式變換W(ejω)的主瓣決定了H(ejω)過渡帶寬。W(ejω)的旁瓣大小和多少決定了H(ejω)在通帶和阻帶范圍內(nèi)波動幅度，常用的幾種窗函數(shù)有：

矩形窗

w(n)=RN(n)；

Hanning窗

；

Hamming窗

；

Blackmen窗

；

Kaiser窗。

式中Io(x)為零階貝塞爾函數(shù)。

五．設(shè)計步驟：

5.1錄制女音：

利用MATLAB中的函數(shù)錄制聲音。function nvyin()fs=11025;

%采樣頻率

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

開始錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double');

%錄制聲音3秒

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

錄音結(jié)束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs);

%播放錄音

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放錄音結(jié)束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音');

%聲音的存儲

5.2采樣語音信號并畫出時域波形和頻譜圖

讀取語音信號，畫出其時域波形和頻譜圖，與截取后的語音信號的時域波形和頻譜圖比較，觀察其變化。程序如下：

[x,fs,bits]=wavread('女音.wav');

%讀取聲音

N=length(x);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半 figure(1);subplot(2,2,1);

%把畫圖區(qū)域劃分為2行2列，指定第一個圖 plot(t, x);

%畫出聲音采樣后的時域波形 title('原女音信號的時域波形');

%給圖形加注標簽說明 xlabel('時間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,2,3);

%把畫圖區(qū)域劃分為2行2列，指定第三個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 title('FFT變換后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%添加網(wǎng)格

subplot(2,2,4);

%把畫圖區(qū)域劃分為2行2列，指定第四個圖 if y~=0

%判斷指數(shù)是否為0

plot(k,20*log10(abs(y(n))));

%畫信號頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid;

%添加網(wǎng)格

%實際發(fā)出聲音落后錄制動作半拍的現(xiàn)象的解決 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]);

%截取語音信號 N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

figure(2);subplot(2,2,1);

%把畫圖區(qū)域劃分為2行2列，指定第一個圖 plot(t,x1);

%畫出聲音采樣后的時域波形 title('截取后女音信號的時域波形');

%給圖形加注標簽說明 xlabel('時間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y1=fft(x1);

%對信號做N點FFT變換

subplot(2,2,3);

%把畫圖區(qū)域劃分為2行2列，指定第三個圖 k=(n-1)*f0;

%頻域采樣點

plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 title('FFT變換后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%添加網(wǎng)格

subplot(2,2,4);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 if y1~=0

%判斷指數(shù)是否為0

plot(k,20*log10(abs(y1(n))));

%畫信號頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid;

%添加網(wǎng)格

原女音信號的時域波形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

結(jié)果分析：

由原女音信號的時域波形可知錄取開始時實際發(fā)出聲音大概落后3500個采樣點，我們把前3500點去除即可解決實際發(fā)出聲音落后錄制動作半拍的現(xiàn)象。由原女音的的頻譜圖和截取后聲音的頻譜圖可看出，對聲音的截取并不會影響它們頻譜分布。

5.3采用雙線性變換法設(shè)計IIR濾波器：

人的聲音頻率一般在（1~~4）kHZ之間，則我們只需要設(shè)計一個帶通濾波器即可濾去聲音頻帶以外的無用噪聲，得到比較清晰的聲音。根據(jù)聲音的頻譜圖分析，設(shè)計一個帶通濾波器性能指標如下：

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;

%通帶波紋系數(shù)

Az=100;

%最小阻帶衰減

wp=[0.2 0.6];

%歸一化通帶數(shù)字截止頻率 wz=[0.1 0.7];

%歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估計契比雪夫I型濾波器階數(shù) [b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化

截 %止頻率，Ap通帶波動

[h,w]=freqz(b,a);

%求數(shù)字濾波器的復頻率響應(yīng) figure(1);subplot(2,1,1);plot(w/pi,abs(h));

%繪制數(shù)字濾波器的頻譜圖 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)');subplot(2,1,2);if abs(h)~=0

%判斷指數(shù)是否為0

plot(w/pi,20*log10(abs(h)));

%繪制數(shù)字濾波器頻譜的分貝圖 end grid;xlabel('omega/pi');ylabel('振幅(分貝)');title('契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)');契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)1振幅（幅值）0.5000.10.20.50.60.70.8?/?契比雪夫Ⅰ型帶通濾波器的幅頻響應(yīng)0.30.40.910振幅(分貝)-200-400-60000.10.20.30.40.5?/?0.60.70.80.91

5.4窗函數(shù)法設(shè)計FFR濾波器

線性相位FIR濾波器通常采用窗函數(shù)法設(shè)計。窗函數(shù)法設(shè)計FIR濾波器的基本思想是：根據(jù)給定的濾波器技術(shù)指標，選擇濾波器長度N和窗函數(shù)ω（n），使其具有最窄寬度的主瓣和最小的旁瓣。其核心是從給定的頻率特性，通過加窗確定有限長單位脈沖響應(yīng)序列h(n)。工程中常用的窗函數(shù)共有6種，即矩形窗、巴特利特（Bartlett）窗、漢寧（Hanning）窗、漢明（Hamming）窗、布萊克曼（Blackman）窗和凱澤（Kaiser）窗。

這次設(shè)計我采用的是布萊克曼來設(shè)計給定數(shù)字帶通濾波器的參數(shù)如下： wp1=0.3pi, wp2=0.6pi, wz1=0.2pi, wz2=0.7pi, Ap=1dB, Az=70dB 程序如下：

Ap=1;

%通帶波紋系數(shù) 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為奇數(shù) hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%理想帶通濾波器 wdWindow=blackman(N);

%布拉克曼窗 hr=wdWindow.*hdWindow';

%點乘

n=0:N-1;

%階數(shù) subplot(2,2,1);stem(n,wdWindow);

%繪制布拉克曼窗時域波形 xlabel('時間');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1);

%求濾波器頻率響應(yīng) subplot(2,2,3);plot(W/pi,abs(H))

%繪制濾波器頻域波形 xlabel('omega/pi');ylabel('振幅');title('FIR帶通濾波器幅頻特性');subplot(2,2,4);

if abs(H)~=0

%判斷指數(shù)是否為0

plot(W/pi,20*log10(abs(H)));

%畫濾波器頻譜的分貝圖 end xlabel('omega/');ylabel('振幅/分貝');title('FIR帶通濾波器幅頻特性');grid;

%添加網(wǎng)格 %---ideallp（）函數(shù)（非系統(tǒng)自有函數(shù)）在系統(tǒng)安裝目錄的WORK子目錄ideallp.m function hd = ideallp(wc,N);% 理想低通濾波器的脈沖響應(yīng)子程序 % hd = 點0 到 N-1之間的理想脈沖響應(yīng) % wc = 截止頻率(弧度)% N = 理想濾波器的長度

tao =(N-1)/2;

% 理想脈沖響應(yīng)的對稱中心位置 n = [0:(N-1)];

% 設(shè)定脈沖響應(yīng)長度 m = n-tao + eps;

% 加一個小數(shù)以避免零作除數(shù)

hd = sin(wc*m)./(pi*m);

% 理想脈沖響應(yīng)

布拉克曼窗1振幅0.500406080時間FIR帶通濾波器幅頻特性500振幅/分貝20FIR帶通濾波器幅頻特性1.51振幅-50-100-15000.5?/10.5000.5?/?1

5.5用IIR濾波器對信號進行濾波

用自己設(shè)計的IIR濾波器分別對采集的信號進行濾波，在Matlab中，IIR濾波器利用函數(shù)filter對信號進行濾波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半 y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 title('濾波前女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%iir帶通的代碼：

Ap=1;

%通帶波紋系數(shù)

Az=100;

%最小阻帶衰減

wp=[0.2 0.6];

%歸一化通帶數(shù)字截止頻率 wz=[0.1 0.7];

%歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估計契比雪夫I型濾波器階數(shù)

[b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化截止頻率，Ap通帶波動 x1=filter(b,a,x);

%對聲音濾波 wavplay(x1)wavwrite(x1,'IIR濾波后女音.wav');N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x1);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,2);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 title('l濾波后女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

濾波前女音的頻譜特性300幅值/A***030004000頻率/Hz濾波后女音的頻譜特性500060006040幅值/A***0頻率/Hz400050006000

結(jié)果分析：

由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3400Hz的頻譜成分?；胤耪Z音信號，由于低頻和高頻成分被濾除，聲音變得較低沉。

5.6用FIR濾波器對信號進行濾波

用自己設(shè)計的FIR濾波器分別對采集的信號進行濾波，在Matlab中,FIR濾波器利用函數(shù)fftfilt對信號進行濾波 程序如下：

[x,fs,bits]=wavread('女音.wav');N=length(x);

%計數(shù)讀取信號的點數(shù)

t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 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為奇數(shù) 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);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x1);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,2);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 title('l濾波后女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

濾波前女音的頻譜特性300200幅值/A***004000頻率/?l濾波后女音的頻譜特性500060006040幅值/A20005001000***03000頻率/Hz***0

結(jié)果分析：

由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3500Hz的頻譜成分。和用IIR濾波器濾波一樣，回放語音信號，由于低頻和高頻成分被濾除，聲音變得較低沉。5.7男女聲語音信號頻譜特點分析

換一個男音錄制同樣一段語音內(nèi)容，分析兩段內(nèi)容相同的語音信號頻譜之間有什么特點。程序如下：

[x,fs,bits]=wavread('女音.wav');

%讀取聲音

N=length(x);

%計數(shù)讀取信號的點數(shù) 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變換

%頻域采樣點

%把畫圖區(qū)域劃分為2行1列，指定第一個圖%繪制原始語音信號的幅頻響應(yīng)圖

%給圖形加注標簽說明

%添加網(wǎng)格

%讀取聲音

%計數(shù)讀取信號的點數(shù)

%信號的時域采樣點

%采樣間隔

%取信號的一半

%對信號做N點FFT變換

%頻域采樣點

%把畫圖區(qū)域劃分為2行1列，指定第二個圖%繪制原始語音信號的幅頻響應(yīng)圖

%給圖形加注標簽說明

%添加網(wǎng)格

axis([0 6000 0 300])；

%改變橫縱坐標便于比較頻譜圖

FFT變換后女音的頻譜特性300200幅值/A***00頻率/?FFT變換后男音的頻譜特性***200幅值/A***00頻率/?400050006000

結(jié)果分析：

通過比較上面女音頻譜圖和男音頻譜圖可知，男音的頻譜集中在低頻部分，高頻成分底，譜線較平滑，聲音聽起來低沉。5.8有背景噪聲的信號分析

從硬盤中把一段噪聲（頻譜能量集中在某個小范圍內(nèi)）疊加到語音信號中，分析疊加前后信號頻譜的變化，設(shè)計一個合適的濾波器，能夠把該噪聲濾除； 程序如下：

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);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 title('加噪聲后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格

%iir帶通濾波器的代碼：

Ap=1;

%通帶波紋系數(shù)

Az=70;

%最小阻帶衰減

wp=[0.2 0.7];

%歸一化通帶數(shù)字截止頻率 wz=[0.1 0.8];

%歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估計契比雪夫I型濾波器階數(shù)

[b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，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);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 k=(n-1)*f0;

%頻域采樣點

plot(k,abs(y1(n)));

%繪制原始語音信號的幅頻響應(yīng)圖 title('濾除噪聲后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格

加噪聲后聲音的頻譜特性3000幅值/A***0008000頻率/Hz濾除噪聲后聲音的頻譜特性***030幅值/A***000頻率/Hz80001000012000

結(jié)果分析

觀察加噪聲后聲音的頻譜圖可知，噪音頻率主要在4000Hz處，只要我們設(shè)計一個，濾波器濾除大概在4000Hz的頻譜即可，回放濾波后的語音信號，可證噪音基本濾除。

六．心得體會：

通過這次課程設(shè)計，讓我對MATLAB的基本應(yīng)用有了更深的了解，還有數(shù)字信號處理在MATLAB中的一些函數(shù)的用法。通過理論推導得出相應(yīng)結(jié)論，并利用MATLAB作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應(yīng)用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。

在這次實驗中，也遇到了很多問題，比如畫信號頻譜的分貝圖時（20*log10(abs(y))）指數(shù)為零時的處理。濾波器的設(shè)計也花了好大的功夫，剛開始不會設(shè)計參數(shù)，一頭霧水，通過同學的指導和討論，得知通過觀察信號的頻譜圖，看噪音頻率集中在那一部分，設(shè)計濾波器把其濾除即可。可反復設(shè)置參數(shù)直到濾波后語音信號的效果好為止。

七．參考文獻：

（1）《MATLAB LabVIEW SystemView》翁劍楓 葉志前 編著, 機械工業(yè)出版社；

（2）《MATLAB及在電子信息課程中的應(yīng)用》陳懷琛 吳大正 高西全編著，電子工業(yè)出版社；

（3）《MATLAB在數(shù)字信號處理中的應(yīng)用》（弟2版）薛年喜 編著，清華大學出版社；

（4）《MATLAB擴展編程》何強 何英

編著，清華大學出版社；（5）《MATLAB7簡明教程》吳清 曹輝林 編著，清華大學出版社；（6）MATLAB5.3精要.編程及高級應(yīng)用》程衛(wèi)國 馮峰 王雪梅 劉藝 編著，機械工程出版社。