應(yīng)用Matlab對(duì)圖像信號(hào)進(jìn)行頻譜分析及濾波

1、應(yīng)用Matlab對(duì)圖像信號(hào)進(jìn)行頻譜分析及濾波實(shí)驗(yàn)?zāi)康?. 鞏固所學(xué)的數(shù)字信號(hào)處理理論知識(shí)，理解信號(hào)的采集、處理、傳輸、顯示和存儲(chǔ)過(guò)程；2. 綜合運(yùn)用專(zhuān)業(yè)及基礎(chǔ)知識(shí)，解決實(shí)際工程技術(shù)問(wèn)題的能力；3. 學(xué)習(xí)資料的收集與整理，學(xué)會(huì)撰寫(xiě)課程設(shè)計(jì)報(bào)告。實(shí)驗(yàn)環(huán)境1. 微型電子計(jì)算機(jī)（PC）；2. 安裝Windows 10操作系統(tǒng)，MATLAB7.0，F(xiàn)ormatfactory，繪畫(huà)板等開(kāi)發(fā)工具。實(shí)驗(yàn)原理 在Matlab軟件平臺(tái)下，讀取和顯示彩色圖像數(shù)據(jù)的相關(guān)函數(shù)和調(diào)用方法如下： 實(shí)驗(yàn)內(nèi)容和任務(wù)要求1. 選取一張彩色圖像（注意不能出現(xiàn)雷同，否則記為0分），提取圖像的灰度值，并顯示出灰度圖像。2. 在圖像中增

2、加正弦噪聲信號(hào)(自己設(shè)置幾個(gè)頻率的正弦信號(hào))，顯示出加入噪聲信號(hào)后的灰度圖像。3. 給定濾波器的性能指標(biāo)，分別設(shè)計(jì)FIR和IIR數(shù)字濾波器，并畫(huà)出濾波器的幅頻響應(yīng)曲線(xiàn)。4. 用自己設(shè)計(jì)的濾波器對(duì)含噪聲圖像信號(hào)進(jìn)行濾波，顯示出濾波后的灰度圖像。5. 對(duì)原始灰度圖像、加入噪聲信號(hào)的灰度圖像和濾波后的灰度圖像進(jìn)行頻譜分析和對(duì)比，分析信號(hào)的變化。實(shí)驗(yàn)分析本實(shí)驗(yàn)要求用Matlab軟件完成對(duì)圖像信息的處理。對(duì)于任務(wù)1，這里采用了一張jpg格式的張學(xué)友新專(zhuān)輯醒著做夢(mèng)的封面圖片，保存在Matlab的work文件夾下。采用imread()函數(shù)讀取，并利用rgb2gray()函數(shù)將其轉(zhuǎn)化為二維的灰度圖像（原始的數(shù)

3、據(jù)類(lèi)型是unit8型，需要將其轉(zhuǎn)化為可用于計(jì)算的double類(lèi)型），并利用imshow()函數(shù)將其顯示出來(lái)；對(duì)于任務(wù)2，在加入噪聲前，需要先將二維數(shù)據(jù)利用循環(huán)嵌套語(yǔ)句轉(zhuǎn)化為一維數(shù)據(jù)，然后加入三個(gè)高頻噪聲，再利用循環(huán)嵌套語(yǔ)句轉(zhuǎn)化為二維的數(shù)據(jù)，利用imshow()函數(shù)顯示出來(lái)；對(duì)于任務(wù)3，這里分別設(shè)計(jì)了滿(mǎn)足一定指標(biāo)的IIR低通濾波器（巴特沃斯）和FIR低通濾波器（哈明窗）并對(duì)其相關(guān)指標(biāo)進(jìn)行了分析。對(duì)于任務(wù)4，利用任務(wù)3中設(shè)計(jì)好的兩個(gè)濾波器分別對(duì)加噪后的灰度圖像進(jìn)行濾波（filter()函數(shù)），并分別顯示濾波后的灰度圖像；對(duì)于任務(wù)5，利用快速傅里葉變換算法（FFT）對(duì)各階段數(shù)據(jù)分別進(jìn)行頻譜分析，并將

4、它們的頻譜繪制在同一張圖上作為對(duì)比。Matlab代碼clc;close all;clear all;%圖像的讀取以及轉(zhuǎn)換x=imread('hh.jpg');%讀取jpg圖像x1=rgb2gray(x);%生成M*N的灰度圖像矩陣M,N=size(x1);%求圖像規(guī)模%生成原始序列并求頻率響應(yīng)x2=im2double(x1);x3=zeros(1,M*N);%初始化for i=1:M for j=1:N x3(N*(i-1)+j)=x2(i,j); endend %將M*N維矩陣變成1維矩陣fs=1000;%掃描頻率1kHzT=1/fs; %掃描時(shí)間間隔L=length(x3)

5、;%計(jì)算序列長(zhǎng)度n=0:L-1;Xk3=fft(x3);%快速傅里葉變換Xf3=fftshift(Xk3);%中心對(duì)稱(chēng)變換f=(n/L-1/2)*fs;%橫坐標(biāo)變幻%生成帶有噪聲的序列并求頻率響應(yīng)fz1=356;fz2=383;fz3=427;%三個(gè)噪聲頻率xz=0.4*sin(2*pi*fz1*n*T)+0.7*sin(2*pi*fz2*n*T)+0.5*sin(2*pi*fz3*n*T);%噪聲序列x4=x3+xz;%加入噪聲的序列x5=zeros(M,N);for i=1:M for j=1:N x5(i,j)=x4(N*(i-1)+j); endend %一維變M*N矩陣Xk4=fft

6、(x4);Xf4=fftshift(Xk4);%設(shè)計(jì)IIR濾波器并分析相關(guān)指標(biāo)wp=250*2/fs;ws=300*2/fs;Rp=3;Rs=20;%通帶截止頻率250Hz,阻帶截止頻率300Hz,通帶衰減3dB,阻帶衰減20dBNm,Wc=buttord(wp,ws,Rp,Rs);%計(jì)算滿(mǎn)足指標(biāo)的最小階數(shù)以及3dB截止頻率b,a=butter(Nm,Wc); %計(jì)算濾波器的分子分母系數(shù)H=freqz(b,a,f*2*pi/fs);%計(jì)算濾波器頻率響應(yīng)mag=abs(H);pha=angle(H);%幅度響應(yīng)和相位響應(yīng)mag1=20*log(mag+eps)/max(mag); %將幅頻響應(yīng)轉(zhuǎn)

7、化為dB的形式%用IIR濾波器對(duì)帶噪序列進(jìn)行濾波并求頻率響應(yīng)x6=filter(b,a,x4);%用IIR濾波Xk6=fft(x6);Xf6=fftshift(Xk6);x7=zeros(M,N);for i=1:M for j=1:N x7(i,j)=x6(N*(i-1)+j); endend%設(shè)計(jì)FIR濾波器并分析相關(guān)指標(biāo)wc=280*2/fs; %6dB截止頻率280kHzfx=0 wc wc 1;m=1 1 0 0;%理想頻幅響應(yīng)b1=fir2(40,fx,m,hamming(41);%計(jì)算FIR濾波器多項(xiàng)式系數(shù)（不聲明窗默認(rèn)為Hamming窗）H1=freqz(b1,1,f*2*pi

8、/fs);%計(jì)算濾波器頻率響應(yīng)mag2=abs(H1);pha1=angle(H1);%幅度響應(yīng)和相位響應(yīng)mag3=20*log(mag2+eps)/max(mag2); %將幅頻響應(yīng)轉(zhuǎn)化為dB的形式%用FIR濾波器對(duì)帶噪序列進(jìn)行濾波并求頻率響應(yīng)x8=filter(b1,1,x4);%用FIR進(jìn)行濾波grd=grpdelay(b1,1,f*2*pi/fs);%計(jì)算群延時(shí)K=round(grd(1);%修正因子(冒號(hào)操作做索引時(shí)，需要整型數(shù)操作)x8=x8(K+1):L),x8(1:K);%對(duì)群延遲進(jìn)行修正Xk8=fft(x8);Xf8=fftshift(Xk8);x9=zeros(M,N);f

9、or i=1:M for j=1:N x9(i,j)=x8(N*(i-1)+j); endend%繪制圖像figure(1);subplot(2,2,1);imshow(x2);title('原始灰度圖像');subplot(2,2,2);imshow(x5);title('加入噪聲后灰度圖像');subplot(2,2,3);imshow(x7);title('IIR濾波器濾波后灰度圖像');subplot(2,2,4);imshow(x9);title('FIR濾波器濾波后灰度圖像');figure(2);subplot(4,

10、1,1);plot(f,abs(Xf3)*2/L,'r-');title('原始灰度圖像幅度譜');subplot(4,1,2),plot(f,abs(Xf4)*2/L,'r-');title('加入噪聲后灰度圖像幅度譜');subplot(4,1,3);plot(f,abs(Xf6)*2/L,'r-');title('IIR濾波器濾波后灰度圖像幅度譜');subplot(4,1,4);plot(f,abs(Xf8)*2/L,'r-');title('FIR濾波器濾波后灰度圖

11、像幅度譜');figure(3);subplot(3,2,1);plot(f,mag);grid;title('IIR濾波器幅度響應(yīng)');xlabel('f/Hz');ylabel('幅度');subplot(3,2,2);plot(f,mag2);grid;title('FIR濾波器幅度響應(yīng)');xlabel('f/Hz');ylabel('幅度');subplot(3,2,3);plot(f,pha);grid;title('IIR濾波器相位響應(yīng)');xlabel(

12、9;f/Hz');ylabel('相位');subplot(3,2,4);plot(f,pha1);grid;title('FIR濾波器相位響應(yīng)');xlabel('f/Hz');ylabel('相位');subplot(3,2,5);plot(f,mag1);grid;title('IIR濾波器幅度響應(yīng)(dB)');xlabel('f/Hz');ylabel('幅度/dB');subplot(3,2,6);plot(f,mag3);grid;title('FIR濾波

13、器幅度響應(yīng)(dB)');xlabel('f/Hz');ylabel('幅度/dB');實(shí)驗(yàn)結(jié)果及分析從實(shí)驗(yàn)結(jié)果來(lái)看，濾波效果還是相當(dāng)不錯(cuò)的。在滿(mǎn)足相同的指標(biāo)情況下FIR濾波器所需要的階數(shù)遠(yuǎn)高于IIR濾波器，但是FIR濾波器的相位響應(yīng)是線(xiàn)性的，濾波器通帶群延時(shí)是常數(shù)，而IIR濾波器是非線(xiàn)性的，濾波器通帶群延時(shí)非常數(shù)。FIR濾波器產(chǎn)生的群延遲樣本數(shù)可以計(jì)算出來(lái)，并進(jìn)行調(diào)整。IIR濾波器則不可以。調(diào)整代碼為：grd=grpdelay(b1,1,f*2*pi/fs);%計(jì)算群延時(shí)K=round(grd(1);%修正因子(冒號(hào)操作做索引時(shí)，需要整型數(shù)操作)x8=x8(K+1):L),x8(1:K);%對(duì)群延遲進(jìn)行修正FIR濾波器濾波后圖像群延遲調(diào)整前后對(duì)比：將兩種濾波器濾波后的圖片放大后對(duì)比：對(duì)比可發(fā)現(xiàn)IIR濾波器濾波后使圖片產(chǎn)生了重影，而FIR濾波器沒(méi)有明顯失真。對(duì)于語(yǔ)音系統(tǒng)，対相位要求不是主要的，因此，選用IIR濾波器較為合適，可以充分發(fā)揮其經(jīng)濟(jì)和高效的特點(diǎn)；圖像信號(hào)和數(shù)據(jù)傳輸?shù)纫圆?/p>