現代圖像、語聲、資料通信對線性相位的要求是普遍的。正是此原因，使得具有線性相位的FIR數字濾波器得到大力發展和廣泛應用。
在實際進行數字信号處理時，往往需要把信号的觀察時間限制在一定的時間間隔内，隻需要選擇一段時間信号對其進行分析。取用有限個資料，即将信号資料截斷的過程，就等于将信号進行加窗函數操作。這樣操作以後，常常會發生頻譜分量從其正常頻譜擴充開來的現象，即所謂的“頻譜洩漏”。當進行離散傅立葉變換時，時域中的截斷是必需的，是以洩漏效應也是離散傅立葉變換所固有的，必須進行抑制。而要對頻譜洩漏進行抑制，可以通過窗函數權重抑制DFT的等效濾波器的振幅特性的副瓣，或用窗函數權重使有限長度的輸入信号周期延拓後在邊界上盡量減少不連續程度的方法實作。
數字帶通濾波器是一種用來過濾時間離散信号的數字系統，通過對抽樣資料進行數學處理來達到頻域濾波的目的。根據其機關重新開機響應函數的時域特性可分為兩類：無限沖擊響應濾波器(IIR)，有限沖擊響應濾波器(FIR)。與IIR濾波器相比，FIR的實作是遞歸的，總是穩定的；更重要的是，FIR濾波器在滿足幅頻響應要求的同時，可以獲得嚴格的線性相位特性。是以，它在高保真的信号處理，如信号音頻，圖像處理，資料傳輸等領域得到廣泛的應用。
數字fir濾波器的設計方法有很多種。如窗函數法設計，頻率采樣設計法和最優化設計法等。
2.程式設計及運作結果
2.1 常見基本運算
2.1.1 極限的計算
MATLAB提供的指令limit()可以完成極限運算，其調用格式如下: limit(F,x,a,'left')
該指令對表達式F求極限，獨立變量x從左邊趨近于a,函數中除F外的參數均可省略，'left'可換成'right'。
例：求極限S= [1+1/x]
代碼如下：
clear;
F=sym('1+1/x')
limit(F,'x',inf,'left')
運作結果：
2.1.2微分的計算
MATLAB提供的函數diff()可以完成對給定函數求導函數的運算，其調用格式如下：diff(fun,x,n)，其意義是求函數fun關于變量x的n階導數，n為1時可省略。這裡的fun用上例的後一種方式來定義較為妥當。
例：求函數y=log的一階導。
代碼如下：
clear;
syms x
y=log(1/x);
dy=diff(y,x)
運作結果：
2.1.3積分的計算
MATLAB中主要用int進行符号積分，用trapz,dblquad,quad,quad8等進行數值積分。
R=int(s,v)%對符号表達式中指定的符号變量v計算不定積分，表達式R隻是表達式函數s的一個原函數，後面沒有帶任何常數C。
R=int(s)%對符号表達式s的确定的符号變量計算不定積分。
R=int(s,a,b)%符号表達式s的定積分，a,b分别為積分的上，下限。
trapz(x,y)梯形積分法，x時表示積分區間的離散化向量，y是與x同維數的向量，表示被積函數，z傳回積分值。
fblquad(‘fun',a,b,c,d)矩形區域二重積分，fun表示被積函數的M函數名，a,b分别為x的上下限，c,d分别為y的上下線。
例1：（不定積分）用符号積分指令int計算
代碼如下：
clear;
syms x
int(x^2*sin(x))
運作結果：
例2：（定積分）計算數值積分
①用符号積分指令int，代碼如下：
clear;
syms x;
int(x^4,x,-2,2)
運作結果：
②改用梯形積分法指令trapz計算積分，代碼如下：
clear;
x=-2:0.01:1;
y=x.^2;
trapz(x,y)
運作結果：
2.1.4級數的計算
MATLAB中主要用symsun,taylor求級數的和及Taylor展開。其中Symsum(s,v,a,b)表達式s關于變量v從a到b求和Taylor(f,a,n)将函數f在a點展開為n-1階Taylor多項式
例：用Symsum計算
代碼如下：
clear;clc;
syms x y
z=1/(x^3);
symsum(z,x,1,9)
運作結果：
2.1.5求解代數方程
MATLAB指令輸人格式：solve('equ1','equ2',...'equN')，其中eqni表示第i個方程。
例1：求解方程 x^2+b*x+c=0
代碼如下：
solve('a*x^2+b*x+c')
運作結果：
2.1.6求解常微分方程
MATLAB中主要用dsolve求符号解析解，ode45,ode23,ode15s求數值解。s=dsolve(‘方程1', ‘方程2',…,'初始條件1'，'初始條件2' …,'自變量') 用字元串方程表示，自變量預設值為t。導數用D表示，2階導數用D2表示，以此類推。S傳回解析解。在方程組情形，s為一個符号結構。
[tout,yout]=ode45(‘yprime',[t0,tf],y0) 采用變步長四階Runge-Kutta法和五階Runge-Kutta-Felhberg法求數值解，yprime是用以表示f(t,y)的M檔案名，t0表示自變量的初始值，tf表示自變量的終值，y0表示初始向量值。輸出向量tout表示節點(t0,t1, …,tn)T,輸出矩陣yout表示數值解，每一列對應y的一個分量。若無輸出參數，則自動作出圖形。
例：求解微分方程的特解
在初始條件下的特解：y(0)=1
代碼如下：
dsolve('(x^2-1)*Dy+2*x*y-cos(x)=0','y(0)=1','x')
運作結果如下：
2.2 矩陣基本計算
2.2.1矩陣的最大值
求一個向量X的最大值的函數有兩種調用格式，分别是：
(1) y=max(X)：傳回向量X的最大值存入y，如果X中包含複數元素，則按模取最大值。
(2) [y,I]=max(X)：傳回向量X的最大值存入y，最大值的序号存入I，如果X中包含複數元素，則按模取最大值。
例：求出矩陣[12,52,-25,63,45,1,0,12]中的最大值及其所在位置
代碼如下：
x=[12,52,-25,63,45,1,0,12];
[y,l]=max(x) %求向量x中的最大值及其該元素的位置
運作結果：
2.2.2矩陣的最小值
求一個向量X的最小值的函數有兩種調用格式，分别是：
(1) y=min(X)：傳回向量X的最小值存入y，如果X中包含複數元素，則按模取最大值。
(2) [y,I]=min (X)：傳回向量X的最小值存入y，最小值的序号存入I，如果X中包含複數元素，則按模取最小值。
例：求出矩陣[-5,85,74,42,-2,3,0,-568]中的最小值及其所在位置
代碼如下：
x=[-5,85,74,42,-2,3,0,-568];
[y,l]=min(x) %求向量x中的最小值及其該元素的位置
運作結果：
2.2.3矩陣的均值
mean(A,1)表示對列取平均，mean(A,2）表示對行取平均，mean(A）則預設為mean(A,1）
例：求[45,47,41,42,46,48,42,45,44]均值
代碼如下：
A=[45,47,41;42,46,48;42,45,44]
mean(A,1)
mean(A,2)
運作結果如下：
2.2.4矩陣的方差
var(a); % 預設來求
var(a, 0); % 預設的公式（用N-1）
var(a, 1); % 另外的公式（用N）
var(a, 0, 1); % 對每列操作
var(a, 0, 2); % 對每行操作
var(a'); % 檢驗
var(a(:)); % 通過直接通路矩陣的存儲，來對矩陣進行操作
例：求方差
代碼如下：
a=rand(3,4) %(0,1)之間産生随機數
b=var(a)
c=var(a(:))
運作結果：
2.2.5矩陣的轉置
矩陣轉置用符号“`”來表示和實作, 如故Z是複數矩陣，則Z`為它們的複數共轭轉置矩陣，非共轭轉置矩陣使用Z.`或conj(Z`)。
例：求[33,34,33,36;32,33,31,33;38,37,30,37;34,36,37,34]的轉置。
代碼如下：
A=[33,34,33,36;32,33,31,33;38,37,30,37;34,36,37,34]
B=A'
運作結果：
2.2.6矩陣的逆
MATLAB中求矩陣的逆用函數inv(A)
例：求矩陣A=[13,14,15,16;12,13,11,13;18,17,10,17;14,16,17,14]的逆
代碼如下：
A=[13,14,15,16;12,13,11,13;18,17,10,17;14,16,17,14]
B=inv(A)
運作結果：
2.2.7矩陣的行列式
調用函數det格式：d = det(X)
例：求矩陣[13,4,1,2;1,2,1,3;8,7,2,7;14,6,7,1]的行列式
代碼如下：
A=[13,4,1,2;1,2,1,3;8,7,2,7;14,6,7,1]
B=det(A)
運作結果：
2.2.8矩陣的特征值計算
MATLAB中用函數eig()來進行求特征值和特征向量，其調用格式如下：d=eig(A),傳回矩陣A的所有特征值。
[V,D] = eig(A)，傳回矩陣A 的特征值和特征向量，它們滿足如下關系：A*V = V*D。
[V,D] = eig( A ,'nobalance') ，在求解特征值和特征向量時不采用初期的平衡步驟。一般來說平衡步驟對輸入矩陣進行調整，這使得計算出的特征值和特征向量更加準确。然而如果輸入矩陣中确實含有值很小的元素（可能會導緻截斷誤差），平衡步驟有可能加大這種誤差，進而得到錯誤的特征值和特征向量。
d = eig(A,B)，傳回矩陣A 和B 的廣義特征值。
[V,D] = eig(A,B)，傳回矩陣A 和B 的廣義特征值和廣義特征向量。
[ V,D] = eig( A,B,flag)，flag 有'chol'和‘qz' 兩種值。當flag ='chol' 時，計算廣義特征值采用B 的Cholesky 分解來實作。當flag ='qz'時，無論矩陣的對稱性如何，都采用QZ算法來求解廣義特征值。
例：求矩陣[33,24,15,26;12,23,31,33;38,27,20,37;14,16,22,14]的特征值和特征向量
代碼如下：
A=[33,24,15,26;12,23,31,33;38,27,20,37;14,16,22,14]
[V D]=eig(A)
eig(A)
運作結果：
2.2.9矩陣的相乘
例：求[1,2,3,4;5,6,7,8;9,1,2,3]和[2;3;6;7]的乘積。
代碼如下：
A=[1,2,3,4;5,6,7,8;9,1,2,3]
B=[2;3;6;7]
C=A*B
運作結果：
2.2.10矩陣的右除和左除
MATLAB中左除和右除分别用“\”和“/”表示
例：求[11,12,13;14,15,16;17,18,19]和[15,14,16;17,16,18;12,14,19]的左除和右除
代碼如下：
A=[11,12,13;14,15,16;17,18,19];
B=[15,14,16;17,16,18;12,14,19];
C=A\B
D=A/B
運作結果：
2.2.11矩陣的幂運算
例：求矩陣[1,2,3,4;5,6,7,8;7,5,4,3]的3次幂
代碼如下：
A=[1,2,3,4;5,6,7,8;7,5,4,3;4,5,6,7]
B=3;
C=A^B
運作結果：
2.3 多項式基本計算
2.3.1多項式加減運算
例：計算[31,12,13,16,87]和[5,9,6,6,3]的加減乘除
程式設計如下：
A=[31,12,13,16,87];
B=[5,9,6,6,3];
C=A+B
D=A-B
運作結果：
2.3.2多項式乘除運算
conv(A,B)用于多項式的乘法；
Decnov(A,B)用于多項式相除求餘式。
例：
P1=[3,4,4,2,3,4,5]
P2=[6,7,8]
P=conv(P1,P2)
P3=deconv(P,P1)
運作結果：
2.3.3多項式求導
polyder(P) %對多項式P進行微分運算
例：
代碼如下
P=[1,2,3,4,5]
P1=poly2sym(P)
Q=polyder(P)
運作結果如下：
2.3.4求根和求值運算
polyval(P) %按數組計算規則計算多項式P的值
polyvalm(P) %按矩陣計算規則計算多項式P的值
例：求多項式
在x=2.5時的值。
代碼如下：
P=[1,4,4]
Pv=polyval(P,0.5)
運作結果：
2.3.5多項式的部分分式展開
調用函數residue格式：[r,p,k] = residue(b,a)，[b,a] = residue(r,p,k)
其中，b與a分别為兩個多項式。其中a為分式之分母，b該分式之分子，兩者均為向量型式；而[r,p,k]等則為行向量，分别代表展開後之分子、分母及餘數之系數，若能除盡則k項應為零。r與p分别為餘數與極數，其個數應相等，但應比a之個數少一。
代碼如下：
b=[3,4,5];
a=[8,-6,3];
[r,p,k]=residue(b,a)
運作結果：
2.3.6多項式的拟合
MATLAB軟體提供了基本的曲線拟合函數的指令．多項式函數拟合：a=polyfit(xdata,ydata,n)其中n表示多項式的最高階數，xdata，ydata為将要拟合的資料，它是用數組的方式輸入．輸出參數a為拟合多項式 y=a1xn+...+anx+a n+1的系數
例：已知五個資料點：[1,5.5],[2,4.3],[3,12.8],[4,2.9],[5,4.98]，試畫出這五個點拟合的三次曲線。
代碼如下：
x=[1,2,3,4,5];y=[5.5,4.3,1.28,0.29,0.498];
p=polyfit(x,y,3)
x2=1:0.1:5;y2=polyval(p,x2);
plot(x,y,'o',x2,y2) ,grid on
運作結果：
2.3.7插值運算
yi = interp1(x,Y,xi,method) %用指定的算法計算插值：
'nearest'：最近鄰點插值，直接完成計算；
'linear'：線性插值（預設方式），直接完成計算；
'spline'：三次樣條函數插值。
'cubic'： 分段三次Hermite插值。
對于超出x範圍的xi的分量，使用方法'nearest'、'linear'、'v5cubic'的插值算法，相應地将傳回NaN。對其他的方法，interp1将對超出的分量執行外插值算法。
yi = interp1(x,Y,xi,method,'extrap')
yi = interp1(x,Y,xi,method,extrapval) %确定超出x範圍的xi中的分量的外插值extrapval，其值通常取NaN或0。
ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method計算二維插值：
'linear'：雙線性插值算法（預設算法）；
'nearest'：最臨近插值；
'spline'：三次樣條插值；
'cubic'：雙三次插值。
例：用二維三次插值函數模拟函數z(x,y)=sin(x)+cos(y)。
代碼如下：
clear;
clc;
x=linspace(0,2*pi,6);
y=linspace(0,2*pi,6);
[X,Y]=meshgrid(x,y);
Z=cos(X)+sin(Y);
surf(X,Y,Z)
x1=linspace(0,2*pi,36);
y1=x1;
[X1,Y1]=meshgrid(x1,y1);
Z1=interp2(X,Y,Z,X1,Y1,'cubic');
figure(2)
surf(X1,Y1,Z1)
運作結果：
3.基于MATLAB的圖像濾波設計
3.1讀入圖像并分别加入高斯噪聲、椒鹽噪聲和乘性噪聲，并比較結果
函數使用imnoise函數進行圖檔的噪聲加入，其調用格式如下：
J=imnoise(I,type,parameters)其中，type是噪聲的類型，有高斯噪聲，椒鹽噪聲，乘性噪聲。類型名分别為：guassian,salt & pepper,speckle，parameters可表示噪聲密度D。
代碼如下：
clear;
I=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
subplot(2,2,1);
imshow(I);
title('原始圖像');
J2=imnoise(I,'gaussian',0.1);
subplot(2,2,3)
imshow(J2)
title('高斯噪聲圖像');
J1=imnoise(I,'salt & pepper',0.1);
subplot(2,2,2)
imshow(J1)
title('椒鹽噪聲圖像');
J3=imnoise(I,'speckle',0.1);
subplot(2,2,4)
imshow(J3)
title('乘性噪聲圖像');
運作結果，對比如下：
噪聲密度D不同對圖像幹擾程度對比，程式如下：
高斯噪聲：
clear;
I=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
J1=imnoise(I,'gaussian',0.03);
subplot(1,3,1);
imshow(J1);
title('d=0.03時的圖像');
J2=imnoise(I,'gaussian',0.15);
subplot(1,3,2)
imshow(J2)
title('d=0.15時的圖像');
J3=imnoise(I,'gaussian',0.5);
subplot(1,3,3)
imshow(J3)
title('d=0.5時的圖像');
運作結果對比如下：
椒鹽噪聲：
clear;
I=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
J1=imnoise(I,'salt & pepper',0.03);
subplot(1,3,1);
imshow(J1);
title('d=0.03時的圖像');
J2=imnoise(I,'salt & pepper',0.15);
subplot(1,3,2)
imshow(J2)
title('d=0.15時的圖像');
J3=imnoise(I,'salt & pepper',0.5);
subplot(1,3,3)
imshow(J3)
title('d=0.5時的圖像');
運作結果對比如下：
乘法噪聲：
I=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
J1=imnoise(I,'speckle',0.03);
subplot(1,3,1);
imshow(J1);
title('d=0.03時的圖像');
J2=imnoise(I,'speckle',0.15);
subplot(1,3,2)
imshow(J2)
title('d=0.15時的圖像');
J3=imnoise(I,'speckle',0.5);
subplot(1,3,3)
imshow(J3)
title('d=0.5時的圖像');
運作結果對比如下：
3.2設計巴特沃斯低通濾波對圖像進行低通濾波處理，顯示結果
巴特沃斯濾波器的特點是通頻帶内的頻率響應曲線最大限度平坦，沒有起伏，而在阻頻帶則逐漸下降為零。 在振幅的對數對角頻率的波得圖上,從某一邊界角頻率開始,振幅随着角頻率的增加而逐漸減少,趨向負無窮大。
3.2.1疊加椒鹽噪聲的巴特沃斯低通濾波
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg'); I=rgb2gray(i);
I1=imnoise(I,'salt & pepper',0.02);
f=double(I1);
g=fft2(f);
g=fftshift(g);
[N1,N2]=size(g);
n=3; %階次設為3
d0=100; %此處d0為截止頻率
n1=fix(N1/2);
n2=fix(N2/2);
for i=1:N1
for j=1:N2
d=sqrt((i-n1)^2+(j-n2)^2);
h=1/(1+0.414*(d/d0)^(2*n));
result(i,j)=h*g(i,j);
end
end
result=ifftshift(result);
subplot(1,2,1),imshow(I1);
title('受椒鹽噪聲污染的圖像');
subplot(1,2,2),imshow(J1);
title('截止頻率為100HZ的巴特沃斯低通濾波處理後');
運作結果：
3.2.2疊加高斯噪聲的巴特沃斯低通濾波
clear;
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg'); I=rgb2gray(i);
I1=imnoise(I,'gaussian',0.02);
f=double(I1);
g=fft2(f);
g=fftshift(g);
[N1,N2]=size(g);
n=3; %階次設為3
d0=100; %此處d0為截止頻率
n1=fix(N1/2);
n2=fix(N2/2);
for i=1:N1
for j=1:N2
d=sqrt((i-n1)^2+(j-n2)^2);
h=1/(1+0.414*(d/d0)^(2*n));
result(i,j)=h*g(i,j);
end
end
result=ifftshift(result);
subplot(1,2,1),imshow(I1);
title('受高斯噪聲污染的圖像');
subplot(1,2,2),imshow(J1);
title('截止頻率為100HZ的巴特沃斯低通濾波處理後');
運作結果：
3.2.3疊加乘性噪聲的巴特沃斯低通濾波
clear;
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg'); I=rgb2gray(i);
I1=imnoise(I,'speckle',0.02);
f=double(I1);
g=fft2(f);
g=fftshift(g);
[N1,N2]=size(g);
n=3; %階次設為3
d0=100; %此處d0為截止頻率
n1=fix(N1/2);
n2=fix(N2/2);
for i=1:N1
for j=1:N2
d=sqrt((i-n1)^2+(j-n2)^2);
h=1/(1+0.414*(d/d0)^(2*n));
result(i,j)=h*g(i,j);
end
end
result=ifftshift(result);
J1=uint8(real(X2));
subplot(1,2,1),imshow(I1);
title('受乘性噪聲污染的圖像');
subplot(1,2,2),imshow(J1);
title('截止頻率為100HZ的巴特沃斯低通濾波處理後');
運作結果：
3.3用MATLAB實作高斯高通濾波器對圖像的處理
I1=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
subplot(1,2,1);
imshow(I1);
title('原圖');
f=double(I1);
g=fft2(f);
g=fftshift(g);
[M,N]=size(g);
d0=100;
m=fix(M/2); n=fix(N/2);
result=g;
for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
subplot(1,2,2);
imshow(J2);
title('高斯高通濾波(d0=100)');
運作結果：
3.4維納濾波和中值濾波對圖像進行處理
随機信号或随機過程(random process)是普遍存在的。一方面，任何确定性信号經過測量後往往就會引入随機性誤差而使該信号随機化；另一方面，任何信号本身都存在随機幹擾，通常把對信号或系統功能起幹擾作用的随機信号稱之為噪聲。噪聲按功率譜密度劃分可以分為白噪聲(white noise)和色噪聲(color noise)，我們把均值為0的白噪聲叫純随機信号(pure random signal)。是以，任何其它随機信号都可看成是純随機信号與确定性信号并存的混合随機信号或簡稱為随機信号。要差別幹擾(interference)和噪聲( noise)兩種事實和兩個概念。非目标信号(nonobjective signal)都可叫幹擾。幹擾可以是确定信号，如國内的50Hz工頻幹擾。幹擾也可以是噪聲，純随機信号(白噪聲)加上一個直流成分(确定性信号)，就成了最簡單的混合随機信号。醫學數字信号處理的目的是要提取包含在随機信号中的确定成分，并探求它與生理、病理過程的關系，為醫學決策提供一定的依據。例如從自發腦電中提取誘發腦電信号，就是把自發腦電看成是幹擾信号，從中提取出需要的資訊成分。是以我們需要尋找一種最佳線性濾波器，當信号和幹擾以及随機噪聲同時輸入該濾波器時，在輸出端能将信号盡可能精确地表現出來。
維納濾波和卡爾曼濾波就是用來解決這樣一類問題的方法：從噪聲中提取出有用的信号。實際上，這種線性濾波方法也被看成是一種估計問題或者線性預測問題。對于在C中的圖檔而言，其進行均值濾波和維納濾波的MATLAB的指令如下：
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
I=rgb2gray(i);
J1=imnoise(I,'gaussian',0,0.1); %高斯噪聲 均值0 方差為0.1
subplot(2,3,4),imshow(K1,[]),title('對高斯噪聲圖像中值濾波');
K2=medfilt2(J2,[3,3]);
subplot(2,3,2),imshow(J2),title('椒鹽噪聲');
subplot(2,3,5),imshow(K2,[]),title('對椒鹽噪聲圖像中值濾波');
K3=medfilt2(J3,[3,3]);
subplot(2,3,3),imshow(J3),title('乘性噪聲');
subplot(2,3,6),imshow(K3,[]),title('對乘性噪聲圖像中值濾波');
運作結果：
以前面經常用的圖檔為例，再進行維納濾波：
代碼如下：
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
I=rgb2gray(i);
subplot(2,4,1),imshow(I); title('原始圖像');
J1=imnoise(I,'gaussian',0,0.005);
subplot(2,4,2),imshow(J1);
title('高斯噪聲圖像');
K1=wiener2(J1);
subplot(2,4,6),imshow(K1);
title('維納濾波');
J2=imnoise(I,'salt & pepper',0.005);
subplot(2,4,3),imshow(J2);
title('椒鹽噪聲圖像');
subplot(2,4,7),imshow(K2);
title('維納濾波');
J3=imnoise(I,'speckle',0.005);
subplot(2,4,4),imshow(J3);
title('乘性噪聲圖像');
K3=wiener2(J3);
subplot(2,4,8),imshow(K3);
title('維納濾波');
運作結果：
參考文獻
[1] 劉泉 《數字信号處理原理與實作》 電子工業出版社
[2] 劉泉 《信号與系統》 高等教育出版社
[3] 周博等編《MATLAB科學計算》 機械工業出版社
[4] 賀興華等著.《MATLAB7.X圖像處理》 人民郵電出版社
[5] 王洪元主編.《MATLAB語言以及在電子資訊工程中的應用》 清華出版社.
[6] 張汗靈. 《MATLAB在圖像進行中的應用》 清華大學出版社
[7] 董霖編 《MATLAB使用詳解——基礎、開發及工程應用》電子工業出版社
[8] Edward W.Kamen.《應用Web和MATLAB的信号與系統基礎》.電子工業出版社