数字图像处理_平均滤波与中值滤波(含MATLAB代码)[小编推荐]

第一篇：数字图像处理_平均滤波与中值滤波(含MATLAB代码)[小编推荐]

数字图像处理实验二

15生医

一、实验内容

产生教材104页题图4.18（右图）所示的二值图像（白为1，黑为0），编程实现习题4.18所要求的处理（3x3的平均滤波和中值滤波）功能（图像四周边界不考虑，处理结果按四舍五入仍取0或1），显示处理前后的图像，比较其异同。

二、运行环境 MATLAB R2014a

三、运行结果及分析

1.四种不同的窗的3x3平均滤波

4邻域平均滤波后图像8邻域平均滤波后图像4邻域加权平均滤波后图像8邻域加权平均滤波后图像

①在MATLAB图形窗界面进行放大可以看出四者之间的差别： 4领域与8邻域之间没有明显差别，但是加权与未加权之间的差别较为明显，体现在：加权后每个矩形块的四个尖角部分都被保留了下来（图像四周边界不考虑），而未加权的尖角处黑色变为白色。②原因分析：

加权后尖角处原来白色的点（1）进行计算3/5=0.6四舍五入后值为1，保持白色，原来黑色的点（0）进行计算2/5=0.4四舍五入后值为0，保持黑色；而未加权尖角处无论原来是黑色还是白色，进行计算2/4=0.5四舍五入后值为1，所以原先的黑色（0）也变成了白色（1）。③下图为放大后的截图：

2.中值滤波与原图像的对比 原图像中值滤波后图像

①在MATLAB图形窗界面进行放大后可观察出：

使用3x3方形中值滤波模板的效果与4领域、8领域加权平均滤波的效果相同，每个矩形块的四个尖角部分都被保留了下来（图像四周边界不考虑）。②原因分析：

套用3x3方形中值滤波模板后，尖角处原来白色的点（1）在窗内1多于0，取中值后仍保持白色，原来黑色的点（0）在窗内0多于1，取中值后仍保持白色。③下图为放大后的截图：

四、心得体会

通过MATLAB编程更加理解了课后题的计算结果，直观地看出了黑白像素点灰度值变化前后的取值异同。同时，对MATLAB实现像素点灰度值的替换机理也有所掌握，比如后边附的程序中会提到的“%”标注的思考。

五、具体程序

% 生成黑白块图像

unit=zeros(64,64);f=zeros(256,256);for i=1:1:32 for j=1:1:32 unit(i,j)=1;% 1/4白块 end end for i=33:1:64 for j=33:1:64 unit(i,j)=1;% 1/4白块 end end for i=1:64:256 for j=1:64:256 f(i:i+63,j:j+63)=unit;end end

% 对原图像进行3x3平均滤波 for i=2:1:255 for j=2:1:255 fave4(i,j)=(f(i-1,j)+f(i+1,j)+f(i,j-1)+f(i,j+1))/4;

fave8(i,j)=(f(i-1,j-1)+f(i-1,j)+f(i-1,j+1)+f(i,j-1)+f(i,j+1)+f(i+1,j-1)+f(i+1,j)+f(i+1,j+1))/8;

fave4_weighted(i,j)=(f(i-1,j)+f(i+1,j)+f(i,j-1)+f(i,j+1)+f(i,j))/5;

fave8_weighted(i,j)=(f(i-1,j-1)+f(i-1,j)+f(i-1,j+1)+f(i,j-1)+f(i,j+1)+f(i+1,j-1)+f(i+1,j)+f(i+1,j+1)+f(i,j))/9;end end fave4=round(fave4);%平均后灰度值有可能是小数，要取整 fave8=round(fave8);fave4_weighted =round(fave4_weighted);fave8_weighted =round(fave8_weighted);subplot(2,2,1);imshow(fave4);title('4邻域平均滤波后图像');subplot(2,2,2);imshow(fave8);title('8邻域平均滤波后图像')subplot(2,2,3);imshow(fave4_weighted);title('4邻域加权平均滤波后图像')subplot(2,2,4);imshow(fave8_weighted);title('8邻域加权平均滤波后图像')

4邻域平均滤波后图像8邻域平均滤波后图像4邻域加权平均滤波后图像8邻域加权平均滤波后图像

% 对原图像进行3x3方形中值滤波 for i=2:1:255 for j=2:1:255

a=[f(i-1,j-1),f(i-1,j),f(i-1,j+1),f(i,j-1),f(i,j),f(i,j+1),f(i+1,j-1),f(i+1,j),f(i+1,j+1)];b=sort(a);% 排序函数

fmid(i,j)=b(5);% 9个数排序的中值 end end subplot(1,2,1);imshow(f);title('原图像');subplot(1,2,2);imshow(fmid);title('中值滤波后图像')

原图像中值滤波后图像

第二篇：数字图像处理_领域平均滤波_中值滤波

东华大学实验报告

课程 数字图像处理 名称 数字图像变换

实验名称: 邻域平均法(box模板)和中值滤波处理

一、实验目的

图像变换是数字图像处理中的一种综合变换，如直方图变换、几何变换等。通过本实验，使得学生掌握两种变换的程序实现方法。

二、实验任务

请设计程序，分别用邻域平均法，其模板为：

111

111*1 9111 和中值滤波法对testnoise图像进行去噪处理（中值滤波的模板的大小也设为3×3）。

三、实验环境

本实验在Windows平台上进行，对内存及cpu主频无特别要求，使用VC或者MINGW（gcc）编译器均可。

四、设计思路

介绍代码的框架结构、所用的数据结构、各个类的介绍（类的功能、类中方法的功能、类的成员变量的作用）、各方法间的关系

试验要求中以给出大致的编程思路和源代码以及代码注释，只有黑框部分需要自己填写。在此不进行赘述。

五、具体实现

实现设计思路中定义的所有的数据类型，对每个操作给出实际算法。对主程序和其他模块也都需要写出实际算法。

注意：源代码中要加上注释。代码：(红色为重点代码)<邻域平均法>(3*3)/*------利用第一次实验课提供的 dhc.h 和 dhc.c文件以获取位图的高 宽 以及从文件头到实际的位图数据的偏移字节数，从而实现对位图实际数据的操作。------*/ #include #include #include #include “hdr.h” /*------定义结构指针------*/ struct bmphdr *hdr;//定义用于直方图变量

unsigned char *bitmap,*count,*new_color;/*------main()函数编写------*/ int main(){ //定义整数 i，j 用于函数循环时的，nr_pixels为图像中像素的个数

int i, j ,nr_pixels,nr_w,nr_h;//定义两个文件指针分别用于提取原图像的数据和生成直方图均衡化后的图像

FILE *fp, *fpnew;//定义主函数的参数包括：输入的位图文件名和输出的位图文件名,此处内容可以不要，在DOS下执行命令的时候再临时输入也可，为了方便演示，我这里直接把函数的参数确定了。

// argc=3;// argv[1]=“test.bmp”;// argv[2]=“testzf.bmp”;

//参数输入出错显示

/* if(argc!= 3){ printf(“please input the name of input and out bitmap filesn”);exit(1);}*/ // 获取位图文件相关信息

// hdr = get_header(argv[1]);hdr = get_header(“testnoise.bmp”);if(!hdr)exit(1);//以二进制可读方式打开输入位图文件

fp = fopen(“testnoise.bmp”, “rb”);if(!fp){

printf(“File open error!n”);

exit(1);} // 文件指针指向数据区域

fseek(fp, hdr->offset, SEEK_SET);

//计算位图像素的个数 nr_pixels = hdr->width * hdr->height;nr_w = hdr->width;nr_h = hdr->height;bitmap = malloc(nr_pixels);new_color = malloc(nr_pixels);count = malloc((nr_w+2)*(+nr_h+2));//读取位图数据到bitmap中

fread(bitmap, nr_pixels, 1, fp);fclose(fp);

//因为图像边缘无法使用邻域平均，所以根据邻近颜色填补图像的周围一圈，存入count[]数组中

//中心图像存入count[] for(i=nr_w+3;i<(nr_w+2)*(+nr_h+2)-nr_w-3;i++){

j=i/(nr_w+2);

if(i%(nr_w+2)!=0&&(i+1)%(nr_w+2)!=0)count[i]=bitmap[i-nr_w-1-2*j];} //填补第一排

for(i=1;i

count[i]=bitmap[i-1];} //填补最后一排

for(i=1;i

count[(nr_w+2)*(nr_h+1)+i]=bitmap[nr_w*(nr_h-1)+i-1];} //填补左边一排

for(i=0;i

count[i*(nr_w+2)]=count[i*(nr_w+2)+1];} //填补右边一排

for(i=0;i

count[(i+1)*(nr_w+2)-1]=count[(i+1)*(nr_w+2)-2];}

//邻域平均3*3 for(j=nr_w+3,i=0;j<(nr_w+2)*(+nr_h+2)-nr_w-3;j++){

if(j%(nr_w+2)!=0&&(j+1)%(nr_w+2)!=0)

new_color[i]=(count[j]+count[j-1]+count[j+1]+count[j-nr_w-2]+count[j-1-nr_w-2]+count[j+1-nr_w-2]+count[j+nr_w+2]+count[j-1+nr_w+2]+ count[j+1+nr_w+2])/9,i++;}

//结果存入bitmap[]中

for(i = 0;i < nr_pixels;i++)

bitmap[i]=new_color[i];

// 打开一个以输出文件名命名的文件,设为可写的二进制形式

fpnew = fopen(“test_lynoise.bmp”, “wb+”);

//由于位图文件的头部信息并没有因直方图均衡化而改变，因此输出图像的头部信息从原位图文件中拷贝即可：

fwrite(hdr->signature, 2, 1, fpnew);fwrite(&hdr->size, 4, 1, fpnew);fwrite(hdr->reserved, 4, 1, fpnew);fwrite(&hdr->offset, 4, 1, fpnew);fwrite(&hdr->hdr_size, 4, 1, fpnew);fwrite(&hdr->width, 4, 1, fpnew);fwrite(&hdr->height, 4, 1, fpnew);fwrite(&hdr->nr_planes, 2, 1, fpnew);fwrite(&hdr->bits_per_pixel, 2, 1, fpnew);fwrite(&hdr->compress_type, 4, 1, fpnew);fwrite(&hdr->data_size, 4, 1, fpnew);fwrite(&hdr->resol_hori, 4, 1, fpnew);fwrite(&hdr->resol_vert, 4, 1, fpnew);fwrite(&hdr->nr_colors, 4, 1, fpnew);fwrite(&hdr->important_color, 4, 1, fpnew);if(hdr->offset > 54)

fwrite(hdr->info,(hdr->offset54), 1, fpnew);

//直方图均衡化的数据(bitmap)赋值 fwrite(bitmap, nr_pixels, 1, fpnew);//关闭

fclose(fpnew);//释放内存（优化程序必需）free(hdr);free(bitmap);free(new_color);free(count);return 0;

得出实验结果图像后，比较这两种方法去噪的效果好坏，并分析具体原因。

通过比较，邻域平均法在降低噪声的同时，会使图像产生模糊；中值滤波在降低噪声的同时，能够较好的保持图像边缘，不过我感觉有点水粉画的效果，让图像有点失真。

完成上述工作后，使用程序进行验证分析：使用邻域平均法时，3×3和5×5模板大小对图像进行处理的效果有何差别？并分析原因。

程序及源代码都放在文件夹中，不放进报告了。上面是对比图，显然5×5模板比3×3模版模糊，因为使用邻域平均法时，模版尺寸越大，则图像模糊程度越大。

六、心得体会

这次实验是对两种去噪方法的比较。而在书本中，我们已经看过两种去噪方法的效果了，通过上机验证，果然如此。

在编写代码中，我把算法设计的比较繁琐，因为图像每一个像素的信息是由一维数组bitmap记录的，所以在算均值与中值时花了很大的力气。主要是添补周围一圈的图像比较繁琐，还有就是在一维数组中找模版中心和模版中的其他点。虽然我觉得用二维数组会使算法比较简单，但是我忘了二维数组如何动态分配空间，所以只能硬着头皮用一维数组做。

第三篇：浅谈 中值滤波

浅谈中值滤波

1.中值滤波的现状

在数字信号处理和数字图像处理的早期研究中,线性滤波是主要的处理手段。线性滤波简单的数学表达式以及某些理想特性使其很容易设计和实现。然而，当信号中含有非叠加性噪声时,例如非线性引起的噪声或非高斯噪声等,线性滤波的处理效果就很难令人满意。在处理图像时,线性滤波将破坏边缘,而且不能有效滤除脉冲噪声。为了克服线性滤波方法的局限性,研究非线性滤波的方法为数字信号处理重要课题之一。

非线性滤波基于对输入信号序列的一种非线性映射关系,常可把某一特定的噪声近似映射为零而保留信号的重要特征,因而可以在一定程度上克服线性滤波的不足。

1971年著名学者J.W.Tukey在他的开拓性论文中提出了中值滤波的概念并用作时间序列平滑。中值滤波一出现就因其具有对尖脉冲的良好抑制能力，在平滑加性噪声时能保持信号的边缘特征等优点而备受瞩目。

常用的中值滤波是非线性滤波的代表，由于经典的中值滤波算法在滤除噪声的同时会使信号中重要的细节信息受损,因此,许多改进的中值滤波算法相继被提出。2.中值滤波

在数据处理中我们经常使用的是滑动中值滤波，即取定中值滤波的跨度N(一般 N 为奇数)，在数据序列中顺次取得 N 个数据，然后将该数据列的中值作为中心位置的值输出以形成新的数据序列，在滤波中应将原数据序列的两个边界各补充(N-1)/2(N为奇数时)个等于边界的点以使滤波后的新数据序列长度与原始的数据序列长度一致。2.1一般中值滤波

2.1．1一般中值滤波的基本原理

设有一个序列：x1,x2,x3,x4,x5，将它们按照绝对值大小重新排列此序列

x3, x5,x2, x4,x1

重排以后的中值是x2，此值就作为滤波的输出。显然，x2不能表示成输入数据和滤波系数的褶积的线性组合。其主要特点有：(1)一般中值滤波绝对阻止噪声峰值，因为中值滤波只取中位数，绝对不会取异常数。例如有一组数(x1,x2,x3,x4,x5)正常数 − a≤xn≤a，n=2,3,4,5 异常数 x1 >>a a 表示一个数，将以上数组自小到大排列后为(x3, x5,x2, x4,x1)取中位数x2，决不会取异常数x1。

(2)一般中值滤波是低通滤波器，中值滤波取中值为序列的输出，可以看作是对数据序列进行局部平滑，这种局部平滑实质就是低通滤波。

(3)一般中值滤波不改变阶越函数在空间、时间上的位置，这一性质对于信号处理中的保护边缘有着重要的作用。

(4)当中值滤波的滤波窗口足够长时，有限宽度的三角波和矩形波可以被完全平滑。

(5)中值滤波由于没有统计效应，对随机出现的小的振幅值有时不能完全平滑，所以通常信号在中值滤波处理以后需要再进行带通滤波。2.1.2一般中值滤波(MF)的数学基础

中值滤波对数字序列有平滑作用，平滑也就是数据逼近，这样则存在误差，如何利用误差最小来确定平滑参数，一般常见的有两种准则：(1)使误差的平方和达到最小；(2)使误差的绝对值和达到最小。平均值平滑的数学原理应用准则(1)，即符合误差的平方和最小。中值滤波则是利用准则(2)来实现对数据序列的平滑。

设x是 n 个数据序列的中位数，xi 表示一组序列。x与xi之差的绝对值和为：

Qxxi（3.1.1）

i1n要使Q最小，则

Q0（3.1.2）xnnxxiQn2即 (xxi)Sign(xxi)0 xxi1i1xxii1式中：Sign——符号函数。

当xi >x 时，Sign 为负； 当xi

这样在选择x时，使得在 n 个数中，有 n/2 个xi大于x,同样有 n/2 个xi小于x，中间的xi即为x；如果 n 为偶数，则取中间的两个xi的平均值为x。2.2 加权中值滤波

2.2.1 加权中值滤波(WM)的基本原理

由上可以看出通过改变加权系数，完全可以改变中值滤波的性质，来达到我们的要求。

2.3 一维中值滤波对信号作用的结果分析

由于中值滤波是一种特殊的非线性滤波手段，它对脉冲的响应为零，（在一个输入上施加一个脉冲函数引起的时间响应。）所以在傅氏域没有“真正的”振幅谱和相位谱。我们只能通过它对已知信号及其频谱特征的响应来分析其各种滤波特性。虽然中值滤波的理论比较完善，但是由于多数情况处理的信号是对称信号，所以并没有人注意到中值滤波对信号相位的影响。

2.3.1一般中值滤波对对称信号相位的影响

（1）在频谱图中，一般中值滤波引入了假高频成分，并且在子波的频带范围内，滤波后子波的主要频带向低频方向移动，此特点在数据处理时应该着重注意，要根据数据处理时的具体要求来判断，同时也成为选择滤波长度的一个条件。

（2）经过一般中值滤波后对称信号的相位不发生移动，这使得我们在处理由对称信号(例如雷克子波、奥姆斯比子波等)作为子波的合成地震记录时，不需要考虑相移问题。但由此就得出结论说中值滤波处理后的所有类型的信号的相位都不发生移动则是片面甚至错误的。2.3.2 一般中值滤波对非对称信号相位的影响

一般中值滤波对非对称信号的处理效果不同于处理对称信号，如果 用处理对称信号的规律来对待非对称信号则往往不能达到预期的效果。对应滤波后的频谱同样向低频方向移动，但假高频现象却并不如对称信号滤波后明显。处理非对称信号的同时必须注意选择的滤波点数是否使相位的改变在要求的范围内。从滤波后和滤波前的最大振幅平方比来看滤波前后的能量变化，发现在同等情况下，一般中值滤波对非对称信号的衰减能力大于对对称信号的衰减。2.3.3 加权中值滤波对对称信号相位的影响

（1）在频谱图中，加权中值滤波也引入了假高频成分；并且滤波后的子波的主要频带向低频方向移动，说明了加权中值滤波的低通滤波特性。

(2)同样加权中值滤波对对称信号的相位不产生影响。2.3.4 加权中值滤波对非对称信号相位的影响

当滤波长度大于一定的子波宽度时，波形已经失去了原有的形态，但是在波形失去原有形态之前，经过加权中值滤波后的子波表现出较好的分辨率特性；在频谱上加权中值滤波仍然表现出低通特性；信号的相位也因滤波而产生了畸变。

2.3.5 一般、加权中值滤波对不同信号作用的比较 一般中值滤波和加权中值滤波对于同一种信号表现出相似的特性：二者在处理对称信号时，都起到了衰减的作用，并且对信号的相位都不产生影响，同时使信号的频谱中掺入了假高频成分，还表现出了中值滤波的低通特性；在处理非对称信号时，除了对信号产生衰减作用外，还使信号的相位发生了畸变。

尽管一般和加权中值滤波有相似之处，但是它们还是存在着较大的差异：

在处理对称信号时，一般和加权中值滤波分别对同一信号进行滤波以后，信号的峰值很接近，但是加权中值滤波比一般中值滤波更有利于提高信号的分辨率，在频谱上加权中值滤波比一般中值滤波表现出更加严重的假高频现象。

在处理非对称信号时，加权中值滤波比一般中值滤波表现出的更好的提高分辨率的性质。而在相位谱分析中，尽管加权和一般中值滤波都使信号的相位发生改变，但是在滤波长度较小的情况下，经过一般中值滤波得到的信号的相位曲线虽然已经发生改变，但是仍然与原始信号的相位曲线有相同的趋势，并没有偏离太多；而此时即使在滤波长度较小的情况下，经过加权中值滤波得到的信号的相位曲线已经变得不可辨认了。

经过以上的讨论，我们可以认识到在实际地震资料的处理中，应用中值滤波除了应该考虑信噪比和分辨率以外，更加不容忽视的就是信号的相移问题，这对于资料的可信度起着至关重要的作用。由于在实际中，经常用到的一维中值滤波是不加权的，所以常常把一维不加权中值滤波简说成一般中值滤波，但是随着对处理手段的进一步要求，加权中值滤波的地位日益突出，并且毕竟不加权的情况只是一种特殊的加权中值滤波，所以一般中值滤波的概念也应该扩充为加权中值滤波。

通过对两种不同加权中值滤波（一般、加实数权）的讨论，总结出了一些关于一维中值滤波方面的经验：

(1)通过不同权系数的选取，中值滤波表现出不同的特性，我们可以根据对实际情况的分析来选取不同的权系数以适应各自的需要。

(2)本次只是选取了两种特殊的权系数来分析，而在实际中存在着更多的权系数的选取方法，但是不管权系数的形式如何，都可以仿照本文的方法加以研究。

(3)虽然中值滤波可以满足一定的要求，但是我们同时也应该注意到它们存在的问题：①中值滤波会引起信号形态上的畸变，而且畸变程度和滤波长度有关；②中值滤波会引入假高频，因此信号在经过中值滤波后可以根据情况做一次低通滤波；③中值滤波对非对称信号进行处理时，会引起相位畸变，因此在使用中值滤波之前应该试验相位畸变是否在处理的允许的范围内；④虽然选取适当的权系数后，加权中值滤波可以使信号提高分辨率，但是同时带来“小台阶”效应，因此经过加权中值滤波处理后的信号推荐做一次平滑处理。

由于中值滤波是一类特殊的滤波方法，因此我们利用它进行信号处理时应该格外注意。为了得到预期的效果，处理之前做一下试验以确定最佳的滤波长度是非常必要的。3中值滤波在地震资料处理中的应用 3.1中值滤波在井间地震资料处理中的应用

中值滤波是一种简便有效且信号失真较小的信号处理方法。在不同的道集域下，井间地震资料中的直达波、一次反射波和多次反射波在相邻道间的时差具有不同的表现形式，利用这一特点，应用中值滤波在不同道集域内对井间地震资料进行滤波处理，可以得到很好的效果。

对于井间地震资料，我们所需要的有效反射波是来自于激发点与接收点下方的一次反射波(上行反射)和来自于激发点与接收点上方的一次反射波(下行反射)，其它波均视为相关干扰或无效信息。中值滤波是以正常时差不同为基础的多道滤波技术，在井间地震特殊的观测系统中，中值滤波可以发挥其自身的优点。

通过对井间地震不同道集域下道间时差的分析可以知道，仅运用中值滤波即可达到较好的波场分离效果。为了验证不同道集域下中值滤波对数据处理的效果，进行直达波与多次波的衰减、一次反射波的增强以及上下行波场的分离。3.1.1直达波和多次波的衰减

首先对数据进行带通滤波，消除有效频带之外的噪音干扰，将共炮点道集重排为共偏移距道集。在共偏移距道集下，根据(1)式和(3)式可知，直达波和多次波除了受速度影响外，其相邻道间时差为0，通过共偏移距道集对初至时间拉平排齐，在一定程度上消除了速度的影响，然后选择适当的时窗参数，采用中值滤波消除相邻道时差为0的波组记录，使直达波和多次反射波得到衰减。3.1.2反射波的增强

把衰减了直达波和多次波的数据体重新抽道组成共中心点道集，对于共中心点道集，由于△s=-△g，根据(2)式可知，一次反射波在不受速度影响的情况下其相邻道间的时差为0，通过共中心点道集对反射波时间拉平排齐，消除速度的影响，再次做中值滤波处理，本次中值滤波是为了保留相邻道时差为0的波组记录，而相邻道时差不为0的波组记录将被减弱，因而一次反射波同相轴得以增强，而其它 波场(如直达波和多次波)再次得到衰减(图4，虚线圈)。3.1.3上、下行波场的分离

首先对上行反射波进行拉平(图5中的②)，然后通过中值滤波使上行反射增强而下行反射减弱(图5中的③)，最后返回原始时间剖面得到上行反射波场(图5中的④)；反之，得到下行反射波场(图5中的⑤)。图6为通过中值滤波最终获得的上行和下行反射波场。

在波场分离中，对于资料相对较好的地震数据，仅应用中值滤波即可达到较好的波场分离效果；对于信噪比较低的资料，可以用中值滤波技术使资料的有效波场加强，并得到上、下行反射信息，然后再利用中值滤波进行波场分离。

井间观测系统所具有的特殊性，使得同一种地震波在不同道集域下的时差表现形式不同，因此可以在井间地震处理中利用中值滤波技术实现直达波和多次波的衰减，反射波的增强以及上、下行反射波的分离。同时，由于中值滤波处理对地震原始资料畸变程度较小，因此可以提高井间地震资料反射成像的质量。3.2多道中值滤波在分离VSP波场中的应用？

多道中值滤波处理流程图

4.中值滤波特点

中值滤波是一个非线性过程，最大优点是算法简单且去噪效果明显。中值滤波具有 如下特点:(1)中值滤波绝对阻止噪声峰值。

(2)中值滤波不改变阶跃函数在空间、时间上的位置。（3）消除尖峰波以及增强部分有效波；

（4）对野外原始地震资料信息的畸变和负面影响较小等优点。(5)中值滤波平滑三角波,其平滑作用随着中值滤波长度N的增加而增加,当其达到一定长度时,可将三角波平滑为具有相等幅度的理想的直流分量。

(4)中值滤波平滑矩形波,若中值滤波足够长时,矩形波被完全平滑。

(6)中值滤波由于沒有像计算均值那样的统计效应,对随机出现的小的振幅值有时不能完全平滑。中值滤波实际上是一个平滑滤波,经过其处理之后,主频向低频移动,高频成分受到损害,正是由于其平滑作用,使处理后的地震数据波形过于一致而显呆板,有一些信息不可能客观反映出来。

5.总结

中值滤波器是一种特殊的非线性滤波器，与线性滤波器不同。线性滤波器的振幅、相位谱完全决定了滤波器在频率域和时间域的特征。与之相比，中值滤波器对脉冲的响应是零，在频率域没有“真正”的振幅谱和相位谱。所以，尽管中值滤波器原理很简单，但了解其特性比了解线性滤波器特性要困难得多。

虽然中值滤波技术在信号处理领域得到很大的重视，特别在非平稳信号的处理中取得了较大的成功，然而中值滤波的一个严重不足是引起相对滤波窗口而言较为“细小”的信号细节结构的破坏和丢失，在图象处理中，中值滤波的这一缺憾要比在一维信号的处理中更加显著。原因主要来自两个方面：

第一，二维信号几乎没有根信号，也就是说几乎所有的二维信号经中值滤波以后都要受到不同程度的破坏；

第二，图象中的某些诸如细线，拐角等细节结构往往包含重要的信息，这些结构的破坏或丢失往往比噪声本身更为不可接受。所以，保护细节的中值类滤波器的研究成为非线性滤波器研究的一个重要方面。多级中值滤波则正是人们在努力寻求的兼有细节保护和噪声抑制的优良特性的滤波器结构。研究表明，长度较小的窗口能够较好地保护信号的细节信息，但却不能有效地滤除随机噪声；而长度较大的滑动窗口能更好地抑制噪声，同时却严重地损失重要信息。根据噪声性质自动改变滤波窗口长度的滑动加权中值滤波器，更好地适应去噪的需求。

第四篇：Matlab频域滤波练习

%% fftshift 对数变换,所应用的图片本身很简单，就只有黑白2种颜色 clc clear

f = imread('.imagesdipum_images_ch04Fig0403(a)(image).tif');imshow(f)

title('原始图像')

imfinfo('.imagesdipum_images_ch04Fig0403(a)(image).tif');%此处如果用Imfinfo(f)就会报错fft

%没有居中的傅里叶频谱

F=fft2(f);%进行二维快速傅里叶变换，其结果和DFT的一样，只是计算机的计算速度变快了而已，因而叫fft S=abs(F);%求傅里叶变换后的幅值 figure,subplot(121),imshow(S,[])

title('傅里叶频谱图像1');%title函数一定要放在坐标显示的下一句才有效。

subplot(122),imshow(S),title('傅里叶频谱图像2')

%居中的傅里叶频谱

Fc=fftshift(F);%将频谱图像原点移至图像矩形中间 S1=abs(Fc);

figure,subplot(121),imshow(S1,[]);%加了第二个参数后显示的图像正常

subplot(122),imshow(S1);%当没有第二个参数时，显示的图像为竖线加一些孤立的黑点，why？

%使用对数后视觉增强后的傅里叶频谱 S2=log(1+S1);imshow(S2,[]);

%%用fftshift对数变换显示稍复杂的图片

clc clear

f=imread('.imagesdipum_images_ch04Fig0409(a)(bld).tif');imshow(f);

f=double(f);%其实这句不要试过对后面的变换结果也没有影响 F=fft2(f);

Fc=fftshift(F);S=abs(Fc);

S2=log(1+S);%如果没有这句的话，那么根本看不到细节的图，所以一定要用对数压缩增加对比度

figure,imshow(S2,[])

%%试一下ifft功能

clc clear

f=imread('.imagesdipum_images_ch04Fig0409(a)(bld).tif');imshow(f);

f=double(f);f(1:8,1:8)F=fft2(f);

f=ifft2(F);%取傅里叶变换后的反傅里叶变换,直接是整数，并不需要像某些代码一样取real部分 f(1:8,1:8)

%%先0扩充再滤波 clc clear f=imread('.imagesdipum_images_ch04Fig0405(a)(square_original).tif');imshow(f);

[M N]=size(f);

F=fft2(f);%没有经过0扩充直接计算fft sig=10;%高斯滤波参数

H=lpfilter('gaussian',M,N,sig);

G=H.*F;%加了.号的乘法表示对应每个元素都相乘，没加.号的乘法说明是真正的矩阵乘法

g=real(ifft2(G));

figure,imshow(g,[]);

PQ=paddedsize(size(f));%PQ为0扩展的面积，这里设置为与图像的大小相同？这里size(f)的大小为256*256

Fp=fft2(f,PQ(1),PQ(2));%如果fft2有3个参数的话，则是进行了0扩展了

Hp=lpfilter('gaussian',PQ(1),PQ(2),2*sig);%构造高斯滤波器 Gp=Hp.*Fp;%对0扩展后的图像进行高斯滤波 gp=real(ifft2(Gp));figure,imshow(gp,[])'

gpc=gp(1:size(f,1),1:size(f,2));%取gp图像的1到f的第一维的行，1到f第二维的列部分，即取与图像大小相同的部分

figure,imshow(gpc,[]);%此时显示的结果应该与g图像一样。

%直接使用空间域滤波

h=fspecial('gaussian',15,7);gs=imfilter(f,h);

figure,imshow(gs,[])

%% clc clear

f=imread('.imagesdipum_images_ch04Fig0405(a)(square_original).tif');f=double(f);imshow(f);

PQ=paddedsize(size(f));%size(f)的大小为256*256

sig=10;

H = lpfilter('gaussian',PQ(1),PQ(2),2*sig);% PQ=[512 512] figure,mesh(abs(H(1:10:512,1:10:512)));%画出滤波器的三维空间图 g=dftfilt(f,H);

figure,imshow(g,[])

%%空间滤波和频域滤波的比较

clc clear

f=imread('.imagesdipum_images_ch04Fig0409(a)(bld).tif');imshow(f);

F=fft2(f);

S=fftshift(log(1+abs(F)));

S=gscale(S);%用gscale来进行将数据缩放到默认值0~255 figure,imshow(S);

h=fspecial('sobel');%构造sobel空间滤波器 figure,freqz2(h);%查看滤波器的图形

PQ=paddedsize(size(f));

H=freqz2(h,PQ(1),PQ(2));%将空间滤波器h转换成频域滤波器H，但是滤波器的原点在矩阵的中心处

H1=ifftshift(H);%原点迁移到左上角

figure,mesh(abs(H1(1:20:1200,1:20:1200)));figure,subplot(121),imshow(abs(H),[]);subplot(122),imshow(abs(H1),[]);

gs=imfilter(double(f),h);%空间滤波，其实就是边缘检测 gf=dftfilt(f,H);%频域滤波

subplot(121),imshow(gs,[]),subplot(122),imshow(gf,[]);%将2种滤波结果显示出来

figure,imshow(abs(gs)>0.2*abs(max(gs(:))));%只显示最大值的20%的边缘

figure,imshow(abs(gf)>0.2*abs(max(gf(:))));

d=abs(gs-gf);a=max(d(:))b=min(d(:))

%%构造低通滤波器

clc clear

f=imread('.imagesdipum_images_ch04Fig0413(a)(original_test_pattern).tif');

imshow(f)

F1=fft2(f);

imshow(log(1+abs(fftshift(F1))),[]);%直接显示取对数后的傅里叶变换

PQ=paddedsize(size(f));

[U V]=dftuv(PQ(1),PQ(2));%这个函数还真不明白什么意思，help中解释的是计算网格频率矩阵U和V，这2个矩阵是用来频率滤波的 D0=0.05*PQ(2);

F=fft2(f,PQ(1),PQ(2));%用0扩充大小为PQ(1)*PQ(2)，然后fft变换 figure,imshow(log(1+abs(fftshift(F))),[]);%显示0扩充后的fft变换图

H=exp(-(U.^2+V.^2)/(2*(D0^2)));%构造频域滤波算子 figure,imshow(fftshift(H),[]);%显示平移后滤波器算子 g=dftfilt(f,H)

figure,imshow(g,[]);

%%绘制一个低通滤波器的mesh图 clc clear

H1=lpfilter('gaussian',500,500,50)%构造一个中心点在左上角的高斯低通滤波器,size为500*500

H2=fftshift(H1);%将中心点平移至矩阵中心

mesh(H1(1:10:500,1:10:500)),figure,mesh(H2(1:10:500,1:10:500));%绘出2者的mesh图

axis([0 50 0 50 0 1]);%xy轴范围0~50，z轴范围0~1

figure,subplot(121),imshow(H1,[]),subplot(122),imshow(H2,[]);%绘出2者的灰度图

%%绘制一个低通滤波器的mesh图 clc clear

H=fftshift(hpfilter('ideal',500,500,100));%构造理想的高通滤波器 mesh(H(1:10:500,1:10:500));axis([1 50 1 50 0 1]);colormap([0 0 0]);%mesh网格全部用黑色画 axis off%关掉坐标轴显示 grid off%关掉网格显示 figure,imshow(H,[])

%%使用高通滤波器

clc clear

f=imread('.imagesdipum_images_ch04Fig0413(a)(original_test_pattern).tif');

imshow(f)

title('原始图像')

PQ=paddedsize(size(f));D0=0.05*PQ(1);

H=hpfilter('gaussian',PQ(1),PQ(2),D0);g=dftfilt(f,H);

figure,imshow(g,[])

%%将高通滤波和直方图均衡结合起来

clc clear

f=imread('.imagesdipum_images_ch04Fig0419(a)(chestXray_original).tif');imshow(f)

title('原始图像')

PQ = paddedsize(size(f));D0 = 0.05*PQ(1);

HBW=hpfilter('btw',PQ(1),PQ(2),D0,2);%构造高通Butterworth滤波器，截断频率为D0 H=0.5+2*HBW;

gdw=dftfilt(f,HBW);

gbw=gscale(gdw);%将灰度值自动缩放到0~255之间 figure,subplot(121),imshow(gdw,[]);title('高通滤波后图像');

gbe=histeq(gbw,256);%直方图均衡化 subplot(122),imshow(gbe,[]);title('高通滤波且直方图均衡化后图像');

ghf=dftfilt(f,H);%H是HBW的线性变换 ghf=gscale(ghf);

figure,subplot(121),imshow(ghf,[]);title('高频强调滤波后图像');

ghe=histeq(ghf,256);

subplot(122),imshow(ghe,[]);title('高频强调且均衡化后图像');

%%fftshift和ifftshift的加深理解 clc clear

A=[2 0 0 1

0 0 0 0

0 0 0 00 0 4]

B=fftshift(A)%中心点平移

C=fftshift(fftshift(A))%还原成A

D=ifftshift(fftshift(A))%也同样还原成A E=ifftshift(A)%平移结果和B一样

第五篇：基于matlab的语音信号滤波处理——数字信号处理课程设计

数字信号处理课程设计

题目： 学院： 专业： 班级： 学号： 姓名： 指导教师：

基于matlab的语音信号滤波处理

物理与电子信息学院电子信息工程

摘要：

语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴学科，是目前发展最为迅速的学科之一，通过语音传递信息是人类最重要，最有效，最常用和最方便的交换信息手段，所以对其的研究更显得尤为重要。

Matlab语言是一种数据分析和处理功能十分强大的计算机应用软件，它可以将声音文件变换成离散的数据文件，然后用起强大的矩阵运算能力处理数据。这为我们的本次设计提供了强大并良好的环境！

本设计要求自己录制一段自己的语音后，在MATLAB软件中采集语音信号、回放语音信号并画出语音信号的时域波形和频谱图。再在Matlab中分别设计不同形式的FIR数字滤波器。之后对采集的语音信号经过不同的滤波器（低通、高通、带通）后，观察不同的波形，并进行时域和频谱的分析。对比处理前后的时域图和频谱图，分析各种滤波器对于语音信号的影响。最后分别收听进行滤波后的语音信号效果，做到了解在怎么样的情况下该用怎么样的滤波器。

目录

1.设计内容……………………………………………………………4 2.设计原理……………………………………………………………4 2.1语音信号的时域分析…………………………………………4 2.2语音信号的频域分析…………………………………………5 3.设计过程……………………………………………………………5 3.1实验程序源代码………………………………………………6

3.1.1原语音信号时域、频域图………………………………6

3.1.2低通滤波器的设计………………………………………6

3.1.3高通滤波器的设计………………………………………7

3.1.4带通滤波器的设计………………………………………8

3.1.5语音信号的回放………………………………………9 3.2调试结果描述…………………………………………………10 3.3所遇问题及结果分析…………………………………………15

3.3.1所遇主要问题…………………………………………16

3.3.2结果分析………………………………………………16 4.体会与收获…………………………………………………………17 5.参考文献……………………………………………………………17

1.设计内容：

1．首先录制好一段自己的语音。

2．用Matlab分别设计好3种类型的滤波器（指标自己确定）：低通型、高通型、带通型。3．用Matlab将语音信号进行采样，并分别将其通过所设计的3种滤波器。4．用Matlab自带的语音返回函数收听滤波后的语音信号，分析并比较其与原语音信号的差异。

2.设计原理：

语音信号时一种非平稳的时变信号，它带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。语音信号分析的目的就在于方便有效的提取并表示语音信号所携带的信息。语音信号处理可以分为时域和变换域等处理方法，其中时域分析是最简单的方法，直接对语音信号的时域波形进行分析，崎岖的特征参数主要有语音的短时能量，短时平均过零率，短时自相关函数等。2.1语音信号的时域分析

信号提取：通过图形用户界面上的菜单功能按键采集电脑上的一段音频信号，完成音频信号的频率，幅度等信息的提取，并得到该语音信号的波形图。

信号调整：在设计的用户图形界面下对输入的音频信号进行各种变化，如变化幅度、改变频率等操作，以实现对语音信号的调整。

2.2语音信号的频域分析

信号的傅里叶表示在信号的分析和处理中起着重要的作用。因为对于线性系统来说，可以很方便地确定其对正弦或复指数和的响应，所以傅里叶分析方法能完善地解决许多信号分析和处理问题。另外，傅里叶表示使信号的某些特性变得更明显，因此，它能更深入地说明信号的各项红物理现象。

由于语音信号时随着时间变化的，通常认为，语音是一个受准周期脉冲或随机噪声源激励的线性系统的输出。输出频谱是声道系统频率响应与激励源频谱的乘积。身份到系统的频率响应及激励源都是随时间变化的，因此一般标准的傅里叶表示虽然适用于周期及平稳随机信号的表示，但不能直接用于语音信号。由于语音信号可以认为在短时间内，近似不变，因而可以采用短时分析法。

1.信号变换：在用户图形界面西啊对采集的语音信号进行Fourier等变换，并画出变换前后的频谱图和倒谱图。

2.信号滤波：滤除语音信号中的噪音部分，可以采用抵用滤波、高通滤波、带通滤波，并比较各种滤波后的效果。

3．设计过程：

3.1实验程序源代码（原语音信号存放在e:下）：

3.1.1．原语音信号的时域、频域图

[x1,fs,bits]=wavread('e:txwz.wav');%sound(x1,fs,bits);figure(1);

plot(x1);%做原始语音信号的时域图形 title('原始语音信号');xlabel('时间 t');ylabel('音量 n');figure(2);y1=fft(x1);%做length(x1)点的FFT y1=fftshift(y1);%平移，是频率中心为0 derta_fs = fs/length(x1);%设置频谱的间隔，分辨率 plot([-fs/2:derta_fs: fs/2-derta_fs],abs(y1));%画出原始语音信号的频谱图

title('原始语音信号的频谱');grid on;3.1.2低通滤波器的设计

%低通滤波：截止频率4000，阻带衰减20dB，过渡带宽0.1π

fc1=4000;N1=2*pi*0.9/(0.1*pi)wc1=2*pi*fc1/fs;if rem(N1,2)==0 N1=N1+1;end Window= boxcar(N1+1);%长度为N1的矩形窗Window b1=fir1(N1,wc1/pi,Window);

figure(3);freqz(b1,1,512);title('低通滤波器的频率响应');x1_low = filter(b1,1,x1);%对信号进行低通滤波 figure(4);plot(x1_low);title('信号经过低通滤波器(时域)');figure(5);plot([-fs/2:derta_fs:fs/2-derta_fs],abs(fftshift(fft(x1_low))));title('信号经过低通滤波器（频域）');3.1.3高通滤波器的设计

%高通滤波：截止频率4000，阻带衰减40dB,过渡带宽0.1π

fc2=4000;N2=2*pi*3.1/(0.1*pi)wc2=2*pi*fc1/fs;N2=N2+mod(N2,2);Window=hanning(N2+1);b2=fir1(N2,wc2/pi,'high',Window);figure(6);freqz(b2,1,512);%数字滤波器频率响应 title('高通滤波器的频率响应');

x1_high = filter(b2,1,x1);%对信号进行高通滤波 figure(7);plot(x1_high);title('信号经过高通滤波器(时域)');figure(8);plot([-fs/2:derta_fs:fs/2-derta_fs],abs(fftshift(fft(x1_high))));title('信号经过高通滤波器（频域）')3.1.4带通滤波器的设计

%带通滤波：下截止频率4000，上截止频率8000，阻带衰减20dB,过渡带宽度0.1π

f1=4000;f2=8000;%带通滤波器的通带范围 w1=2*pi*f1/fs;w3=w1+0.1*pi;w2=2*pi*f2/fs;w4=w2-0.1*pi;w=[(w1+w3)/2,(w2+w4)/2];B=0.1*pi;N3=ceil(2*0.9*pi/B);N3=N3+mod(N3,2);Window=boxcar(N3+1);b3=fir1(N3,w/pi,'stop',Window);%带通滤波器

figure(9);freqz(b3,1,512);%数字滤波器频率响应 title('带通滤波器的频率响应');x1_daitong = filter(b3,1,x1);%对信号进行带通滤波 figure(10);plot(x1_daitong);title('信号经过带通滤波器(时域)');figure(11);plot([-fs/2:derta_fs:fs/2-derta_fs],abs(fftshift(fft(x1_daitong))));title('信号经过带通滤波器（频域）');3.1.5语音信号的回放（分别执行）

sound(x1,fs,bits);%原始信号

sound(x1_low,fs,bits);%经过低通滤 sound(x1_high,fs,bits);%经过高通滤波 sound(x1_daitong,fs,bits);%经过带通滤波

3.2调试结果描述

原始语音信号的时域图形：

原始语音信号频谱：

低通滤波器的频率响应：

信号经过低通滤波后的时域波形：

信号经过低通滤波后的频域波形

高通滤波器的频率响应：

信号经过高通滤波后的时域波形：

信号经过高通滤波后的频域波形：

带通滤波器的频率响应：

信号经过带通滤波后的时域波形：

信号经过带通滤波后的频域波形：

3.3所遇问题及结果分析

3.3.1所遇主要问题

1．在高通与带通滤波器的设计时老是报错，但同样的用法在低通滤波器中就可以实现

b2=fir1(N2,wc2/pi,'high',Window);??? Error using ==> fir1 The window length must be the same as the filter length.其要求在fir函数中所选用的窗长要和滤波器长度一致。但在参考书上指出，滤波器阶数必须为窗长加1。经上网查询后，原来高通、带阻滤波器的阶数应该控制为奇数，因为如果阶数为偶数，则在π点必有一零点，这对于高通带阻来说是不允许的，故取阶数为奇数，而你FIR1滤波器阶数为M+1阶，所以你的M必须为偶数，所以可以将程序中去窗长算法由原程序的: N2=2*pi*0.9/(0.1*pi);if rem(N2,2)==0 N2=N2+1;End 和: N3=2*pi*0.9/(0.1*pi);if rem(N3,2)==0 N3=N3+1;End 改为了：

N2=N2+mod(N2,2);和: B=0.1*pi;（B为过渡带宽）N3=ceil(2*0.9*pi/B);N3=N3+mod(N3,2);2.在设计高通滤波器时先是使用的矩形窗，用矩形窗验证出来的结果中低频语音分量依旧很强，不能将其全部抑制在0，之后换窗，选着了最小衰减可以到达53dB的海明窗，再次试验，非常成功！3.3.2结果分析

经过回放三个不同类型滤波器输出的语音信号，并与原语音信号对比得到了如下结论。

语音高频成分音质非常尖锐，齿音中，声音有些暗淡。语音低频成分音质沉稳，空间感觉强，语音浑厚。语音中频成分音质有力度，有通透感。

4．体会与收获

以往都是通过课本来感性的认知语音信号，通过本次的课程设计，让我对语音信号有了一个较为实际的认识。于此同时，让我再次把数字信号处理及数字滤波器的设计方法重新进行了复习和学习。而对于Matlab，也再次让我感受到了其功能的强大。最为重要的是，本次课 17

程设计让我重新审视了学习的过程：只去做实验是不行的，首先还是要思考，遇到了问题查书籍，百度搜索也只是一种手段，更加重要的是想，再理解，只有这样才能真正的做好实验。

5.参考文献

《数字信号处理》钱同惠、百度文库