1、應用Matlab對圖像信號進行頻譜分析及濾波實驗目的1. 鞏固所學的數字信號處理理論知識，理解信號的采集、處理、傳輸、顯示和存儲過程；2. 綜合運用專業及基礎知識，解決實際工程技術問題的能力；3. 學習資料的收集與整理，學會撰寫課程設計報告。實驗環境1. 微型電子計算機（PC）；2. 安裝Windows 10操作系統，MATLAB7.0，Formatfactory，繪畫板等開發工具。實驗原理 在Matlab軟件平臺下，讀取和顯示彩色圖像數據的相關函數和調用方法如下： 實驗內容和任務要求1. 選取一張彩色圖像（注意不能出現雷同，否則記為0分），提取圖像的灰度值，并顯示出灰度圖像。2. 在圖像中增

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

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

4、它們的頻譜繪制在同一張圖上作為對比。Matlab代碼clc;close all;clear all;%圖像的讀取以及轉換x=imread('hh.jpg');%讀取jpg圖像x1=rgb2gray(x);%生成M*N的灰度圖像矩陣M,N=size(x1);%求圖像規模%生成原始序列并求頻率響應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; %掃描時間間隔L=length(x3)

5、;%計算序列長度n=0:L-1;Xk3=fft(x3);%快速傅里葉變換Xf3=fftshift(Xk3);%中心對稱變換f=(n/L-1/2)*fs;%橫坐標變幻%生成帶有噪聲的序列并求頻率響應fz1=356;fz2=383;fz3=427;%三個噪聲頻率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);%設計IIR濾波器并分析相關指標wp=250*2/fs;ws=300*2/fs;Rp=3;Rs=20;%通帶截止頻率250Hz,阻帶截止頻率300Hz,通帶衰減3dB,阻帶衰減20dBNm,Wc=buttord(wp,ws,Rp,Rs);%計算滿足指標的最小階數以及3dB截止頻率b,a=butter(Nm,Wc); %計算濾波器的分子分母系數H=freqz(b,a,f*2*pi/fs);%計算濾波器頻率響應mag=abs(H);pha=angle(H);%幅度響應和相位響應mag1=20*log(mag+eps)/max(mag); %將幅頻響應轉

7、化為dB的形式%用IIR濾波器對帶噪序列進行濾波并求頻率響應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%設計FIR濾波器并分析相關指標wc=280*2/fs; %6dB截止頻率280kHzfx=0 wc wc 1;m=1 1 0 0;%理想頻幅響應b1=fir2(40,fx,m,hamming(41);%計算FIR濾波器多項式系數（不聲明窗默認為Hamming窗）H1=freqz(b1,1,f*2*pi

8、/fs);%計算濾波器頻率響應mag2=abs(H1);pha1=angle(H1);%幅度響應和相位響應mag3=20*log(mag2+eps)/max(mag2); %將幅頻響應轉化為dB的形式%用FIR濾波器對帶噪序列進行濾波并求頻率響應x8=filter(b1,1,x4);%用FIR進行濾波grd=grpdelay(b1,1,f*2*pi/fs);%計算群延時K=round(grd(1);%修正因子(冒號操作做索引時，需要整型數操作)x8=x8(K+1):L),x8(1:K);%對群延遲進行修正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濾波器幅度響應');xlabel('f/Hz');ylabel('幅度');subplot(3,2,2);plot(f,mag2);grid;title('FIR濾波器幅度響應');xlabel('f/Hz');ylabel('幅度');subplot(3,2,3);plot(f,pha);grid;title('IIR濾波器相位響應');xlabel(

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

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