精华内容
下载资源
问答
  • 频率域法把图像看成一种二维信号进行二维傅里叶变换的信号增强,采用低通滤波法可以去掉图像的噪声;采用高通滤波法,则可增强边缘等高频信号,使模糊的图像变得清晰。滤波是一种重要的图像预处理手段,目前,...

    摘要:图像的频域滤波是图像增强的一种方法。图像增强是图像处理的方法之一,有频率域法和空间域法。频率域法把图像看成一种二维信号,对其进行二维傅里叶变换的信号增强,采用低通滤波法可以去掉图像的噪声;采用高通滤波法,则可增强边缘等高频信号,使模糊的图像变得清晰。滤波是一种重要的图像预处理手段,目前,图像的频域滤波应用非常广泛,例如,基于频域滤波的织物疵点的检测、颤振信号处理中的时频域滤波方法等。

    本文将主要研究常用图像滤波方法分析、频域滤波的特点及使用场合以及其如何用MATLAB实现图像的增强处理。

    关键词:图像增强、频域滤波、MATLAB

    目录

    摘要

    Abstract

    1绪论-1

    1.1数字图像处理的研究现状及发展方向-1

    1.2研究内容和研究方法-2

    1.2.1研究内容-2

    1.2.2研究方法-2

    2数字图像处理-2

    2.1数字图像的基本概念-2

    2.2数字图像处理软件-4

    2.2.1 MATLAB简介-4

    2.2.2其他图像处理软件的介绍-4

    3图像的频率域变换-5

    3.1连续傅里叶变换-5

    3.1.1一维连续傅里叶变换-5

    3.1.2二维连续傅里叶变换-6

    3.2离散傅里叶变换-6

    3.3快速傅里叶变换(FFT)-7

    3.4本章小结-7

    4图像增强-8

    4.1图像噪声-8

    4.2频率域图像增强概念-8

    4.3 空间域滤波和频率域滤波的比较-9

    5频率域滤波增强-9

    5.1一些基本的滤波器介绍-10

    5.2频率域低通滤波器-10

    5.2.1理想低通滤波器-10

    5.2.2巴特沃斯(Butterworth)低通滤波器-10

    5.2.3高斯低通滤波器-11

    5.3频率域高通滤波器-11

    5.3.1理想高通滤波器-11

    5.3.2巴特沃斯型高通滤波器-11

    5.3.3高斯高通滤波器-12

    5.4同态滤波器的简单介绍-12

    5.5本章小结-12

    6频域滤波增强应用实例-13

    6.1理想型低通滤波器的应用-13

    6.2巴特沃斯(Butterworth)型低通滤波器的应用-13

    6.3高斯型低通滤波器的应用-15

    6.4理想型高通滤波器的应用-16

    6.5巴特沃斯型高通滤波器的应用-16

    6.6高斯型高通滤波器的应用-17

    6.7频域滤波应用于舰船热尾流红外图像-18

    7结论-20

    致谢-20

    参考文献-21

    附录-21

    展开全文
  • matlab简单分析频域滤波和时域滤波

    千次阅读 2020-04-06 19:42:33
    对信号进行低通滤波的一种快速方法,低频信号幅度设为1,负频率镜像过去,把虚部设为0 Fs = 2048; dt = 1.0/Fs; T = 1; N = T/dt; t = linspace(0,T,N); x1 = sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*...

    华中科技大学《数字信号分析理论实践》第六单元 信号数字滤波的概念 学习总结记录

    频域滤波

    在这里插入图片描述

    • 对信号进行低通滤波的一种快速方法,低频信号幅度设为1,负频率镜像过去,把虚部设为0
      在这里插入图片描述
    Fs = 2048;
    dt = 1.0/Fs;
    T = 1;
    N = T/dt;
    t = linspace(0,T,N);
    x1 = sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
    subplot(411);
    plot(t,x1);
    axis([0,0.1,-2,2]);
    P = fft(x1,N);
    Pyy = 2*sqrt(P.*conj(P))/N;
    f = linspace(0,Fs/2,N/2);
    subplot(412)
    plot(f,Pyy(1:N/2));
    P1(N) = 0;
    for k = 1:N
    	P1(k) = real(P(k))+1i*imag(P(k));
    end
    for k = 200:N-200
    	P1(k) = 0;
    end
    Pyy = 2*sqrt(P1.*conj(P1))/N;
    subplot(413)
    plot(f,Pyy(1:N/2));
    x2 = ifft(P1);
    subplot(414);
    plot(t,real(x2));
    axis([0,0.1,-2,2]);
    

    在这里插入图片描述

    • 长数据滤波
      在这里插入图片描述

    数字时域滤波器——关键得到Z变换形式的滤波器

    • 数字差分——简单的高通滤波器
      x(n)x(n+1)x(n1)2ΔtH(z)=az1az1,a=12Δtx'(n)\approx\frac{x(n+1)-x(n-1)}{2\Delta t}\Rightarrow H(z)=az^1-az^{-1},a = \frac{1}{2\Delta t}
    a = 1/(2*dt);
    For K = 1 To N-1
    	x(k) = a*x(k+1)-a*x(k-1)
    Next
    

    在这里插入图片描述

    Fs = 5000;
    dt = 1/Fs;
    T = 1;
    N = T/dt;
    t = linspace(0,T,N);
    a = 5;
    f = 2;
    y = a*sin(2*pi*f*t)+0.3*sin(2*pi*50*t);
    subplot(211)
    plot(t,y);
    x(N) = 0;
    for i = 2:N-1
    	x(i) = (y(i+1)-y(i-1))/2;
    end
    subplot(212)
    plot(t,x);
    
    • 数据平滑——简单的低通滤波器
      在这里插入图片描述
    Fs = 5000;
    dt = 1/Fs;
    T = 1;
    N = T/dt;
    t = linspace(0,T,N);
    a = 5;
    f = 2;
    y = a*sin(2*pi*f*t)+0.8*sin(2*pi*500*t);
    subplot(211)
    plot(t,y);
    x = y;
    for i = 1:N-5
    	x(i) = (y(i)+y(i+1)+y(i+2)+y(i+3)+y(i+4))/5;
    end
    subplot(212)
    plot(t,x);
    

    在这里插入图片描述

    • 高斯平滑滤波器——将系数改进,中间值权重最大
      在这里插入图片描述
    • 数字滤波器设计——求取滤波器系数
      • 时域数字滤波就是将信号和滤波器的单位脉冲响应 (Z域滤波器系数) 进行卷积分,寻找所需的滤波器的单位脉冲响应的过程称为滤波器设计,信号卷积分过程称为滤波
        在这里插入图片描述
    • 注意对滤波器起始点进行处理
      在这里插入图片描述
    展开全文
  • 1、频率滤波图像的空间域滤波:用各种模板...图像频域滤波,先把图像转换到频域空间,然后不同的频率点进行滤波,使用信号处理的技术,图像实现滤波。比如实现图像的轮廓提取,在空间域滤波中我们使用一个拉普拉...

    1、频率滤波

    图像的空间域滤波:用各种模板直接与图像进行卷积运算,实现对图像的处理,这种方法直接对图像空间操作,操作简单。图像处理不仅可以在空间域进行还可以在频率域进行,把空间域的图像开窗卷积形式,变换得到频率域的矩阵点乘形式得到比较好的效果。图像频域滤波,先把图像转换到频域空间,然后对不同的频率点进行滤波,使用信号处理的技术,对图像实现滤波。比如实现图像的轮廓提取,在空间域滤波中我们使用一个拉普拉斯模板就可以提取,而在频域内,我们使用一个高通滤波模板,可以实现轮廓的提取。

    图像特征与像素点数值的关系:

    图像尺寸418*564,如果把每一行所有像素,一行564个点的灰度作为一维向量画图,取前三行画成三条曲线,就得到了下面的图形。

    subplot(2,1,1);imshow(img), title('原始图像');line1 = img(1, :);line2 = img(2, :);line3 = img(3, :);subplot(2,1,2);hold onplot(line1, 'r');plot(line2, 'g');plot(line3, 'b');

    输出结果如下图。

    40a38144f5e63708b619c0b17d7faa10.png

    图像前几行的特征如下图(为了方便观察提取了不止三行):

    e7038d8d5e049e221f3d598333eb4d97.png

    从图像矩阵前几行的像素特征可见,图像平坦的地方像素曲线也平坦,图像亮的地方就是图像像素剧烈变化的地方,图像像素值发生较大差异的地方也是图像发生突变的地方,这些位置一般就是图像轮廓。图像的频率体现了图像中灰度变化剧烈程度,是灰度在平面空间上的梯度。所以前面讲的空域滤波,可以使用平滑来滤除噪声实现平滑,从上面曲线图像上看,可以按照信号处理思想来理解,平滑滤波就是频率低通。对应的,频率高通滤波就是空域的提取边界。

    2、图像傅里叶变换

    傅里叶变换可以将一个时域信号转换成在不同频率下对应的振幅及相位,其频谱就是时域信号在频域下的表现,而反傅里叶变换可以将频谱再转换回时域的信号。

    imgPath = 'E:opencv_picsrc_picpic4.bmp';img = imread(imgPath);img=rgb2gray(img); f=fft2(im2double(img)); %FFTF=f; %FFT频谱T=log(F+1); %频谱对数变换subplot(1,2,1),imshow(img),title('原始图像');subplot(1,2,2),imshow(T,[]),title('原始图像其频谱图');

    输出结果如下,幅度图在四个角出现亮光。因为实信号以fs为采样速率的信号在 fs/2处混叠,所以实信号fft的结果中前半部分对应[0, fs/2],后半部分对应[ -fs/2, 0]。横向和纵向都需要把频率转换到[ -fs/2, fs/2]的区间,可以使用fftshift函数。

    cade4e26d36177cd8100bc03118ddff9.png

    在数字图像处理中,常常需要将F(u,v)的原点移到N*N频域的中心,以便能清楚地分析傅里叶谱的情况,平移前空域、频域原点均在左上方。而fftshift的作用就是这样,将0频谱移到正中心。

    Y = fftshift(X) 通过将零频分量移动到数组中心,重新排列傅里叶变换 X。

    如果 X 是向量,则 fftshift 会将 X 的左右两半部分进行交换。

    如果 X 是矩阵,则 fftshift 会将 X 的第一象限与第三象限交换,将第二象限与第四象限交换。

    如果 X 是多维数组,则 fftshift 会沿每个维度交换 X 的半空间。

    e617afffa846458fba921c3bd8112f77.png
    a =     1     2     3     4     6     7     8     9>> fftshift(a)ans =     8     9     6     7     3     4     1     2

    平移后的图:

    de5ab14a27417810c911981c2aca3969.png

    傅里叶变换的能量集中在频率很小的圆内,当D0增大时能量衰减很快,高频部分虽然携带的能量很少,但是包含丰富的边界和细节信息,所以当截止频率D0变小时,虽然亮度足够(因能量损失不大),但图像变模糊。

    5f5e11bd85452c6538c58792c5565050.png

    3、频域滤波

    傅里叶变换可以把图像从空域变换到频域,而傅里叶反变换可以将图像的频谱变换为空域图像。可以利用图像空域和频域之间的对应关系,尝试将空域卷积滤波变换为频域滤波,而后再将频域滤波处理后的图像反变换回空间域,从而达到图像增强的目的。

    G=imnoise(img,'gaussian', 0, 0.05);%模拟均值为0方差为0.05的高斯噪声,H=fft2(im2double(G)); %FFTH=fftshift(H); %FFT频谱平移T=log(abs(H)); %频谱对数变换subplot(2,2,1), imshow(G),title('添加高斯噪声图像');subplot(2,2,2),imshow(T, []),title('频谱图');
    b774569bf6debe6c8348e0beb0157831.png

    使用傅里叶变换得到图像的频域幅度,应用傅里叶反变换得到空域图像.

    img2 = (ifft2(ifftshift(H))); %===频域的图反变换到空域img3 = im2uint8(mat2gray(img2)); %===取其灰度图subplot(2,2,3),imshow(img3);title('anti-Fourier');

    反变换效果如下图。

    ee05c55bcedcbe56150e0874e0c1ec39.png

    那么就可以在频率对图像进行滤波,然后反变换回空域就得到滤波后的图像。频率图可见在图像中心,是频率为0位置,所以如果限定只有低频通过则为低通,如果抑制频谱图中心位置,则是高通。可见图像在频域滤波是非常直观的。

    1)低通滤波器:

    cb3ca6d450a673b4aa17c0b508e85910.png

    其中D(u,v)为频率域上(u,v)点到中心的距离,D0由自己设置.

    4eee85d0ab72bc96c369c1df537d2dc3.png

    2)高通滤波器:

    高通滤波器同低通滤波器非常类似,只不过二者通过的波正好是相反的.

    8e0ad37168a2fdb806c9ad9e6ab03318.png
    7b6657b4275c15a99e41b38459382e7f.png

    3)高斯低通滤波

    fc3f58417e051f37e4b449e3b45c4140.png
    e564144f1554ee3a0604ac902053a8f5.png

    4)高斯高通滤波

    7a8bf77099e02433f4b351031c742536.png
    8b1098e7a5a7a4789dd88af407b98c96.png

    高斯频率低通滤波:

    S=fftshift(fft2(G)); %G是添加高斯噪声的图像[M,N]=size(S);                                          d0=5; %GLPF滤波,d0=5,15,30                    n1=floor(M/2);                          n2=floor(N/2);                           for i=1:M     for j=1:N        d=sqrt((i-n1)^2+(j-n2)^2);                 h=1*exp(-1/2*(d^2/d0^2));          S(i,j)=h*S(i,j);                       endend S=ifftshift(S);                           S=uint8(ifft2(S));    subplot(2,2,4),imshow(S),title('高斯低通滤波图像');

    输出图像如下:

    35753fc5150ae459f70779f41be27d33.png

    高斯滤波函数标准差变大,则高斯函数的傅里叶变换图像频率范围变宽,允许通过的频率范围变大,对应空域允许通过的图像变换更加多,可以允许的噪声更多。相反,如果标准差变小,频域范围变窄,允许通过的频率范围变小,对应空域对低频要求更加严格被抑制的高频更多,所以图像更加平滑,可能出现模糊现象。下面展示了不同的σ标准差情况下的滤波效果。

    05dbf223e4c405da71dde83b83ac6ebf.png

    改为高通滤波,结果如下图:

    7b4095e99814369bbc42329bfb6d1cd7.png
    展开全文
  • 使用Matlab对给定的信号和预期需求,设计了几种连续时间滤波器,包括四种最常见的连续滤波器:巴特沃斯(Butterworth)、切比雪夫滤波器I型和II型(chebyshev)、考尔滤波器(Cauer);验证了滤波器性能,并进行了分析和...

    0754c77e2c8bc208e5424be32280b1c0.png
    使用Matlab对给定的信号和预期需求,设计了几种连续时间滤波器,包括四种最常见的连续滤波器:巴特沃斯(Butterworth)、切比雪夫滤波器I型和II型(chebyshev)、考尔滤波器(Cauer);验证了滤波器性能,并进行了分析和对比,给出了代码,并将滤波器设计过程常用操作简单封装成函数,简化编程。

    待处理信号

    e97b0c3a42c5c4b1ff2f9425e5f21c8d.png
    不同音所在的频率

    这里使用gamme.m函数,生成一段依次从Do、 Re、Mi、...一直到升Do的音频,函数输入是每个频率保持的长度 和 采样频率,输出是一段二维数组,即生成的音频。

    function [gam,t] = gamme(duree,fe)
    
    % [do re mi fa sol la si doo]
    freqnotes = [262 294 330 349 392 440 494 523];
    t = 0:1/fe:duree;
    
    gam = [];
    for note=1:8, 
        gam = [gam, sin(2*pi*freqnotes(note)*t)];
    end
    N = length(gam);
    t = (0:N-1)/fe;

    调用函数,生成音频,并使用matlab的soudnsc函数试听

    clear all
    fe = 8192 ; %采样频率
    duree = 1; %每个音符持续1秒
    [sig,t] = gamme(duree,fe);
    soundsc(sig, fe); % listen

    对原始信号进行频率分析

    fmin = 250;
    fmax = 550;
    
    f = [0:length(sig)-1]*(fe/length(sig));
    S = fft(sig)/fe;
    
    figure()
    plot(f, abs(S));
    axis([fmin fmax 0 0.5])

    b616d5bd90742653686566ff8b7de83f.png
    八个音,每个的幅值都是一样的,是因为它们在信号中的比例(时长)相同

    连续时间滤波器

    滤波器指标

    我们希望除去后三个音,即La So 升Do,对照这个表,得到滤波器指标

    e97b0c3a42c5c4b1ff2f9425e5f21c8d.png

    类型:低通

    通带容许的最小增益:-3dB

    阻带增益:-40dB

    过渡带中心频率:420Hz (在Sol和La之间)

    过渡带宽度:100Hz (这里100Hz有些宽,为的是让不同的滤波器的对比更加明显,我们会在后面看到它的影响)

    e63413902e908322b2262de44a8cbb4c.png
    实际的低通滤波器模型,从左到右:通带、过度带、阻带
    • 之所以通带并不是完全通过,是因为大多滤波器会有抖动存在,并且即使没有抖动,滤波器的频域响应也不会是水平直线,因此设计时应当留有一定的裕度

    9077adbffb68b1f25bbbc10acaac8d8c.png
    一个带有抖动的低通滤波器的频域响应

    b543b8eb5111849a4e474ddae5ff8efa.png
    三个没有抖动的低通滤波器

    滤波器设计

    1、常量定义

    Ap = 3;
    Aa = 40; %这里的增益不用加负号
    Deltaf = 100;
    fc = 420;
    fa = fc + Deltaf/2;
    fp = fc - Deltaf/2;

    2、计算滤波器阶数

    [nc1 wnc1]= buttord(2*pi*fp,2*pi*fa,Ap,Aa,'s')
    [nc2 wnc2]= cheb1ord(2*pi*fp,2*pi*fa,Ap,Aa,'s')
    [nc3 wnc3]= cheb2ord(2*pi*fp,2*pi*fa,Ap,Aa,'s')
    [nc4 wnc4]= ellipord(2*pi*fp,2*pi*fa,Ap,Aa,'s') %cauer滤波器

    - 结果分析

    f035548214cd4bc8546eedb1733b01c5.png
    • nc代表需求的阶数,wnc代表对应滤波器的【由通带到过度带】的频率,这个频率对于不同的滤波器来说是不一样的,直观上可以理解为,这是因为每个滤波器的频域响应的”形状“都不一样,因此它们与【-3dB的水平线】的交点的横坐标不一样。

    d1fc1de83f2cadb78a96404ee3fd9fe4.png
    不同滤波器与【-3dB的水平线】的交点的横坐标不一样,因此计算出来的wnc也不一样
    • 命名的‘c’代表连续时间,即模拟滤波器
    • 此外,可以看到,Butterworth的阶数最高,是20,两个Chebyshev阶数相同,是8,而Cauer滤波器只需要5阶就可以
      • 通常情况下,滤波器的阶数对应了滤波器的复杂程度和成本,阶数越高,成本越高,为了达到相同的性能指标,Butterworth需要更高的阶数,Cauer需要最少的,但是有失就有得,Butterworth的好处在于它在通带和阻带都没有抖动,是光滑的,Chebyshev I型在通带有抖动,Chebyshev II型在阻带有抖动,Cauer在通带和阻带都有抖动。
      • 这体现了一种tradeoff的思想,如果想要滤波器在通带和阻带没有抖动,那么代价就是需要更高的滤波器阶数,如果想要更简单的滤波器结构,即阶数低,那么代价是滤波器的频域响应会有抖动
      • 这个抖动的影响,我们会在后面的仿真中看到

    3、计算低通滤波器参数

    [numc1, dend1] = butter(nc1,wnc1,'low','s');
    [numc2, dend2] = cheby1(nc2,Ap,wnc2,'low','s');
    [numc3, dend3] = cheby2(nc3,Aa,wnc3,'low','s');
    [numc4, dend4] = ellip(nc4,Ap,Aa,wnc4,'s');
    • 值得注意的是,四个函数的参数是不同的,butter默认只需要输入【阶数,ord函数计算出来的频率】,但cheby1需要输入【阶数,通带增益,ord函数计算出来的频率】,cheby2则是把通带增益换成了阻带增益,ellip,即cauer滤波器则需要两个都输入。这与滤波器的频域响应抖动是对应的:Cheby1 通带有抖动,Cheby2阻带有抖动,Cauer 通带阻带都有抖动。

    对比四种滤波器的频域响应

    f = linspace(fmin*0.1,fmax*2,1000);
    Hc1 = freqs(numc1,dend1,2*pi*f);
    Hc2 = freqs(numc2,dend2,2*pi*f);
    Hc3 = freqs(numc3,dend3,2*pi*f);
    Hc4 = freqs(numc4,dend4,2*pi*f);
    figure()
    hold on
    plot(f,20*log10(abs(Hc1)))
    plot(f,20*log10(abs(Hc2)))
    plot(f,20*log10(abs(Hc3)))
    plot(f,20*log10(abs(Hc4)))
    axis([fmin*0.1 fmax*2 -70 5])
    grid on
    legend('Butterworth',' Chebyshev de type 1' ,' Chebyshev de type 2' ,' Cauer');

    bbecdb31b3684b394afe5c1b8d0a97e9.png
    四种滤波器的频域响应
    • 除了抖动,还可以看到过渡带(从通带到阻带的部分)的宽度也不同,紫色和橙色下降地比较迅速,过渡带比较窄,蓝色和黄色下降地比较缓慢,过渡带比较宽,理论上来讲,过渡带越窄,滤波器越接近理想的滤波器模型,滤波性能越好。

    应用滤波器处理信号

    使用lsim命令得到滤波后的信号

    sigf1 = lsim(tf(numc1,dend1),sig,t);
    sigf2 = lsim(tf(numc2,dend2),sig,t);
    sigf3 = lsim(tf(numc3,dend3),sig,t);
    sigf4 = lsim(tf(numc4,dend4),sig,t);

    性能分析

    可以看一下得到的信号的频域与原信号有什么不同

    • 频域上
    figure()
    hold on
    subplot(3,2,1)
    FFTplot(sig,fe,fmin,fmax) %这里的FFT是手写函数,用于绘制信号的频谱,很简单,详细定义见文末
    subplot(3,2,3)
    FFTplot(sigf1,fe,fmin,fmax)
    subplot(3,2,4)
    FFTplot(sigf2,fe,fmin,fmax)
    subplot(3,2,5)
    FFTplot(sigf3,fe,fmin,fmax)
    subplot(3,2,6)
    FFTplot(sigf4,fe,fmin,fmax)

    f7743a47808d7de80a2f5f464a935eb3.png
    左上角:原信号,其余四个对应四种不同滤波器
    • 时域上
    figure()
    hold on
    subplot(3,2,1)
    plot(t,abs(sig))
    title('Raw signal')
    subplot(3,2,3)
    plot(t,abs(sigf1))
    title('Butterworth')
    subplot(3,2,4)
    plot(t,abs(sigf2))
    title('Chebyshev de type 1')
    subplot(3,2,5)
    plot(t,abs(sigf3))
    title('Chebyshev de type 2')
    subplot(3,2,6)
    plot(t,abs(sigf4))
    title('Cauer')

    bb0c4f763244610378d6575a28915035.png
    • 原始信号每秒依次发出Do,Re,Mi....幅值都是1,所以是个大方块
    • 频域和时域的图很像,这是因为这个例子的特殊性:时间上多个频率是独立的,即这段信号中,每个频率的成分都是依次出现的,不存在多个频率同时发生的情况。
    • 可以看到:
      • 第七个和第八个音都基本被很好的消除了,第六个音差一点
      • 我们设计阻带的增益是-40dB,根据Bode的计算公式,这意味着“滤波后的该频率的幅值是原先的
        倍”
      • 第五个音并不是我们想要消除的,但是都受到了不同程度的影响,这是因为我们的滤波器指标有误,其中的过渡带太宽了,导致第五个音受到了影响
      • 不同滤波器在通带的抖动情况也不一样,Butter和Cheby2在通带无抖(但依旧在第四个音有可见的衰减!),Cheby1和Cauer在通带有抖,这影响了最终得到的信号质量

    重新设计性能指标

    从上面的仿真可以看到,因为指标的错误,第五个音被意外地滤掉了一大部分,并且四个滤波器的通带表现也不好,要么是有抖,要么是有可见的衰减,为了得到更好的滤波效果,我们更改滤波器指标,让它更加严格。

    Ap = 1; %通带增益从-3到-1
    Aa = 40;
    Deltaf = 20; %过渡带带宽从100到20
    fc = 420;
    fa = fc + Deltaf/2;
    fp = fc - Deltaf/2;

    重新计算四种滤波器,可以得到性能更加好但阶数也更高的滤波器

    值得一提的是,此时Butterworth滤波器的系数计算结果会发散,我怀疑这意味Butterworth不能满足对应的频域指标

    我也把计算四种连续时间低通滤波器并画图的代码封装成了函数,只需要一行,就可以得到下图,函数放在了文末

    [numc1,denc1]=FiltrecoePBA(Ap,Aa,fc,Deltaf,'c1',1); %c1代表cheby1型,最后的1代表输出频域对比图

    20bb751ed70fcaa1012f8cfc66c53c8f.png

    新得到的这几种滤波器都能更好的完成滤波任务,这体现在:第五个音没有被滤掉,后三个音都被很好地滤掉。

    2ec46d541a7206f5ba46af74c08db65e.png
    上:原始信号,下:滤波后的信号

    手动封装的函数

    为了便于对比各个滤波器的性能,让程序更加简洁,封装了几个常用操作:

    • 输出信号傅里叶图谱
    • 输入预期指标,计算四种连续低通滤波器,并输出四种滤波器频域对比图
    • 输入信号和滤波器,输出滤波后信号,并显示滤波前后频域分析

    信号频域分析

    • 函数说明和例程
    Draw the Fourier transform of a signal.
    Inputs
    sig: signal
    fe: signal sampling frequency
    Optional input
    [min_ max _] represents the frequency domain interval.
    -------------------------------------------------------------------------
    example: %复制之后直接跑即可
        fe = 1000;            % Sampling frequency                    
        T = 1/fe;             % Sampling period       
        L = 7000;             % Length of signal
        t = (0:L-1)*T;        % Time vector
        sig = 0.7*sin(2*pi*50*t) + sin(2*pi*128*t) + 0.2*randn(size(t)); 
        FFTplot(sig,fe)
    

    例程中生成了一段带有噪声的两个正弦信号的叠加,输出是该段信号的频域图

    a8fa2e749c5b87f627136d06a329f33f.png
    • 函数正文
    function no_output = FFTplot(sig, fe,min_, max_)
    
    if nargin == 2
        min_ = 0 ;
        max_ = fe/2;
    elseif nargin < 2
        disp('input error!')
        return 
    end
    
    f = [0:length(sig)-1]*(fe/length(sig));
    S = fft(sig)/fe;
    plot(f, abs(S),'LineWidth',1.5);
    
    axis([min_ max_ 0 max(abs(S))*1.2])
    grid on
    
    end
    

    计算四种连续低通滤波器

    • 函数说明和例程
    Calculating Transfer Function Parameters for Four Continuous-Time Low-Pass Filters
    Input: 
            Ap,  ondulation bande passante  
            Aa,  attenuation bande attenuee
            fc, 
            Deltaf,
            type, 'b'/'c1'/'c2'/'e'
            plots = 1 : show bode plot  
                    else : without plot 
    output: num/den of filter
    -----------------------------------------------------------------------------------------------
    example: %复制之后直接跑即可
            Ap = 3;
            Aa = 40;
            Deltaf = 100;
            fc = 420;
            type = 'c1'
            plots = 1
            [numc,denc] = FiltrecoePBA(Ap,Aa,fc,Deltaf,type,plots); % calculate filter

    例程定义了对应的低通滤波器性能指标,调用该函数计算cheby1类型滤波器,返回滤波器分子分母参数,并画出四种滤波器频域响应

    e3bbb08b087605a7194e2c358e698aca.png
    • 函数正文
    function [numc,denc] = FiltrecoePBA(Ap,Aa,fc,Deltaf,type,plots)
    fa = fc + Deltaf/2;
    fp = fc - Deltaf/2;
    
    [nc1 wnc1]= buttord(2*pi*fp,2*pi*fa,Ap,Aa,'s');
    [nc2 wnc2]= cheb1ord(2*pi*fp,2*pi*fa,Ap,Aa,'s');
    [nc3 wnc3]= cheb2ord(2*pi*fp,2*pi*fa,Ap,Aa,'s');
    [nc4 wnc4]= ellipord(2*pi*fp,2*pi*fa,Ap,Aa,'s');
    
    [numc1, denc1] = butter(nc1,wnc1,'low','s');
    [numc2, denc2] = cheby1(nc2,Ap,wnc2,'low','s');
    [numc3, denc3] = cheby2(nc3,Aa,wnc3,'low','s');
    [numc4, denc4] = ellip(nc4,Ap,Aa,wnc4,'s');
    
    
    
    if type == 'b'
        numc_ = numc1;
        denc_ = denc1;
    elseif type == 'c1'
        numc_ = numc2;
        denc_ = denc2;
    elseif type == 'c2'
        numc_ = numc3;
        denc_ = denc3;
    elseif type == 'e'
        numc_ = numc4;
        denc_ = denc4;
    else 
        disp(type)
        disp('not define!')
        return 
    end
    
    if sum(isnan(numc_))>0 || sum(isnan(denc_))>0
        disp(['Error: filter of type ',type,' diverge!'])
        return 
    else 
        numc = numc_;
        denc = denc_;
    end
    
    
    
    
    if plots==1
        f = linspace(fc/10,fc*4,1000);
        Hc1 = freqs(numc1,denc1,2*pi*f);
        Hc2 = freqs(numc2,denc2,2*pi*f);
        Hc3 = freqs(numc3,denc3,2*pi*f);
        Hc4 = freqs(numc4,denc4,2*pi*f);
        figure()
        hold on
        plot(f,20*log10(abs(Hc1)))
        plot(f,20*log10(abs(Hc2)))
        plot(f,20*log10(abs(Hc3)))
        plot(f,20*log10(abs(Hc4)))
        axis([fc/10,fc*4 -2*Aa 5])
        grid on 
        legend('Butterworth',' Chabyshev type 1' ,' Chabyshev type 2' ,' Cauer');
        hold off
    end
    
    
    
    end

    输入信号和滤波器,输出滤波后信号,并显示滤波前后频域分析

    • 函数说明及例程
    Given a signal and a continuous filter, draw the spectrum of the raw signal 
    and the filtered spectrum, return the filtered signal
    Inputs
    sig: Signal
    fe: Sampling Frequency
    numc : the molecular coefficient of the continuous filter transfer function.
    denc: The denominator of the continuous filter transfer function.
    Optional input
    min_ max_ : Frequency Interval
    Ouput = filtered signal
    ----------------------------------------------------------------------------------------
    Example:  %复制之后直接跑即可
        fe = 1000;            % Sampling frequency                    
        T = 1/fe;              % Sampling period       
        L = 7000;             % Length of signal
        t = (0:L-1)*T;        % Time vector
        sig = 0.7*sin(2*pi*50*t) + sin(2*pi*128*t) + 0.2*randn(size(t));  % generate signal 
        Ap = 3; Aa = 40; Deltaf = 50; fc = 128; type = 'c1';, plots = 0;
        [numc,denc] = FiltrecoePBA(Ap,Aa,fc,Deltaf,type,plots); % calculate filter
        FiltreA(sig,fe,numc,denc,30,150);

    例程中生成了一段带有噪声的两个正弦信号的叠加,并生成了一个c1类型的低通滤波器,输出滤波后的信号,并显示滤波前后的频域对比

    84ceefe6f24cf5523c5be8584c1abfb1.png
    • 函数正文
    function [sigf] = FiltreA(sig, fe, numc, denc, min_, max_)
    
    t = [0:1/fe:1/fe*(length(sig)-1)];   
    sigf = lsim(tf(numc,denc),sig,t);
    
    subplot(2,1,2)
    FFTplot(sigf, fe,min_, max_)
    title('Signal filtered')
    grid on 
    subplot(2,1,1)
    FFTplot(sig, fe,min_, max_)
    grid on 
    title('Raw Signal')
    
    end
    展开全文
  • 应用Matlab对人体的心电信号进行滤波实验目的综合应用信号频谱分析和数字滤波器设计的知识,实现心电信号的滤波。加深理解信号时域和频域分析的物理概念,理解设计指标的工程概念,认识不同类型滤波器的特性和适用...
  • MATLAB 语音信号处理之显示频谱、滤波显示原始语音信号的特征频带(谱峰位置)设计低通滤波器语音信号进行处理 显示原始语音信号的特征频带(谱峰位置) 利用Windows自带的录音机录制了一个16s的语音文件,编写...
  • 数字信号处理课程设计 基于matlab的语音信号处理 摘要 利用所学习的数字信号处理知识设计了一个有趣的音效处理系统首先设计了几种不同的滤波器声音进行滤波处理分析了时域和频域的变化比较了经过滤波处理后的声音...
  • 这是之前使用MATLAB做的一个项目,来语音信号进行分析处理,设计语音信号的导入及播放操作,读取信号长度、画图及生成滤波等操作,MATLAB小白可零基础上手,只要跟着教程一步一步来就可以,非常简单。具体教程已...
  • (频域滤波要从傅里叶的推导开始,所以先跳过,从空域开始)今天废话多了点,因为内容不是很多,滤波的概念其实是频域概念,即对信号频率进行处理,高于或低于截止频率的将被干掉,或者带通带限,就有了高通滤波器,低...
  • 这是之前使用MATLAB做的一个项目,来语音信号进行分析处理,设计语音信号的导入及播放操作,读取信号长度、画图及生成滤波等操作,MATLAB小白可零基础上手,只要跟着教程一步一步来就可以,非常简单。具体教程已...
  • 1.基于matlab的快速Fourier变换part1 转自:http://blog.sina.com.cn/s/blog_a1d5b9ba0102vxj2.html一...FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。那其在实际应用中,有哪些用途呢?1.有些信号在...
  • ##MATLAB之时域及频域的乐器信号分析及处理 ...这是之前使用MATLAB做的一个项目,来语音信号进行分析处理,设计语音信号的导入及播放 操作,读取信号长度、画图及生成滤波等操作,MATLAB小白可零基础上
  • 语音信号的数据滤波和处理) 由于是大二被同学带着做的所以不是很熟悉,像请大神仔细的解说一下,用来复试的问答。 首先:fft是傅里叶变换的函数吗?为什么在y1=fft()之后...
  • 语音信号滤波去噪

    2012-10-17 19:00:55
    滤波前后语音信号波形的变化中,由于我们录制的语音信号噪声不大,所以观察并不明显,但在频域波形中,我们可以明显的看到设计的滤波器语音信号进行滤波处理,将噪声进行了滤除。此次滤波基本达到了要求,完成...
  • 小波分析由于能同时在时域和频域对信号进行分析,所以它能有效地实现对信号的去噪。本文利用MATLAB 将信号分为左右声道,利用小波函数、数字滤波器滤掉左右声道噪声成分。最后综合两个信号即得到理想音质的信号。
  • (频域滤波要从傅里叶的推导开始,所以先跳过,从空域开始)今天废话多了点,因为内容不是很多,滤波的概念其实是频域概念,即对信号频率进行处理,高于或低于截止频率的将被干掉,或者带通带限,就有了高通滤波器,低...
  • 1.基于matlab的快速Fourier变换part1 转自:http://blog.sina.com.cn/s/blog_a1d5b9ba0102vxj2.html一...FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。那其在实际应用中,有哪些用途呢?1.有些信号在...
  • 课题基于Matlab有噪音语音信号处理的设计与实现,综合运用数字信号处理的理论知识加噪声语音信号进行时域、频域分析和滤波。在设计实现的过程中,用巴特沃斯、切比雪夫和双线性变换法设计IIR数字滤波器,并利用...
  • 数字信号处理课程设计 基于 matlab 的语音信号处理 精彩文档 实用标准文案 摘要 利用所学习的数字信号处理知识 设计了一个有趣的音效处理系统 首先设计了几种不 同的滤波器声音进行滤波处理 分析了时域和频域的...
  • 根据给定的一段MIT-BIH心电信号(100号),信号进行时域和频域分析,得知其中噪声成分,分别用双线性变换法设计IIR滤波器和频率抽样法设计FIR滤波器,将心电信号中的干扰滤除,分析滤波前后时域及频域发生的变化...
  • 摘 要 介绍了利用 MATLAB 信号处理工具箱进行 FIR 滤波器设计的三种方法程序设计法 ...1 前言 数字滤波器是一种用来过滤时间离散信号的数字系统通过抽样数据进行数学处理来 达到频域滤波的目的根据其单位冲激响
  • 测力仪软件Dynoware里面可以信号进行零漂补偿、求平均值、滤波等操作,但想要查看信号的频域,进行频域分析,还得将数据导出到matlab中处理。 首先,选择Dynoware中的Export Data按钮,将力信号数据输出为txt...
  • MATLAB语音信号处理系统GUI设计

    千次阅读 2020-11-30 21:50:10
    要求1)语音信号的采集2)语音信号的频谱分析3)设计高通、低通数字滤波器和画出其频率响应4)用滤波器对信号进行滤波5)比较滤波前后语音信号的波形及频谱6)播放原始语音信号7)实现快放、慢...
  • [转]基于MATLAB的FIR滤波器设计与滤波

    万次阅读 2009-02-04 22:54:00
    摘 要 介绍了利用MATLAB信号处理工具箱进FIR滤波器设计...关键词 MATLAB,数字滤波器,有限冲激响应,窗函数,仿真1 前言 数字滤波器是一种用来过滤时间离散信号的数字系统,通过抽样数据进行数学处理来达到频域
  • 摘? 要 介绍了利用TLAB信号处理工具箱进行FI滤波器设计的三种方法... 数字滤波器是一种用来过滤时间离散信号的数字系统,通过抽样数据进行数学处理来达到频域滤波的目的根据其单位冲激响应函数的时域特性可分为两类:
  • 摘 要 介绍了利用MATLAB信号处理工具箱进行FIR滤波器设计的三种方法:程序设计法、FDATool设计法和SPTool设计法,... 数字滤波器是一种用来过滤时间离散信号的数字系统,通过抽样数据进行数学处理来达到频域滤波
  • 本设计先完成语音信号的...把原始语音信号、加噪语音信号和滤波后的信号进行时域变换和频域变换,画出它们的时域波形和频域波形图,从视觉角度比较分析滤波的效果。也可将这3类信号进行播放从听觉角度感受滤波的效果。
  • 首先通过调用matlab中函数读取一段音乐信号,再对此音乐信号分别加上高斯白噪声、单音频噪声、多音频噪声,之后通过双线性变化方法设计无限长数字脉冲响应低通滤波器,并分别所加不同噪声的音乐信号进行滤波,并...
  • 利用matlab对心电信号进行分析,具体为:读取,插值,滤波,时域波形和频域波形分析
  • 自适应信号处理论文-基于matlab的自适应...LMS算法是时域变换,DCT是频域变换,文章采用MATLAB相关函数实现了对信号变换的仿真,并对这两种算法进行了一定的对比。关键词MATLAB,LMS算法,DCT变换1、引言LMS算法是...

空空如也

空空如也

1 2 3 4
收藏数 73
精华内容 29
关键字:

matlab对信号进行频域滤波

matlab 订阅