图像处理汉明窗

2020-04-19 11:49:55 ASKLW 阅读数 52

图像解析力算法—SFR(Spatial Frequency Response)原理分析(一)中,我们已经分析了SFR的前四个步骤,接下来,我们继续讨论以下的五个步骤

4、重新定位ROI,获得ESF

5、对获得的ESF进行四倍超采样

6、通过差分运算获得LSF

7、对LSF应用汉明窗

8、进行DFT运算

 

4、重新定位ROI,获得ESF

这一步其实比较复杂,我也不确定在我的讲述之下,大家是否能够听懂,我尝试用简单的方式讲解一下。

首先,转换坐标轴,将坐标轴转换到计算出来的矩心直线上。如图所示:

5、对获得的ESF进行四倍超采样

然后,我们将每一行中X轴坐标相等的像素值累加起来,然后求均值后得到下面第一行的数组。

 

这里是将整张图片的像素转换为一个长度为rows的数组,然后进行4倍的扩增。然后其他没有像素的地方,是进行一个向前寻找非0的像素值进行替换,这样就获得了ESF。

6、通过差分运算获得LSF

这不用多说,就是得到的ESF数组元素进行差分,获得LSF。一下是函数公式:

7、对LSF应用汉明窗

对上面的LSF数组进行汉明窗处理,这一步主要也是看公式,如下:

8、进行DFT运算

上图是具体的DFT公式表达,下图是代码中的计算过程

最后的到的SFR数组就是空间频域响应的值。

好了,SFR的具体过程就是这么多,介绍到这里。因为这也是博主本人自己对SFR的理解的归纳总结,如果有大佬发现有错漏的地方,希望不吝赐教,指出纠正。有不明白的地方也可以留言,我们一起探讨。

今天就先写到这里吧。接下来,博主就会对SFR中的主要步骤的函数进行分析。

2019-04-18 11:46:56 mary_0830 阅读数 487

(七)约束的最小二乘法(规则化)滤波

说明:另一个容易接受的线性复原方法叫约束的最小二乘法滤波,也称为规则化滤波。两个函数f和h的2D离散卷积数学表达式如下所示:

在这里插入图片描述
在matlab中是用函数deconvreg实现的,基本语法为:

frest = deconvreg(g,PSF,NOISEPOWER,RANGE);  %g是被污染的图像,frest是复原的图像

如果没有NOISEPOWER和RANGE这两个参数,则会产生逆滤波效果。

实验代码及结果如下所示:

f = imread('C:\Users\Public\Pictures\Sample Pictures\horse.jpg');
f1 = im2double(rgb2gray(f));
g = imnoise(f,'poisson');
PSF = fspecial('motion',4,45);
frest1 = deconvreg(g,PSF,4);  %NOISEPOWER的初试估计为4
frest2 = deconvreg(g,PSF,0.4,[1e-7 1e7]);
subplot(121),imshow(frest1,[]),title('参数为4的滤波');
subplot(122),imshow(frest2,[]),title('参数为0.4、范围为[1e-7 1e7]的滤波');

(图一)
在这里插入图片描述

PSF = fspecial('motion',4,45);

(图二)
在这里插入图片描述

PSF = fspecial('motion',2,45);

(图三)
在这里插入图片描述
(图四)
在这里插入图片描述对比上面四幅图,没有看出同一幅图的两幅图有什么差别。通过修改deconvreg的参数噪声能力和范围之后,得到如下结果:

frest1 = deconvreg(g,PSF,0.8);  %NOISEPOWER的初试估计为0.8
frest2 = deconvreg(g,PSF,0.8,[1e-4 1e4]);

(图五)
在这里插入图片描述

结论:
通过修改参数可以使得约束最小二乘法的滤波效果变好;更换图像后,滤波效果也是很明显,图像并没有出现模糊,但还是出现了噪声。调整参数后,虽然滤波结果仍存在很多不足,但至少基本可以看出细节及内容。修改参数后得到的图五,有范围的约束最小二乘法滤波的图像显示,图像没有那么模糊,马匹后面的部门没有虚化和晕染,比前面没有范围的图像清晰一些,但是还是存在噪声很严重。应该是RANGE限制了图像的噪声范围,使得噪声的随机性降低,固定于一个范围内,因此显现的图像会比无范围的更加好一些。

注意: 若复原图像呈现运算中由离散傅里叶变换引入的振铃,在通常情况下,在调用函数deconvreg之前,需要使用函数edgetaper来帮助。

(八)利用露西-查理德森算法的迭代非线性复原

说明:这种算法是一种非线性算法,简称为L-R算法。L-R算法起源于最大似然公司,在这个方程式中,图像是用泊松统计建模的。数学表达式如下:

在这里插入图片描述

L-R算法用deconvlucy函数实现,基本语法是:

f = deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);   %f是复原图像,g是退化图像,PSF是点扩散函数,NUMIT是迭代次数(默认为10) 

其中,DAMPAR是标量,指定了结果图像与原始图像g之间的偏离阈值。当像素值偏离原值的范围在DAMPAR内时,就不用再迭代。这既抑制了图像像素噪声,又保留图像细节。默认值为0(无衰减)。WEIGHT是数组,大小与g相同,它为每个像素都施以权重以反映像素的质量。

注意: 若复原图像呈现运算中由离散傅里叶变换引入的振铃,在通常情况下,在调用函数deconvlucy之前,需要使用函数edgetaper来帮助。

实验代码:

g = checkerboard(8);                  %产生一幅64×64的方形图像
subplot(131),imshow(pixeldup(g,8)),title('原始图像');  %用函数pixeldup将图像大小放大为512×512像素
PSF = fspecial('gaussian',7,10);      %产生7×7大小、标准差为10的高斯PSF
SD = 0.01;
g = imnoise(imfilter(g,PSF),'gaussian',0,SD^2);   %添加高斯噪声
subplot(132),imshow(g,[]),title('高斯噪声污染模糊后的图像');
DAMPAR = 10*SD; %对参数DAMPAR赋10倍于SD的值
LIM = ceil(size(PSF,1)/2);        %设置参数WEIGHT
WEIGHT = zeros(size(g));          %大小为64×64,值为0的4像素宽的边界
WEIGHT(LIM + 1:end - LIM,LIM + 1:end - LIM) = 1;    %其他像素为1
NUMIT = 5;
f5 =deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);
subplot(133),imshow(pixeldup(f5,8),[]),title('用参数为5的迭代复原图像');

实验结果:

在这里插入图片描述
实验代码:

NUMIT = 10;
f10 =deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);
subplot(131),imshow(pixeldup(f10,8),[]),title('用参数为10的迭代复原图像');

NUMIT = 25;
f25 =deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);
subplot(132),imshow(pixeldup(f25,8),[]),title('用参数为25的迭代复原图像');

NUMIT = 200;
f200 =deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);
subplot(133),imshow(pixeldup(f200,8),[]),title('用参数为200的迭代复原图像');

实验结果:
在这里插入图片描述

结论: 即使增加了迭代次数,但在复原结果方面并没有明显的改善。只能对比出迭代次数越多,得到的图像比之前亮一些。非线性迭代法滤波并没有很好的滤波效果,还是使得滤波图像更加棋盘化,这样的效果是由于像素灰度级不够而引起的。也说明了露西-查理德森算法的迭代非线性复原没有很好的还原高斯噪声污染过的图像。

(九)盲去卷积

说明:盲去卷积是对随机噪声污染的量进行估计时采用的最佳策略。
在matlab中用函数deconvblind来执行盲去卷积,基本语法如下:

[fr,PSF] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);

使用最大似然算法对图像g解卷积,返回去模糊图像fr和恢复的点扩散函数PSF。 生成的PSF是与INITPSF相同大小的正数组,归一化,所以它的总和增加到1。PSF的恢复受其初始猜测大小INITPSF的影响较大,而其值较小(一个数组是一个更安全的猜测)。为了改善恢复,可以传入附加参数,如NUMIT,DAMPAR,WEIGHT。

NUMIT(可选)是迭代次数(默认值为10)。DAMPAR(可选)是一个数组,用于指定图像g(根据泊松噪声的标准偏差)的结果图像的阈值偏差,低于此值会发生阻尼。 对于在DAMPAR值内偏离其原始值的像素,迭代被抑制。 这可以抑制这些像素中的噪音,并在其他地方保留必要的图像细节。 默认值为0(无阻尼)。WEIGHT(可选)分配给每个像素以反映相机的拍摄质量。 将一个坏像素分配给零权值,从而排除该像素。 您可以根据平场校正的数量来调整自己的体重,而不是给予好像素的权重。 默认值是与输入图像I大小相同的单位数组。

注意: 若复原图像呈现运算中由离散傅里叶变换引入的振铃,在通常情况下,在调用函数deconvblind之前,需要使用函数edgetaper来帮助。

实验代码:

g = checkerboard(8);                  %产生一幅64×64的方形图像
g = imnoise(imfilter(g,PSF),'gaussian',0,SD^2);   %添加高斯噪声
PSF = fspecial('gaussian',7,10);      %产生7×7大小、标准差为10的高斯PSF
subplot(221),imshow(pixeldup(PSF,73),[]),title('原始PSF');
DAMPAR = 10*SD; %对参数DAMPAR赋10倍于SD的值
LIM = ceil(size(PSF,1)/2);        %设置参数WEIGHT
WEIGHT = zeros(size(g));          %大小为64×64,值为0的4像素宽的边界
WEIGHT(LIM + 1:end - LIM,LIM + 1:end - LIM) = 1;    %其他像素为1
SD = 0.01;
g = imnoise(imfilter(g,PSF),'gaussian',0,SD^2);   %添加高斯噪声
INITPSF = ones(size(PSF));
NUMIT = 5;
[g5,PSF5] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);
subplot(222),imshow(pixeldup(PSF5,73),[]),title('用参数5迭代时PSF的估计值');
NUMIT = 20;
[g20,PSF20] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);
subplot(223),imshow(pixeldup(PSF20,73),[]),title('用参数20迭代时PSF的估计值');
NUMIT = 10;
[g10,PSF10] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);
subplot(224),imshow(pixeldup(PSF10,73),[]),title('用参数10迭代时PSF的估计值');

实验结果:
在这里插入图片描述实验结论:
以上图只是改变了参数迭代的次数,其他的并没有修改。从修改后的看出,迭代次数为5的估计值图像的四个角颜色变深了,图像整体的颜色也变深了;迭代次数为10的虽然整体颜色也有点深,但是没有迭代次数是5那么严重,随着迭代次数的增加,正如参数为20的图像,在这几幅图中接近原图像。

(十)来自投影的图像重建

说明:一维投影的图像重建是医学中的重要应用之一。

在这里插入图片描述
当收集到一副图像各个角度的投影后,并希望通过这些投影的图像重建原图像。

(1)函数radon

函数radon是用来为给定的二维矩阵列产生一组平行射线投影,基本语法是:

R =radon(I,theta);  %I是二维阵列,theta是角度值的一维阵列

更一般的语法是:

[R,xp] = radon(I,theta); %xp包含沿着x’轴的坐标值

用来产生众所周知的图像的某个有用函数语法为:

P = phantom(def,n);  %def是指定产生头部模型类型的字符串,n是行数和列数(默认值为256)

其中,字符串def的可用值为:‘Shepp-Logan’(对比度很低)和’ModifiedShepp-Logan’(改进了对比度)。

实验代码:

f = imread('C:\Users\Public\Pictures\Sample Pictures\horse.jpg ');
g1=rgb2gray(f); 
g1(100:500,250:350) = 1;
g2 = phantom('Modified Shepp-Logan',600);
subplot(221),imshow(g1),title('图一');
subplot(222),imshow(g2),title('图二');
theta = 0:0.5:179.5;
[R1,xp1] = radon(g1,theta);        %执行雷登变换
[R2,xp2] = radon(g2,theta);
R1 = flipud(R1');           %投影图像
R2 = flipud(R2');
subplot(223),imshow(R1,[],'XData',xp1([1 end]),'YData',[179.5 0])      %显示投影图像
axis xy
axis on
xlabel('\rho'),ylabel('\theta')
subplot(224),imshow(R2,[],'XData',xp2([1 end]),'YData',[179.5 0])
axis xy
axis on
xlabel('\rho'),ylabel('\theta')

实验结果:
在这里插入图片描述
在这里插入图片描述

(2)函数iradon

函数iradon从给定的取自不同角度的一组投影来重建一幅图像(切片)。简单说来,就是iradon计算反randon变换。基本语法是:

I = iradon(R,theta,interp,filter,frequency_scaling,output_size)

其中,
R是反投影数据,列是从左到右以角度渐增的函数来组织的一维反投影;
theta描述角度,据此来获取头像;
interp是字符串,定义了用于产生最终重建图像的内插方法;
filter指定了在滤波反投影计算中使用的滤波器;
frequency_scaling是处于(0,1)范围内的标量,通过改变频率轴的比例来修改滤波器,默认值为1。若frequency_scaling小于1,滤波器被压缩,以适应【0,frequency_scaling】范围;
output_size是标量,规定了重建图像的行列数,若没有指定,则由投影的长度来确定:output_size = 2*floor(size(R,1)/(2*sqrt(2)))

实验代码:

theta = 0:0.5:179.5;
R1 = radon(g1,theta);  
R2 = radon(g2,theta);
f1 = iradon(R1,theta,'none');        
f2 = iradon(R2,theta,'none');
subplot(121),imshow(f1,[]),title('不进行滤波的图像一');
subplot(122),imshow(f2,[]),title('不进行滤波的图像二');

实验结果:
在这里插入图片描述
在这里插入图片描述

结论: 由图可知,不进行滤波的效果图最差,而且模糊了图像,只能看到一个大致的轮廓,基本看不到原来的细节(就像上图有马匹的图一样,只能大概看出有一匹马,其他细节完全无法识别)。

实验代码:

f1_ram = iradon(R1,theta);          %默认值滤波
f2_ram = iradon(R2,theta);
subplot(121),imshow(f1_ram ,[]),title('默认滤波的图像一');
subplot(122),imshow(f2_ram ,[]),title('默认滤波的图像二');

实验结果:
在这里插入图片描述
在这里插入图片描述
实验代码:

f1_hamm = iradon(R1,theta,'Hamming');     %汉明窗滤波
f2_hamm = iradon(R2,theta,'Hamming');
subplot(121),imshow(f1_hamm ,[]),title('用汉明窗滤波的图像一');
subplot(122),imshow(f2_hamm ,[]),title('用汉明窗滤波的图像二');

实验结果:
在这里插入图片描述
在这里插入图片描述

实验代码:

f1_hann= iradon(R1,theta,'Hann');     %韩窗滤波
f2_hann = iradon(R2,theta,'Hann');
subplot(121),imshow(f1_hann,[]),title('用韩窗滤波的图像一');
subplot(122),imshow(f2_hann ,[]),title('用韩窗滤波的图像二');

实验结果:
在这里插入图片描述

结论: 默认值滤波器滤波后的效果比较差,汉明窗和韩窗的效果都基本还原了图像,但也是会模糊图像。图像二的三种滤波法看出,默认值滤波器滤波的效果差,不仅模糊图像而且降低对比度;图像一(图像是马匹)对比看出,默认值滤波器模糊图像。

实验代码:

f1_near = iradon(R1,theta,'nearest');    %最近邻内插法
f1_lin = iradon(R1,theta,'linear');      %线性内插法
f1_cub = iradon(R1,theta,'cubic');       %三次内插法
subplot(131),imshow(f1_near,[]),,title('用最近邻内插法滤波的图像一');
subplot(132),imshow(f1_lin,[]),,title('用线性内插滤波的图像一');
subplot(133),imshow(f1_cub,[]),,title('用三次内插滤波的图像一');

实验结果:
在这里插入图片描述

在这里插入图片描述

结论: 第一次用图显示不出三种内插法哪种比较好,之后用了马匹的图像,明显的看出线性内插和三次内插法比最近邻内插法滤波效果好。放大之后发现,最近邻内插法使得图像不仅模糊还产生噪声;线性内插法增加了对比度;三次内插发则细化了细节部分,对比更加明显。

(3)扇形射束的数据处理

给定扇形射束数据,工具箱使用的方法是把扇形射束变换为与之对应的平行射束,然后用早期的平行射束得到反投影。在扇形射束中,射线可以描述成一条直线L(ρ,θ)。平行射束和扇形射束存在某种对应关系:
在这里插入图片描述
其中,在这里插入图片描述是相应的平行射束投影。
matlab中用函数fanbeam产生扇形射束投影,基本语法为:

B = fanbeam(g,D,param),val1,param2,val2,...)   %g是包含被投射物体的图像,D是从扇形射束的顶点到旋转中心的距离

假定旋转中心是图像的中心,规定D大于g的一半:

D = K*sqrt(size(g,1)^2 + size(g,2)^2)/2  %K是大于1的常数

函数fanbeam通过计算对于任何角度覆盖全部图像需要多少射线数来决定传感器的数目。
在这里插入图片描述

实验代码及效果如下:

g1 = zeros(600,600);
g1(100:500,250:350) = 1;
g2 = phantom('Modified Shepp-Logan',600);
D = 1.5*hypot(size(g1,1),size(g2,2))/2;
B1_line = fanbeam(g1,D,'FanSensorGeometry','line','FanSensorSpacing',1,'FanRotationIncrement',0.5);
B1_line = flipud(B1_line');
B2_line = fanbeam(g2,D,'FanSensorGeometry','line','FanSensorSpacing',1,'FanRotationIncrement',0.5);
B2_line = flipud(B2_line');
B1_arc = fanbeam(g1,D,'FanSensorGeometry','arc','FanSensorSpacing',0.08,'FanRotationIncrement',0.5);
B1_arc = flipud(B1_arc');
B2_arc = fanbeam(g2,D,'FanSensorGeometry','arc','FanSensorSpacing',0.08,'FanRotationIncrement',0.5);
B2_arc = flipud(B2_arc');
subplot(221),imshow(B1_line,[]),title('g1直线扇形射束图像');
subplot(222),imshow(B2_line,[]),title('g2直线扇形射束图像');
subplot(223),imshow(B1_arc,[]),title('g1圆弧扇形射束图像');
subplot(224),imshow(B2_arc,[]),title('g2圆弧扇形射束图像');

在这里插入图片描述
我们发现相应的外形十分不同,扇形射束投影出现了“歪斜”,这是扇形直接对于平行射束的结果。直线扇形射束图像比圆弧扇形射束更暗,这是与传感器的形状有关,圆弧扇形射束扫到检测的范围比较宽且呈现弧度,因此显示出的图像比较亮一些。

函数ifanbeam可用来给定一组扇形射束投影得到滤波反投影图像:

I = ifanbeam(B,D,...,param1,val1,param2,val2,..);

实验代码及效果如下:

g = phantom('Modified Shepp-Logan',600);
D = 1.5*hypot(size(g1,1),size(g2,2))/2;
B1 = fanbeam(g,D);
f1= ifanbeam(B1,D);
subplot(131),imshow(f1,[]),title('默认值产生并重建幻影图像');
B2 = fanbeam(g,D,'FanRotationIncrement',0.5,'FanSensorSpacing',0.5);
f2 =  ifanbeam(B2,D,'FanRotationIncrement',0.5,'FanSensorSpacing',0.5,'Filter','Hamming');
subplot(132),imshow(f2,[]),title('旋转和传感器间隔增量为0.5、汉明窗产生并重建幻影图像');
B3 = fanbeam(g,D,'FanRotationIncrement',0.5,'FanSensorSpacing',0.5);
f3 =  ifanbeam(B3,D,'FanRotationIncrement',0.5,'FanSensorSpacing',0.05,'Filter','Hamming');
subplot(133),imshow(f3,[]),title('传感器间隔增量为0.05、汉明窗产生并重建幻影图像');

在这里插入图片描述

扇形和平行投影间转换的函数:

P = fan2para(F,D,param1,val1,param2,val2,...)  %扇形数据转换为平行射束数据
F = para2fan(P,D,param1,val1,param2,val2,...)  %平行数据转换为扇形射束数据

实验代码及效果图:

g1 = phantom('Modified Shepp-Logan',600);
g1(100:500,250:350) = 1;
g2 = phantom('Modified Shepp-Logan',600);
D = 1.5*hypot(size(g1,1),size(g2,2))/2;
B1_line = fanbeam(g1,D,'FanSensorGeometry','line','FanSensorSpacing',1,'FanRotationIncrement',0.5);
B2_arc= fanbeam(g2,D,'FanSensorGeometry','arc','FanSensorSpacing',0.08,'FanRotationIncrement',0.5);
P1_line = fan2para(B1_line,D,'FanRotationIncrement',0.5,'FanSensorGeometry','line','FanSensorSpacing',1,'ParallelCoverage','halfcycle','ParallelRotationIncrement',0.5,'ParallelSensorSpacing',1);
P2_arc = fan2para(B2_line,D,'FanRotationIncrement',0.5,'FanSensorGeometry','arc','FanSensorSpacing',0.08,'ParallelCoverage','halfcycle','ParallelRotationIncrement',0.5,'ParallelSensorSpacing',1);
P1_line = flipud(P1_line');
P2_arc = flipud(P2_arc');
subplot(121),imshow(P1_line,[]),title('扇形射束投影');
subplot(122),imshow(P2_arc,[]),title('平行射束投影');

在这里插入图片描述

2017-06-21 21:08:58 gxiaoyaya 阅读数 11738

为什么要加汉明窗?什么叫加窗?

 在信号处理中,可以说加窗处理是一个必经的过程,因为我们的计算机只能处理有限长度的信号,因此原始信号X(t)要以T(采样时间)截断,即有限化,成为XT(t)后再进一步处理,这个过程序就是加窗处理,但什么时候用什么窗呢?这时我们就要对所需用到的函数窗做一定的了解。在平时,我们用得最多的是矩形窗,这个也很容易理解,好像我们屋子里的窗口一样,透过窗口我们可以看到外面的世界,但在如果我们理窗口远一些的话,我们的看到的范围将减少,越远就越小。实际的信号处理过程中,我们用的矩形窗,但矩形窗在边缘处将信号突然截断,窗外时域信息全部消失,导致在频域增加了频率分量的现象,即频谱泄漏。避免泄漏的最佳方法是满足整周期采样条件,但实际中是不可能做到的。对于非整周期采样的情况,必须考虑如何减少加窗时造成的泄漏误差,主要的措施是使用合理的加窗函数,使信号截断的锐角钝化,从而使频谱的扩散减到最少。

    首先介绍一下为什么要用函数窗:函数窗的主要用于对截断处的不连续变化进行平滑,减少泄漏。此外,加窗处理还有很多其它的原因,如减少噪声干扰、限定测试的持续时间、从频率接近的信号中分离出幅值不同的信号……

常见的几种窗的基本指标:

 

一个窗是否合适:窗谱主瓣宽度就尽可能的窄,且能量集中在主瓣内,以获得较陡的过渡带;窗谱旁瓣与主瓣相比应尽可能的小,旁瓣能量衰减要快,以利于增加阻带衰耗。

汉明窗就是信号窗口的一种,在matlab中执行命令,画出plot(hamming(100))的图如下:


 

它主要部分的形状像sin(x)在0到pi区间的形状,而其余部分都是0.这样的函数乘上其他任何一个函数f,f只有一部分有非零值。

 

为什么汉明窗这样取呢?因为之后我们会对汉明窗中的数据进行FFT,它假设一个窗内的信号是代表一个周期的信号。(也就是说窗的左端和右端应该大致能连在一起)而通常一小段音频数据没有明显的周期性,加上汉明窗后,数据形状就有点周期的感觉了。

 

因为加上汉明窗,只有中间的数据体现出来了,两边的数据信息丢失了,所以等会移窗的时候,只会移1/3或1/2窗,这样被前一帧或二帧丢失的数据又重新得到了体现。

 

简单的说汉明窗就是个函数,它的形状像窗,所以类似的函数都叫做窗函数。希望你能明白。

 

2.加Hanmming窗的作用

现在在看G.723.1,对语音编码刚入门,

发现在对信号进行LPC分析前,对信号乘以一个Hamming 窗,

乘法是:信号直接乘以一个HammingWindowTable中的值,这个加窗有什么作用?

如果是限制带宽的话, 在时域应对信号应做卷积的, 不明白,请赐教

 

因为要处理的是无限长序列中的一段,所以必须对这段序列加窗采集出来。

 

典型的窗口大小是25ms,帧移是10ms。汉明窗函数为

            W(n,α ) = (1 -α ) - α cos(2*PI*n/(N-1)),0≦n≦N-1

    一般情况下,α取0.46 。

谁能解释一下这个函数吗?我实在是不理解,谢谢.

 

由于直接对信号(加矩形窗)截断会产生频率泄露,为了改善频率泄露的情况,加非矩形窗,一般都是加汉明窗,因为汉明窗的幅频特性是旁瓣衰减较大,主瓣峰值与第一个旁瓣峰值衰减可达40db。

 

举例:

a=wavread('jiasiqi.wav');   %将音频信号jiasiqi.wav读入
subplot(2,1,1),                  %分配画布,一幅图上共两个图,这是第一个
plot(a);title('original signal');  %画出原始信号,即前面这个音频信号的原始波形
grid                                    %添加网格线
N=256;                               %设置短时傅里叶变换的长度,同时也是汉明窗的长度
h=hamming(N);                   %设置汉明窗
for m=1:N                       %用汉明窗截取信号,长度为N,主要是为了减少截断引起的栅栏效应等
b(m)=a(m)*h(m)
end
y=20*log(abs(fft(b)))           %做傅里叶变换,取其模值,即幅频特性,然后用分贝(dB)表示
subplot(2,1,2)                     %分配画布,第二副图
plot(y);title('短时谱');            %画出短时谱
grid                                        %添加网格线



2014-02-28 13:12:40 ningyaliuhebei 阅读数 44505

图像处理算法工程师(索贝公司)

一、填空:
1、常用的插值方法有:最近邻插值、双线性插值、立方卷积插值。
2、常用的边缘检测算子有:一阶: Roberts Cross算子, Prewitt算子, Sobel算子, Canny算子, 罗盘算子
二阶: Marr-Hildreth。
3、能够表征一副图像的基本特征有:灰度值、纹理、形状
4、FIR滤波器设计中常用的窗函数:三角形(Bartlett)窗、汉宁(Hanning)窗、汉明(Hamming)窗、

布莱克曼(Blackman)窗
5、视频流处理单元是:音频流处理单元是:
6、(2006)10转换成16进制:7d6
7、X86体系中,常用寄存器中经常用来存储数据的是:
8、C++类中三种存取权限类型:private、public和protected。
9、视频帧播放速度的单位是:PAL制式是——25fps,NTSC是——30fps。
10、mfc中,CFile类最大支持读写——字节,Windows下动态加载一个动态函数名————


汉王机器视觉(软件工程师):

1.以下变量pValue分别是什么类型?并请谈谈你对static和const的理解。

http://bbs.chinaunix.net/thread-143183-1-1.html

贴两个链接供参考:http://www.cnblogs.com/dc10101/archive/2007/08/22/865556.html

http://blog.csdn.net/ccccdddxxx/article/details/7085165

(1)static int(*pValue)[10];静态的指向整形数组的指针(数组指针)

(2)int(*pValue[10])(int);

原题主要两部分,第一部分关于算法基础;第二部分关于图像模式知识,题目如下:

①指针的概念理解【概念理解。指针数组、数组指针、指向函数的指针、指针函数...】

②sizeof计算 【比较简单,网上到处都是...】

③用C编程实现字符串匹配 【这个有点难度,编的有点离谱...】

④用C编程实现一维最大熵阈值分割  【需要知道熵的表达式才好做,刚好笔试前我用过此算法,还算有印象...】

⑤给出图像像素表绘制灰度直方图、图像大小变换后重新绘制图像像素表、阈值化后绘制图像像素表;【图像的一些基本概念】

⑥车辆检测中,常见的去除人物干扰算法;【写了几种自己稍微了解的..】

⑦利用贝叶斯的简单运算分类;【模式识别中最基本的知识,却真记不住,还好,概率论学的还不错,后面回来查了一下,竟然也对了,原来贝叶斯只是一个简单的概率计算而已,以前上课怎这么难呢?...】

⑧一副彩色图像,里面含苹果、菠萝、梨子、香蕉,如何进行特征提取和分类器设计 【说了我熟悉的用哪个Gabor滤波器提取特征,利用纹理特征提取分类...后面work才知道,原来是想考彩色分割相关只是哦...嘿嘿】

 面 试

说明:主管比较厉害,问的好似没什么逻辑,却几乎包含了所有常见的面试题,最糟糕的,完后,你还觉得是在聊天,能记住的真的很少,脑海中还有影响的几个暂时记录如下.

内容:主要是两方面,一是对于简历上的实践经历/项目详细询问(特别注重细节,问的很细很细);另外问道的问题大致如下(零散):

①近2-3年的发展规划?

②如何处理校园职务/活动与科研/学习的关系(两种似乎不同的性格)?

③工资待遇要求?

④如何处理公司实习与学校学习/毕业任务的关系,孰轻孰重?

⑤说一件你印象最深的事情?

⑥为何从原来公司(简历上说明的)辞职?

⑦如何理解责任和道德的?法律约束的行为(应该做和必须做的行为),社会公认的行为

一、逻辑题

1、住在某个旅馆的同一房间的四个人A、B、C、D正在听一组流行音乐,她们当中有一个人在修指甲,一个人在写信,一个人躺在床上,另一个人在看书。
(1)A不在修指甲,也不在看书;
(2)B不躺在床上,也不在修指甲;
(3)如果A不躺在床上,那么D不在修指甲;
(4)C既不在看书,也不在修指甲;
(5)D不在看书,也不躺在床上。
她们各自在做什么呢?

2、如果我们在21的2与1之间添加进去若干个0,使它变成:20…01,现在问:这种20…01的数中,是否有能被21整除的?如果没有,那是为什么?如果有,那么有多少个?

3、一个农夫发现围成正方形的围栏比长方形的节省4个木桩但是面积一样.羊的数目和正方形围栏的桩子的个数一样但是小于36,问有多少羊?说出你的计算过程。

二、编程题

1、如果已经定义了:float x=1.5; inta=1,b=3,c=2;
以下两个switch语句,哪个正确,哪个错误,为什么?

switch(x)
{

case 1.0:printf("*\n");
case 2.0:printf("**\n");

}
switch(a+b)
{

case 1:printf("*\n");
case 2+1:printf("**\n");

}

 

2、设int arr[] = {6, 7, 8, 9, 10};

    int *ptr = arr;

    *(prt++)+=123;

printf("%d,%d",*ptr,*(++ptr));

请问输出结果是什么?

 

3、请问以下程序的输出结果是什么?

int first()

{

int i=1;

return(i++);

}

 

int second()

{

static i=1;

return(i++);

}

 

void main()

{

  inti;

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

{

printf(“first   %d\n”,first());

}

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

{

printf(“second   %d\n”, second());

}

}

 

4、已知strcpy函数的原型是:char* strcpy(char*strDest,const char*strSrc);
(1)不调用库函数,实现strcpy函数。
(2)解释为什么要返回char*。

三、信号处理

1、有一个正弦信号隐藏在高斯白噪声中,请问如何检测出该正弦信号的频率?

 

2、FIR滤波器和IIR滤波器有什么区别,各有什么优缺点?

 

3、A wheel, rotating at 6Hz, is seen in a dark room by means of astrobe light flashing at a rate of 8Hz. Determine the apparent rotational speedand sense of rotation of the wheel. Repeat the question if the flashes occur at12Hz, 16Hz or 24Hz.

 

4、A filter is described by the following sample processing algorithm:

     For each input x(n) do:

            

a) Determine the transfer function H(z)=Y(z)/X(z) of the filter.

   b)Draw the canonical realization form of the system.

四、请回答下面的问题(必答题)

1、在您以前进行的有关图像信号处理的工作中,

[1] 请列举您单独负责的项目(或者作为其中主要负责人),并简述您负责的部分。 

[2] 项目中您认为困难点在哪里?是如何解决的?

 

2、您觉得自己应聘这个职位的优势和不足是什么?您最喜欢从事的工作是什么?如果您加入北阳,您将为未来的工作做出哪些准备?

 

五、选做题

1、汉明窗(Hamming)和矩形窗是信号处理中常用的窗函数,请阐述二者的适用场合。

 

2、Define the difference between an emulator and a simulator. What arethe benefits offered by each and which is the most suitable when developingreat-time DSP software?



2019-01-02 00:07:14 liubing8609 阅读数 976

窗函数在图像处理中的应用

1. 频谱混乱的三角函数图像

下图是一个45度倾斜的单一频率的余弦函数图像,请注意图中的边界都不是均匀过渡到外界的,全是不连续的跳变。

https://img-blog.csdn.net/20180524085021172?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

下面来看看这幅图的频谱会是什么样?

https://img-blog.csdn.net/20180524085232658?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

频谱混乱的原因有两个:

第一个重要的原因就是一个域(不论是时域还是频域)的不连续导致了另一个域的振荡,拖尾,泄漏。所以下图中的频谱就会看起来非常混乱。

第二个原因就是DFT的周期性。当使用MATLAB中的函数FFTFFT2的时候,一定要意识到,所做的计算本质上离散傅里叶级数(为什么是级数而不是变换,我会再写一篇文章说明)Discrete Fourier transformDFT)。要知道,DFT的计算隐含着默认的周期性,这一周期性不论是在时域还是在频域都是无穷的,都是无限循环的,从二维信号上讲,即沿着X方向循环复制,也在Y方向循环复制。下图中用红色方框框出的只是其中的一个周期。

https://img-blog.csdn.net/20180425172531208?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

请注意上图中左图中(空间域)余弦信号图和他相邻的图像在红色边缘处的各种不连续(discontinuity)。这种不连续直接导致了右图中(频率域)及其混乱,复杂的频谱。所以,在处理任何图像时,千万不要忽视了图像的边缘。比如说,FFT之前的0填充,这也是一种会经常用到的避免频谱泄漏的有效方式。

汉宁窗

现在对原始图像进行加窗,这里选择的是汉宁窗(hanning window)。下图为汉宁窗的时间域函数图像及其频谱特征。

总的来说,窗函数都比较光滑,就好像一个高斯滑动均值滤波器磨平了原始信号中边缘的不连续。

https://img-blog.csdn.net/20180426094105120?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

https://img-blog.csdn.net/20180426102704718?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

请注意上图中我用红框框出的几个重要参数,这是衡量窗函数的三个重要指标。

对图像加窗

加窗后的频谱图远比加窗前的更集中,更清晰。

https://img-blog.csdn.net/20180426101333318?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

https://img-blog.csdn.net/20180426102528252?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

2. 著名的Cameraman.tiff

刚才的那个例子中,空间域的不连续是非常明显,显而易见的(TIPS:一个域的不连续导致了另外一个域的频振荡和拖尾,不只是从空间域/时间域到频域是这样)。下面用图像处理中的经典图像来看另外一种情况。在这个例子中,边界上的跳变远没有上一个例子那么明显,但也会引起频谱的泄露,污染了(smear)真实的频谱信息,而这种情况最为常见,最容易被忽视。

下图为cameraman的原图和频谱图,注意频谱图中我用箭头标出的一条白色的中轴线,也是频谱的泄漏造成的。

https://img-blog.csdn.net/20180426103441563?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

产生这根白线的原因主要是,原始图像在DFT周期复制的时候,头顶的蓝天部分较亮,而脚下草坪部分较暗。由于图像横向边缘的灰度值突变(见下图),在频谱图中的Y方向产生了频谱的泄露。当然了,图像的纵向边缘也能看到一些灰度的不平滑过渡,势必也会引起图像频谱在X方向的改变,但在此例中并不明显。

https://img-blog.csdn.net/20180426105521393?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

下面看看加窗后的图像和加窗后频谱的改变。

https://img-blog.csdn.net/2018042611025985?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

https://img-blog.csdn.net/20180426110546479?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

频谱图中间的那条白线,没了!

3. 上三角矩阵和下三角矩阵

理论上讲上三角矩阵和下三角矩阵的频谱都应该只有一条斜线,不管是左斜还有右斜,但是往往大家看到的上/下三角矩阵的频谱图当中还有一个十字架。这都是因为没有加窗出现了DFT周期性的边界不连续。

https://img-blog.csdn.net/20180427213534403?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

注意图中的十字架,这是本不该出现在频谱图中的成分。

https://img-blog.csdn.net/20180523173419890?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

下图是加窗后的真实频谱,见下图。

https://img-blog.csdn.net/20180524084635844?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RhZHV6aW1hbWE=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70

总结,这篇文章用三个简单的例子,粗略的解释了窗函数在DIP中的应用。实际图像处理的工作中,但凡是用到了频谱操作,很多时候正确的使用窗函数往往是大有裨益的!

【转载】https://blog.csdn.net/daduzimama/article/details/80079215

窗函数

阅读数 7457

海明窗滤波

阅读数 2056