数字信号处理实验报告

 

班    级：  1105074143  

姓    名：  路晓冬    

学    号： 1105074143 

  实验一  频谱分析与采样定理

一、实验目的

1、观察模拟信号经理想采样后的频谱变化关系。

2、验证采样定理，观察欠采样时产生的频谱混叠现象。

3、加深对DFT算法原理和基本性质的理解。

4、熟悉FFT算法原理和FFT的应用。

二、实验原理

    根据采样定理，对给定信号确定采样频率，观察信号的频谱。

三、实验内容和步骤

1）实验内容

      在给定信号为：

1．x(t)=cos(100*π*at)

2．x(t)=exp(-at)

3．x(t)=exp(-at)cos(100*π*at)

其中a为实验者的学号，记录上述各信号的频谱，表明采样条件，分析比较上述信号频谱的区别。

2）实验步骤

1．根据采样理论、DFT的定义、性质和用DFT作谱分析的有关内容。

2．根据FFT算法原理和基本思想。

3．确定实验给定信号的采样频率，编制对采样后信号进行频谱分析的程序

四、实验过程

%实验一:频谱分析与采样定理

T=0.0001;               %采样间隔T=0.0001

F=1/T;                 %采样频率为F=1/T

L=0.02;                %记录长度L=0.02 

N=L/T;                

t=0:T:L;            

a=36;

f1=0:F/N:F;

f2=-F/2:F/N:F/2;

%%%%%%%%%%%%%%%%%%%%%%%%%

x1=cos(100*pi*a*t);

y1=T*abs(fft(x1));  % 求复数实部与虚部的平方和的算术平方根  

y11=fftshift(y1);

figure(1),

subplot(3,1,1),plot(t,x1);title('正弦信号');

subplot(3,1,2),stem(y1);title('正弦信号频谱');

subplot(3,1,3),plot(f2,y11);title('正弦信号频谱');

%%%%%%%%%%%%%%%%%%%%%

x2=exp(-a*t);

y2=T*abs(fft(x2));

y21=fftshift(y2);

figure(2),

subplot(3,1,1),stem(t,x2);title('指数信号');

subplot(3,1,2),stem(f1,y2);title('指数信号频谱');

subplot(3,1,3),plot(f2,y21);title('指数信号频谱');

%%%%%%%%%%%%%%%%%%%%%

x3=x1.*x2;

y3=T*abs(fft(x3));

y31=fftshift(y3);

figure(3),

subplot(3,1,1),stem(t,x3);title('两信号相乘');

subplot(3,1,2),stem(f1,y3);title('两信号相乘频谱');

subplot(3,1,3),plot(f2,y31);title('两信号相乘频谱');

正弦信号频谱：

指数信号频谱：

两信号相乘频谱：

五、实验结果及分析

    奈奎斯特抽样定理为抽样频率必须大于或等于信号频谱最高频率的2倍，即。取采样间隔T=0.0001s，满足奈奎斯特抽样定理。实验中正弦信号的频率分量为，指数信号的频率分量为0，两信号相乘后频率分量为1400Hz，所以三幅频谱图中相应出现1400Hz、0Hz、1400Hz的频率分量，与实际计算结果一致。频谱图中显示区域的最高频率即为在该采样频率下的最高可分析频率，等于信号采样频率的一半，即。

六、实验体会

    通过实验一我熟悉了matlab的使用，加深了对奈奎斯特定理的理解。

实验二  卷积定理

一、实验目的

通过本实验，验证卷积定理，掌握利用DFT和FFT计算线性卷积的方法。

二、 实验原理

时域圆周卷积在频域上相当于两序列DFT的相乘，因而可以采用FFT的算法来计算圆周卷积，当满足时，线性卷积等于圆周卷积，因此可利用FFT计算线性卷积。

三、实验内容和步骤

1、给定离散信号和，用图解法求出两者的线性卷积和圆周卷积；

2、编写程序计算线性卷积和圆周卷积；

3、比较不同列长时的圆周卷积与线性卷积的结果，分析原因。

四、实验过程

%实验二:卷积定理

x=[3 0 2 1 3];  %原始序列

y=[3 0 2 1 3];

%直接计算线性卷积

z=conv(x,y);

figure(1),subplot(311),stem(x);axis([1 9 0 4]);

subplot(312),stem(y);axis([1 9 0 4]);

subplot(313),stem(z);axis([1 9 0 30]);

%利用FFT计算

N=10;%N=10时

x1=[x zeros(1,N-length(x))];

y1=[y zeros(1,N-length(y))];

X1=fft(x1);

Y1=fft(y1);

Z1=X1.*Y1;

z1=ifft(Z1);

figure(2),

subplot(321),stem(x1);

subplot(322),stem(real(X1));

subplot(323),stem(y1);

subplot(324),stem(real(X1));

subplot(325),stem(z1);

subplot(326),stem(real(Z1));

N=5;%N=5时

x2=[x zeros(1,N-length(x))];

y2=[y zeros(1,N-length(y))];

X2=fft(x2);

Y2=fft(y2);

Z2=X2.*Y2;

z2=ifft(Z2);

figure(3),

subplot(321),stem(x2);

subplot(322),stem(real(X2));

subplot(323),stem(y2);

subplot(324),stem(real(X2));

subplot(325),stem(z2);

subplot(326),stem(real(Z2));

直接计算线性卷积：

利用圆周卷积计算，取N=10：

利用圆周卷积计算，取N=5：

五、实验结果及分析

时域圆周卷积在频域上相当于两序列DFT的相乘，因而可以采用FFT的算法

来计算圆周卷积，当满足NM+L—1时，线性卷积等于圆周卷积，因此可利用FFT计算线性卷积。可见，对于直接线性卷积结果，对于圆周卷积，取N=10时，所得结果与直接线性卷积计算结果一致，而当N=5<M+L—1，得出 ，与直接线性卷积计算结果不同，产生混叠。

六、实验体会

    通过实验二我学会了用matlab实现卷积算法的方法，要用FFT计算线性卷积，序列长度必须满足NM+L—1。

实验三  IIR滤波器设计实验

一、实验目的

1.学习模拟－数字变换滤波器的设计方法

2.掌握双线性变换滤波器的设计方法

3.掌握实现数字滤波的具体方法。

二、实验要求 

1、用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π, π]频率区间上，最小衰减大于15dB.

2、0.02π为采样间隔，打印出数字滤波器在频率区间[0, π/2]上的频率响应特性曲线。

3、用所设计的滤波器对实际心电图信号采样序列进行仿真滤波处理，观察总结滤波作用与效果

附：心电图采样序列x(n)

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n),其中存在高频干扰。在实验中以x(n)作为输入序列，滤除其中的干扰成分。

{x(n)}={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0}

三、实验过程

%先设计模拟滤波器，再转化数字滤波器   

wp=0.2*pi;

ws=0.3*pi;

Rp=1;

Rs=15;

Ts=0.02*pi;

Fs=1/Ts;

wp1=2/Ts*tan(wp/2);%将模拟指标转变成数字指标

ws1=2/Ts*tan(ws/2);

[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');  %选择滤波器的最小阶数  

[Z,P,K]=buttap(N);%创建butterworth模拟滤波器

[Bap,Aap]=zp2tf(Z,P,K);

[b,a]=lp2lp(Bap,Aap,Wn);   

[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换 

[H,W]=freqz(bz,az,50);%绘制频率响应曲线

L=length(W)/2+1;

figure(1),plot(W(1:L)/pi,abs(H(1:L))),grid,xlabel('角频率(\pi)'),ylabel('频率响应幅度');

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,...

-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,...

-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

y=filter(bz,az,x);            %滤波

figure(2),

subplot(2,1,1),plot(x),title('原始信号');

subplot(2,1,2),plot(y),title('滤波后信号');

数字滤波器在频率区间[0, π/2]上的频率响应特性曲线：

对实际心电图信号采样序列进行仿真滤波处理：

四、实验结果及分析

用双线性变换法设计巴特沃斯低通IIR数字滤波器，先由已知条件求的数字截止角频率，，转化为模拟截止角频率，，有指标得出，，

继而求出与N与所设计的系统函数，以用双线性变换法将其转换为数字滤波器。

五、实验体会

    通过此实验进一步掌握了用双线性变换法设计巴特沃斯低通IIR数字滤波器的方法，双线性变换法不可能产生频率混叠现象。

实验四  FIR滤波器设计实验

一、实验目的 

1.熟悉滤波器的计算机仿真方法

2.掌握用窗函数法设计FIR数字滤波器的原理和方法。

3.了解各种窗函数对滤波特性的影响

二、实验要求 

1.设计一线性相位FIR低通滤波器滤波器，给定抽样频率为Ω_s=3π×10⁴(rad/s),通带截止频率为Ω_p=3π×10³(rad/s),阻带起始频率为Ω_st=6π×10³(rad/s),阻带衰减比大于50dB。

2.选择不同的窗函数设计该滤波器，观察其频率响应函数有什么变化

三、实验过程

wp1=3000*pi;                  %通带截止角频率  

ws1=6000*pi;                  %阻带截止角频率

wsam1=30000*pi;               %采样角频率

fsam1=wsam1/(2*pi);           %采样截止频率

passrad=(wp1+ws1)/2/fsam1;    %截止频率

A=3.3;

wdelta=(ws1-wp1)/fsam1;

N=ceil(2*pi/wdelta*A);         %滤波器的阶数

w1=hamming(N+1);               %用汉明窗实现

L=N/2+1;

n=1:1:N+1;

hd=sin(passrad*(n-L))./(pi*(n-L));            %理想低通滤波器

if(N==ceil(N/2)*2)

   hd(L)=passrad/pi;

end

h1=hd.*w1';                            %加窗

[mag1,rad]=freqz(h1);

plot(rad,20*log10(abs(mag1)));         %频谱显示

四、实验结果及分析

    阻带衰减比大于50dB，所以选择汉明窗能满足要求。汉明窗时域定义式为，设计所要求的过渡带宽，而汉明窗的带宽为，求出有限长序列长度N，且关于偶对称。理想线性低通滤波器单位抽样响应为，则所设计的滤波器单位抽样响应为。

五、实验体会

通过实验掌握用窗函数法设计FIR数字滤波器的原理和方法，深刻理解了各窗函数的特性。

 

第二篇：MATLAB数字信号处

柔性数学形态学在红外图像边缘

检测中的MATLAB 实现

摘　　　要: 　在介绍适用于灰度图像的柔性形态学概念及特性的基础上,结合红外舰船图像实例,用MATLAB 实现了灰度图像的形态学边缘检测。文章给出了MATLAB 实现的主程序,列出了一些柔性形态学处理和传统边缘检测处理的结果。通过仿真过程和讨论可以看出在图像处理过程中,MATLAB 具有直观、简洁的优越性,可以在这个领域得到广泛的应用。

关　键　词:   红外图像; 边缘检测; 柔性数学形态学

引    言

柔性形态变换是Koskinen 等人在经典形态变换基础上提出的一类非线性算子[1 ]。它放宽了经典形态算子的定义,以获得一定程度的鲁棒性,但还保留了经典形态算子的优良特性,因此在有噪声的情况下目标的边缘检测效果比传统的形态算子性能更好。但在运用柔性数学形态学检测图像边缘时,往往要进行大量的数学计算。目前流行用Fortran、C语言等编制计算程序,既需要对有关算法有深刻的了解,还需要熟练地掌握所用语言的语法及编程技巧。对多数科学工作者而言,同时具备这两方面才能有一定困难。MATLAB 语言被称为是一种“演草

纸式的科学计算语言”[2 ],它的强大计算和模拟功能使得在许多应用领域的各种计算、演算、模拟等的工作变得相当简单,是一份实时进行图像处理的有力工具。本文着重讨论了通过MATLAB 运用柔性数学形态学来进行红外图像处理的实现问题,结合红外舰船图像实例用MATLAB 进行了柔性形态学边缘检测和传统算子的边缘检测,给出了MATLAB 实现的

主程序,比较了这些算法,讨论了MATLAB 在这方面的优越性和实际应用。

1 　柔性形态学的基本原理和运算

数学形态学是一门建立在积分几何与随机集论基础上的结构分析方法,积分几何能够得到各种几何参数的间接测量,以及反映图形的体视性质从而反映出图形的结构、形态的性质,而图像中的随机性质(如灰度分布等) 的处理则以随机集论为基础。数学形态学之所以能在图像处理中得到广泛应用,在于它的研究对象是n 维空间里目标的形状,而不是数学指标;它通过对目标的形状变换,能有效地提取出与目标有关的几何信息。数学形态学的主要内容是借助于不同形状的结构元素与图像间的一系列结构变换来处理和分析图像。从广义上看,数学形态学是图像处理的统一理论,因为传统的图像处理中的线性算子和非线性算子都是形态学算子的特例。柔性形态学方法[3 ]是将“顺序统计”的思想注入标准数学形态学中的一种方法,是一种顺序滤波。柔性数学形态学方法用排序加权统计方法代替最小、最大法。权值与结构元素有关,并由核心和软边界两大部分组成。柔性数学形态学具有硬数学形态学相似的代数特性,但具有更强的抗噪声干扰的能力,对加性噪声及微小形状变化不敏感。

1. 1 　柔性形态学算子

在柔性形态变换中,结构元素被分割成“硬核”(相当于标准的结构元素) 和“柔性边缘”两部分,而经典形态算子中的最大最小运算,在柔性形态算子中被排序统计所代替。为了便于用计算机实现形态学运算,下面给出灰度图像形态学基本运算的一种定义方法。设F 表示原始红外灰值图像, F( m , n) 表示图像F 在( m , n) 点的灰度。设集合A , B 为定义在Z2

上的凸集,且使A A B , B 被分为“硬核”A 和“柔性边缘”B \ A 两个子集, 这里“ \ ”代表集合差。A , B均为结构元素。柔性形态学: 首先定义一种重复集(multiset) ,重复集中包含的元素可以重复。元素f ( a) 重复k 次被表示为

{ k ◇f ( a) } = { f ( a) , f ( a) , ?, f ( a) } ( k 次)(1)

式中k 为正整数, 且1 ≤ k ≤ min{ card ( B) / 2 ,

card ( B \ A) }

其中card( B) 代表集合B 的基数(cardinality) ; a ∈B ;

柔性形态腐蚀:

F ß [ B , A , k ] ( x) = ({ k ◇f ( a) | a ∈Ax} ∪

{ f ( b) | b ∈( B \ A) x} ) 中第k 个小的值。

柔性形态膨胀:

F Ý [ B , A , k ] ( x) = ({ k ◇f ( a) | a ∈Ax} ∪

{ f ( b) | b ∈( B \ A) x} ) 中第k 个大的值。

　A 、B 的映像集分别为^A ( x) = { x | x = - a , a

∈A} 和^B ( x) = { x | x = - b , b ∈B} ,则形态学开

闭算子定义如下:

形态学开:

F . [ B , A , k ] = ( F ß [ B , A , k ]) Ý [ ^B , ^A , k ](2)

形态学闭:

F ·[ B , A , k ] = ( F Ý [ B , A , k ]) ß [ ^B , ^A , k ](3)

柔性形态滤波是极值滤波的推广,当K = 1 时,即为极小值滤波;当K =| B | 时,即为极大值滤波,当K 为奇数时, 即为中值滤波。通过采用表示结构元素与图像匹配程度的结构元素参数K, 可形成加权顺序滤波。而常用的加权顺序滤波采用“中心加权”,即对具有对称结构的结构元素, 中心的权高于或低于四周的权。柔性形态的腐蚀和膨胀运算与原图的差可以用

作全方位的边缘检测。用这种结构元素作开运算,可以去掉宽度小于结构矩阵直径的“毛刺”或噪声点;作闭运算可以去掉宽度小于结构矩阵直径的“麻坑”,但可以保留“毛刺”和噪声点。首先用柔性形态开运算抑制噪声,然后利用膨胀与原图的差来检测

图像的边缘。

1. 2 　结构元素

数学形态学对图像的分析处理主要依据于特定的形状变换,而形状变换的特性由具体的运算和结构元素的几何特征决定。并且,形态变换的元素复杂度直接依据于结构元素的形状、大小。因此,运用数学形态学进行图像处理的关键在于结构元素的选择。可以说,数学形态学最大的贡献是在几何变换中引入了“形态滤波器”———结构元素的概念。结构元素具有各向同性,且它的大小代表不同尺度的处理效果。不同的结构元素对同一图像进行相同的形态运算,会导致图像处理的结果有所差异。因此,结构元素并没有固定的形状大小,必须根据不同的目标图像和所需信息的形状特征,设计、选择具体的结构元素。结构元素图像可以是二值图像也可以是灰度图像,当结构元素图像为二值图像时,其所描述的只有目标的形状信息,而无法描述目标的灰度分布特征,从而相应的灰度形态学处理相当于在二值化的灰度图像上进行二值形态操作。当结构元素图像为灰度__图像时,则可利用结构元图像中所包含的灰度信息去获取二值形态操作无法得到的目标的灰度分布特征,达到更好的图像处理效果。由于处理的是海空背景下的舰船红外图像,经过试验后结构元素取为十字形,而且为灰度值图像。

2 　柔性形态学的MATLAB 实现和计

算机仿真结果

　　MATLAB 的图像处理工具箱中含有形态学函数,如:dilate. m、erode. m 和bwmorph. m 等。但这些函数只是针对二值图像处理的,其中的结构元素也是二值图像。实际应用中绝大多数遇到的是灰值图像,如红外图像、遥感图像、显微图像、照片等。一般情况下,可以选取适当阈值将灰度值图像转换为二值图像,再用相应的运算进行处理,但是这种转换势必会损失一些有用的信息,况且由二值图像再恢复为灰值图像又是一件极为困难的工作。为此,在原

有库函数的基础上自行构建灰度图像形态学函数,并加入顺序统计因子K构成柔性形态学函数来处理舰船红外图像。

2. 1 　灰度柔性形态学腐蚀算子的MATLAB程序

clear all ;

F1 = imread(’test256. bmp’) ;

F = double (F1) ;

[ a ,b ] = size (F) ;

SE = [1 1 1 1 1 ;1 1 1 1 1 ;1 1 1 1 1 ] 3 180

[ c ,d ] = size (SE) ;

e = (c + 1) / 2 ;

f = (d + 1) / 2 ;

hardcore = [1 3 ;2 1 ;2 2 ;2 3 ;2 4 ;2 5 ;3 3 ] ;

k = 2 ;

R = zeros (a ,b) ;

for x = 1 :a

　for y = 1 :b

　　sort number = 0 ;

　　sort numberh = 0 ;

　　for i = 1 :c

　　　m = i - e ;

　　　for j = 1 :d

　　　　n = j - f ;

　　　　if (x + m) > = 1 & (x + m) < = a & (y

+ n) > = 1 &(y + n) < = b

sort number = sort number + 1 ;

D(sort number) = F ( (x + m) , (y + n) ) + SE(i ,

j) ;

　if all ( [ i ,j ] = = hardcore (1 , :) ) | all ( [ i ,j ] = =

hardcore (2 , :) )

sort numberh = sort numberh + 1 ;

H(sort numberh) = D(sort number) ;

　　　　　end

　　　　end

　　　end % now the difference set is completed

　　　end

　　　　number calculation = max(size (D) ) ;

　　　　number k = max(size (H) ) ;

T = zeros (1 ,number calculation + (k - 1) 3 num2

ber k) ;

　　　　for repeat number = 1 :k - 1

T( (repeat number - 1) 3 number k + 1 :repeat

number 3 number k) = H;

　　　　　　end

　T( (k - 1) 3 number k + 1 :max ( size (T) ) ) =

D ;

　　T = sort (T) ; %sort the difference set

　　R (x ,y) = T(max ( size ( T) ) - k + 1) ; %the

kth smallest number

　　if R(x ,y) > 255

　　　R(x ,y) = 255 ;

　　end

　　D = 0 ; %2003. 6. 2

　　H = 0 ;

　end

end

上面给出的是柔性形态学腐蚀运算的MATLAB主程序,柔性膨胀程序可以相应给出。由这两种最基本的形态学算子进行不同组合就可以组合成许多形态算子,如开、闭运算等。程序中SE 为结构元素,取为灰度值结构,值的大小由阈值理论确定。传统的边缘检测方法如Robert 算子、Prewitt 算子和Sobel 算子等在图像处理工具箱中都有现成的函数可以调用。调用格式为edge (image ,’Robert’) 、edge (image ,’Prewitt’) 和edge (image ,’sobel’) 。

2. 2 　舰船红外图像的MATLAB仿真结果及讨论

对实际的海天背景下的舰船红外图像使用柔性形态学进行的MATLAB 仿真结果图如下。在试验中k 分别取1 和2 ,结构元素取3 ×3 和5 ×5 的矩阵,内核为十字形。图1 为舰船红外原始图像,大小为256 ×256。图2 是使用3 ×3 的结构元素,内核是十字型, k 取2的柔性形态膨胀运算结果图。图3 是图1 的柔性开运算图, k 为2。比较图1 ,图2 和图3 ,可以看出形态膨胀的结果是强化了亮噪声,提升了整体图像的灰度;而柔性形态开运算则有效地抑制了噪声,而保留了大的灰度区域。图4 是使用5 ×5 的结构元素,内核为十字型, k 取1 的形态开运算图。与图3 比,噪声抑制得更厉害,原来海面亮斑区域出现一些暗点。图5 是柔性算子处理得到的图像边缘检测图,k为2。图6～图8 分别是Robert 、Prewitt 和Sobel 形态算子的检测图。从这4 个处理结果可以明显看出图5 的结果是最好的,满足了图像处理的要求,并且在边缘的连续性及各向同性方面都优于传统方法。上面的仿真过程都是通过MATLAB 软件实现的,可以看出MATLAB 在验证算法,实现一些复杂图像处理方面具有别的软件所无法比拟的简洁性、方便性和快速性。通过MATLAB 仿真可以方便地进行算法结果讨论,大大加快研究进程。

3 　结论

由以上讨论可知,这种柔性数学形态学在红外图像边缘检测中具有明显的优越性。整个算法的验证仿真可以通过简单的MATLAB 命令来完成,而要完成同样的任务其他高级计算机语言需要使用许多编码才能完成。这样在实际的图像处理过程中使用MATLAB 可以大大提高试验的效率,快速验证研究中的新设想,这些性能对于这个领域的研究是有实际的意义的。

参考文献:

[1 ] 　KOSKINEN L , ASTOLA J , NEUVO Y. Soft morphological

filters[J ] . Proc. SPIE Int . Society of Optical Engineering

[C] . SPIE Vol 1568 , 1991 :2622270.

[2 ] 　PART2EANADER E1SJOBERGA. MATLAB5 手册[M] . 北

京:机械工业出版社,20001

[3 ] 　唐常青,吕宏伯1 数学形态学方法及其应用[M] . 北

京:科学出版社,19901

[4 ] 　DONG Y Z, ZHOU X D , SHEN T S , LOU S L. Edge De2

tection of IR Ship Images Based on Soft Mathematical Mor2

phology[A] . International Symposium on Optical Science and

Technology[C] . SPIE , Vol : 5203 , 20## :5812590.

[5 ] 　DONG Y Z,WANGCJ , ZHOU X D. Structuring Element in

Detection of IR Images Based on Soft Mathematical Morpholo2

gy[A] . International Symposium on Optical Science and Tech2

nology [C] . SPIE , Vol : 5203 , 20## :57125811

图1 　舰船红外图像　　     　图2 　k = 2 形态膨胀图

图3 　k = 2 形态开运算　　　    图4 　k = 1 形态开运算

图5 　柔性形态算子　　　      图6 　Robert 形态算子

图7 　Prewitt 形态算子　　　      图8 　Sobel 形态算子