精华内容
下载资源
问答
  • MATLAB高通滤波器程序

    2018-08-29 16:30:56
    本程序是基于MATLAB软件数字高通滤波器
  • 通过对理想滤波器的加矩形窗实现滤波器设计
  • 频域滤波器器的一些相关实验和滤波器的实现(包括理想低通,Butterworth低通,高斯低通,理想高通,Butterworth高通,高斯高通滤波器
  • mATlab自编理想低通滤波器 数字图像处理
  • 精选文档 DSP设计滤波器报告 一低通滤波器的设计 一实验目的掌握 IIR 数字低通滤波器的设计方法 二实验原理 1滤波器的分类 滤波器分两大类经典滤波器和现代滤波器 经典滤波器是假定输入信号 x( n) 中的有用成分和...
  • matlab自定义函数(滤波器部分),dsp运用到的函数,,AFD_BUTT.M,AFD_CHB1.M,AFD_CHB2.M,AFD_ELIP.M,BLACKMAN.M,CAS2DIR.M,CASFILTR.M等等
  • 一、简介我们知道,在一幅图像中,其低频成分对应者图像变化缓慢的部分,对应着图像大致的相貌和轮廓。而其高频成分则对应着图像变化剧烈的部分,对应着图像的细节(图像的噪声也属于高频成分)。...理想低通滤波器...

    一、简介

    我们知道,在一幅图像中,其低频成分对应者图像变化缓慢的部分,对应着图像大致的相貌和轮廓。而其高频成分则对应着图像变化剧烈的部分,对应着图像的细节(图像的噪声也属于高频成分)。

    低频滤波器,顾名思义,就是过滤掉或者大幅度衰减图像的高频成分,让图像的低频成分通过。低频滤波器可以平滑图像,虑去图像的噪声。而与此相反的高频滤波器,则是过滤低频成分,通过高频成分,可以达到锐化图像的目的。

    理想低通滤波器的滤波非常尖锐,而高斯低通滤波器的滤波则非常平滑。Butterworth低通滤波器则介于两者之间,当Butterworth低通滤波器的阶数较高时,接近于理想低通滤波器,阶数较低时,则接近于高斯低通滤波器。

    下面我们所说的都是对应于二维图像处理的情况。


    二、这三种滤波器的相关介绍

    在这三种低通滤波器的表达式中,我们都用D0来表示其截止频率。D(u, v)表示距离频率矩形中心的距离。


    我并没有说明振铃现象产生的原因,只是说明了什么情况下振铃现象较为明显。


    1. 理想低通滤波器(ILPF)

    理想低通滤波器在以原点为圆心、D0为半径的园内,通过所有的频率,而在圆外截断所有的频率。(圆心的频率最低,为变换的直流(dc)分量)。函数如下:



    可以看出,理想低通滤波器的过渡非常急剧,会产生振铃现象。


    2. Butterworth低通滤波器

    函数表达式如下(在有些书中,Butterworth的函数的平方才等于右边的表达式,在这里我们按照课本的写法,计算较为方便),其中n称为Butterworth低通滤波器的阶数:



    从滤波器的函数图中,我们可以看出其过渡没有理想低通滤波器那么剧烈,从图(c)中可以看出,阶数越高,滤波器的过度越剧烈,振铃现象将越明显。


    3. 高斯低通滤波器(GLPF)

    函数表达式如下:



    高斯滤波器的过度特性非常平坦,因此不会产生振铃现象。


    三、频域滤波步骤

    1. 基本步骤

    1. 给定一幅大小为m*n的图像f(x,y)。选择适当的填充参数P和Q,一般令P = 2m,Q = 2n。

    2. 对图像f(x, y)填充0,填充后得到图像大小为P*Q的图像fp(x, y)。

    3. 用(-1)^(x+y)乘以fp(x,y)将其移到变换中心(中心化)。

    4. 计算fp(x, y)的DFT,得到F(u,v)。 

    5. 生成一个实的,对称的滤波函数H(u, v),大小为P*Q,中心在(P/2, Q/2)处。然后相乘(矩阵点乘)得到G(u,v) = H(u,v)F(u,v)。

    6. 对G(u, v)反傅里叶变换,然后取实部,再乘以(-1)^(x+y)进行反中心变换最后得到gp(x,y)。

    7. 提取gp(x,y)左上角的m*n区域,对提取的部分进行标准化处理,得到最终的结果图像g(x,y)。


    2. 相关步骤说明

    1. 对图像进行0填充,得到大小为P*Q的图像,主要是为了避免在循环卷积中出现的缠绕错误。当两个矩阵大小相同时,P≥2m-1,Q≥2n-1时可以避免环绕错误。由于对于偶数尺寸的矩阵计算其傅里叶变换较快,因此P取2m,Q取2n。

    2. 图像乘以(-1)^(x+y),再对图像进行傅里叶变换可以得到将原点移到中心的傅里叶变换。这样,对于F(u,v)来说,中心的频率最低,四周的频率较高。每一点的值表示该频率对于的幅度。在matlab中,也可以不乘以(-1)^(x+y),直接对填充后的图像进行傅里叶变换,之后使用fftshift函数对其傅里叶变换进行中心化。得到的结果是一样的。

    3. 由于图像是一个实函数,所以其傅里叶变换是一个旋转对称的傅里叶变换。

    4. 关于对称中心。课本中指出,对于一个长度为M的一维序列,当M为偶数时,位置0和M/2呈现零的特性,当M为奇数时,只有位置0呈现零的特性。因此,在matlab中,由于矩阵下标是从1开始的,我们的P = 2m,Q = 2n。因此中心点的位置为(P/2+1, Q/2+1),即(m+1, n+1)。

    5. 标准化处理。对最终得到的图像每一点的值进行处理。使其范围变为[0, 255]的uint8。或者[0, 1.0]的double值。


    四、matlab代码实现Butterworth低通滤波器

    在使用matlab代码的实现过程中,对于这三种低通滤波器,只是在实验低通滤波器函数H(u,v)的代码中有部分不同,其他部分一致。因此,在下面中,只给出实现Butterworth低通滤波器的代码,不给出其他两种滤波器的代码。


    1.  Butterworth滤波器的代码如下:

    该函数为Bfilter,输入为需要进行Butterworth滤波的灰度图像,Butterworth滤波器的截止频率D0以及Butterworth滤波器的阶数n。输出为进行滤波之后的图像(图像的值已经归一化到[[0, 255])。

    function [image_out] = Bfilter(image_in, D0, N)
    % Butterworth滤波器,在频率域进行滤波
    % 输入为需要进行滤波的灰度图像,Butterworth滤波器的截止频率D0,阶数N
    % 输出为滤波之后的灰度图像
    
    [m, n] = size(image_in);
    P = 2 * m;
    Q = 2 * n;
    
    fp = zeros(P, Q);
    %对图像填充0,并且乘以(-1)^(x+y) 以移到变换中心
    for i = 1 : m
        for j = 1 : n
            fp(i, j) = double(image_in(i, j)) * (-1)^(i+j);
        end
    end
    % 对填充后的图像进行傅里叶变换
    F1 = fft2(fp);
    
    % 生成Butterworth滤波函数,中心在(m+1,n+1)
    Bw = zeros(P, Q);
    a = D0^(2 * N);
    for u = 1 : P
        for v = 1 : Q
            temp = (u-(m+1.0))^2 + (v-(n+1.0))^2;
            Bw(u, v) = 1 / (1 + (temp^N) / a);
        end
    end
    
    %进行滤波
    G = F1 .* Bw;
    
    % 反傅里叶变换
    gp = ifft2(G);
    
    % 处理得到的图像
    image_out = zeros(m, n, 'uint8');
    gp = real(gp);
    g = zeros(m, n);
    for i = 1 : m
        for j = 1 : n
            g(i, j) = gp(i, j) * (-1)^(i+j);
            
        end
    end
    mmax = max(g(:));
    mmin = min(g(:));
    range = mmax-mmin;
    for i = 1 : m
        for j = 1 : n
            image_out(i,j) = uint8(255 * (g(i, j)-mmin) / range);
        end
    end
    
    end

    2. 测试代码如下

    测试时Butterworth滤波器的阶数为2。截止频率分别为10,30,60,160,460。

    clear all;
    close all;
    clc;
    
    image1 = imread('3.bmp');
    
    image2 = Bfilter(image1, 10, 2);
    image3 = Bfilter(image1, 30, 2);
    image4 = Bfilter(image1, 60, 2);
    image5 = Bfilter(image1, 160, 2);
    image6 = Bfilter(image1, 460, 2);
    
    % 显示图像
    subplot(2,3,1), imshow(image1), title('原图像');
    subplot(2,3,2), imshow(image2), title('D0 = 10, n = 2');
    subplot(2,3,3), imshow(image3), title('D0 = 30, n = 2');
    subplot(2,3,4), imshow(image4), title('D0 = 60, n = 2');
    subplot(2,3,5), imshow(image5), title('D0 = 160, n = 2');
    subplot(2,3,6), imshow(image6), title('D0 = 460, n = 2');


    3. 运行结果如下,可以看出,与课本给出的结果一致。




    展开全文
  • 理想低通滤波器 (ILPF):简单地切断所有距离变换原点指定距离
  • 理想高通滤波 MATLAB

    2017-02-16 15:11:23
    用于图像的理想滤波器高通滤波
  • 常用滤波器Matlab程序设计

    千次阅读 多人点赞 2019-10-31 23:22:00
    常用滤波器Matlab程序设计 (低通滤波器、高通滤波器、带通滤波器、带阻滤波器) 以下四个滤波器都是切比雪夫I型数字滤波器 1.低通滤波器 ​ 低通滤波(Low-pass filter) 是一种过滤方式,规则为低频信号能正常通过,...

    常用滤波器Matlab程序设计

    (低通滤波器、高通滤波器、带通滤波器、带阻滤波器)

    以下四个滤波器都是切比雪夫I型数字滤波器

    1.低通滤波器

    ​ 低通滤波(Low-pass filter) 是一种过滤方式,规则为低频信号能正常通过,而超过设定临界值的高频信号则被阻隔、减弱。但是阻隔、减弱的幅度则会依据不同的频率以及不同的滤波程序(目的)而改变。它有的时候也被叫做高频去除过滤(high-cut filter)或者最高去除过滤(treble-cut filter)。低通过滤是高通过滤的对立。

    低通滤波

    ​ 低通滤波可以简单的认为:设定一个频率点,当信号频率高于这个频率时不能通过,在数字信号中,这个频率点也就是截止频率,当频域高于这个截止频率时,则全部赋值为0。因为在这一处理过程中,让低频信号全部通过,所以称为低通滤波。

    ​ 低通过滤的概念存在于各种不同的领域,诸如电子电路,数据平滑,声学阻挡,图像模糊等领域经常会用到。

    ​ 在数字图像处理领域,从频域看,低通滤波可以对图像进行平滑去噪处理。

    低通滤波器

    对于不同滤波器而言,每个频率的信号的减弱程度不同。当使用在音频应用时,它有时被称为高频剪切滤波器,或高音消除滤波器。

    低通滤波器概念有许多不同的形式,其中包括电子线路(如音频设备中使用的hiss 滤波器、平滑数据的数字算法、音障(acoustic barriers)、图像模糊处理等等,这两个工具都通过剔除短期波动、保留长期发展趋势提供了信号的平滑形式。低通滤波器在信号处理中的作用等同于其它领域如金融领域中移动平均数(moving average)所起的作用;低通滤波器有很多种。其中,最通用的就是巴特沃斯滤波器和切比雪夫滤波器。

    接下来用Matlab程序设计低通滤波器生成m文件如下:

    function y=lowp(x,f1,f3,rp,rs,Fs)
    %低通滤波
    %使用注意事项:通带或阻带的截止频率的选取范围是不能超过采样率的一半
    %即,f1,f3的值都要小于 Fs/2
    %x:需要带通滤波的序列
    % f 1:通带截止频率
    % f 3:阻带截止频率
    %rp:边带区衰减DB数设置
    %rs:截止区衰减DB数设置
    %Fs:序列x的采样频率
    wp=2*pi*f1/Fs;
    ws=2*pi*f3/Fs;
    % 设计切比雪夫滤波器;
    [n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs);
    [bz1,az1]=cheby1(n,rp,wp/pi);
    %查看设计滤波器的曲线
    [h,w]=freqz(bz1,az1,256,Fs);
    h=20*log10(abs(h));
    figure;plot(w,h);title('当前低通滤波器的通带曲线');grid on;
    y=filter(bz1,az1,x);%对序列x滤波后得到的序列y
    end
    

    为了方便起见,专门封装专门画频谱图的函数,函数名plot_fft。后边所有绘制有关函数频谱图将会直接调用,为避免函数调用出现问题,此函数生产m文件保存到程序运行当前文件夹下。以下是plot_fft函数设计实现。

    %画信号的幅频谱和功率谱
    function plot_fft(y,fs,style,varargin)
    %当style=1,画幅值谱;当style=2,画功率谱;当style=其他的,那么花幅值谱和功率谱
    %当style=1时,还可以多输入2个可选参数
    %可选输入参数是用来控制需要查看的频率段的
    %第一个是需要查看的频率段起点
    %第二个是需要查看的频率段的终点
    %其他style不具备可选输入参数,如果输入发生位置错误
    nfft=2^nextpow2(length(y));%找出大于y的个数的最大的2的指数值(自动进算最佳FFT步长nfft)
    %nfft=1024;%人为设置FFT的步长nfft
      y=y-mean(y);%去除直流分量
    y_ft=fft(y,nfft);%对y信号进行DFT,得到频率的幅值分布
    y_p=y_ft.*conj(y_ft)/nfft;%conj()函数是求y函数的共轭复数,实数的共轭复数是他本身。
    y_f=fs*(0:nfft/2-1)/nfft;%T变换后对应的频率的序列
    % y_p=y_ft.*conj(y_ft)/nfft;%conj()函数是求y函数的共轭复数,实数的共轭复数是他本身。
    if style==1
        if nargin==3
           %plot(y_f,abs(y_ft(1:nfft/2)));%论坛上画FFT的方法
           plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));%matlab的帮助里画FFT的方法
           ylabel('幅值');xlabel('频率');title('信号幅值谱');
        else
           f1=varargin{1};
           fn=varargin{2};
           ni=round(f1 * nfft/fs+1);
           na=round(fn * nfft/fs+1);
           plot(y_f(ni:na),abs(y_ft(ni:na)*2/nfft));
           ylabel('幅值');xlabel('频率');title('信号幅值谱');
        end
    elseif style==2
               plot(y_f,y_p(1:nfft/2));
               ylabel('功率谱密度');xlabel('频率');title('信号功率谱');
        else
           subplot(211);plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));
           ylabel('幅值');xlabel('频率');title('信号幅值谱');
           subplot(212);plot(y_f,y_p(1:nfft/2));
           ylabel('功率谱密度');xlabel('频率');title('信号功率谱');
    end
    end
    
    低通滤波器案例设计:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %低通滤波器案例设计
    clear all; clc; close all;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % rp=0.1;rs=20;%通带边衰减DB值和阻带边衰减DB值(根据实际情况设定参数此为测试默认值)
    % Fs=2000;%采样率
    fs=2000;%采样频率
    t=(1:fs)/fs;%采样时间
    ff1=100;%信号频率100,400
    ff2=400;
    x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t);%带测试的x序列
    figure;
    subplot(211);plot(t,x);%原序列图象。
    subplot(212);plot_fft(x,fs,1);%原图像对应频谱
    %低通测试
    % y=filter(bz1,az1,x);
    y=lowp(x,300,350,0.1,20,fs);%低通滤波器函数测试。
    figure;
    subplot(211);plot(t,y);%低通滤波器输出序列
    subplot(212);plot_fft(y,fs,1);%低通滤波器输出序列对应频谱
    %plot_fft()函数已在上面给出主需要调用即可。
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-DrKPScht-1572535339774)(G:\研究生\项目小组任务\程序设计\第三周任务\lowp_result.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-46kcurke-1572535339778)(G:\研究生\项目小组任务\程序设计\第三周任务\band_low.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-D8XG3Piu-1572535339778)(G:\研究生\项目小组任务\程序设计\第三周任务\lowp.bmp)]

    总结:

    ​ 通过简单测试,以上三幅图分别是滤波前后的时频图,滤波器的滤波特性曲线图。通过图可以看出低通滤波器成功留下了100Hz的低频成分而把不要的高频成分去除了。实现低通滤波器基本功能。

    2.高通滤波器

    ​ 高通滤波(high-pass filter) 是一种过滤方式,规则为高频信号能正常通过,而低于设定临界值的低频信号则被阻隔、减弱。但是阻隔、减弱的幅度则会依据不同的频率以及不同的滤波程序(目的)而改变。它有的时候也被叫做低频去除过滤(low-cut filter)。高通滤波是低通滤波的对立。

    高通滤波

    ​ 高通滤波是只对低于某一给定频率以下的频率成分有衰减作用,而允许这个截频以上的频率成分通过,并且没有相位移的滤波过程。主要用来消除低频噪声,也称低截止滤波器。

    ​ 高通滤波属于[频率域]滤波,它保留高频,抑制低频,是[图像锐化]的一种方式。

    高通滤波器

    ​ 高通滤波器是一种让某一频率以上的信号分量通过,而对该频率以下的信号分量大大抑制的电容、电感与电阻等器件的组合装置。其特性在时域及频域中可分别用冲激响应及频率响应描述。后者是用以频率为自变量的函数表示,一般情况下它是一个以复变量jω为自变量的的复变函数,以H(jω)表示。它的模H(ω)和幅角φ(ω)为角频率ω的函数,分别称为系统的“幅频响应”和“相频响应”,它分别代表激励源中不同频率的信号成分通过该系统时所遇到的幅度变化和相位变化。可以证明,系统的“频率响应”就是该系统“冲激响应”的傅里叶变换。当线性无源系统可以用一个N阶线性微分方程表示时,频率响应H(jω)为一个有理分式,它的分子和分母分别与微分方程的右边和左边相对应。

    接下来用Matlab程序设计高通滤波器生成m文件如下:

    function y=highp(x,f1,f3,rp,rs,Fs)
    %高通滤波
    %使用注意事项:通带或阻带的截止频率的选取范围是不能超过采样率的一半
    %即,f1,f3的值都要小于 Fs/2
    %x:需要带通滤波的序列
    % f 1:通带截止频率
    % f 2:阻带截止频率
    %rp:边带区衰减DB数设置
    %rs:截止区衰减DB数设置
    %Fs:序列x的采样频率
    wp=2*pi*f1/Fs;
    ws=2*pi*f3/Fs;
    % 设计切比雪夫滤波器;
    [n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs);
    [bz1,az1]=cheby1(n,rp,wp/pi,'high');
    %查看设计滤波器的曲线
    [h,w]=freqz(bz1,az1,256,Fs);
    h=20*log10(abs(h));
    figure;plot(w,h);title('当前高通滤波器的通带曲线');grid on;
    y=filter(bz1,az1,x);
    end
    
    高通滤波器案例设计:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %高通滤波器案例设计测试程序
    clear all; clc; close all;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fs=2000;%采样频率
    t=(1:fs)/fs;%采样时间
    ff1=100;%信号频率100,400
    ff2=400;
    x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t);%带测试的x序列
    figure;
    subplot(211);plot(t,x);%原序列图谱。
    subplot(212);plot_fft(x,fs,1);%原图像对应频谱
    %高通测试
    z=highp(x,350,300,0.1,20,fs);%高通滤波器函数测试。
    figure;
    subplot(211);plot(t,z);%高通滤波器输出序列图谱
    subplot(212);plot_fft(z,fs,1);%高通滤波器输出序列对应频谱
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ARL0tTr0-1572535339782)(G:\研究生\项目小组任务\程序设计\第三周任务\highp.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Gkr6YIt8-1572535339786)(G:\研究生\项目小组任务\程序设计\第三周任务\highp_band.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-tcCtzaxF-1572535339790)(G:\研究生\项目小组任务\程序设计\第三周任务\highp_result.bmp)]

    总结:

    ​ 通过简单测试,以上三幅图分别是滤波前后的时频图,滤波器的滤波特性曲线图。通过图可以看出高通滤波器成功留下了400Hz的高频成分而把不要的低频成分去除了。实现高通滤波器基本功能。

    3.带通滤波器

    ​ 带通滤波器(band-pass filter)是一个允许特定频段的波通过同时屏蔽其他频段的设备。比如RLC振荡回路就是一个模拟带通滤波器。

    定义:

    带通滤波器是指能通过某一频率范围内的频率分量、但将其他范围的频率分量衰减到极低水平的滤波器,与带阻滤波器的概念相对。一个模拟带通滤波器的例子是电阻-电感-电容电路(RLC circuit)。这些滤波器也可以用低通滤波器高通滤波器组合来产生。

    工作原理:

    ​ 一个理想的带通滤波器应该有一个完全平坦的通带,在通带内没有放大或者衰减,并且在通带之外所有频率都被完全衰减掉,另外,通带外的转换在极小的频率范围完成。

    ​ 实际上,并不存在理想的带通滤波器。滤波器并不能够将期望频率范围外的所有频率完全衰减掉,尤其是在所要的通带外还有一个被衰减但是没有被隔离的范围。这通常称为滤波器的滚降现象,并且使用每十倍频的衰减幅度的dB数来表示。通常,滤波器的设计尽量保证滚降范围越窄越好,这样滤波器的性能就与设计更加接近。然而,随着滚降范围越来越小,通带就变得不再平坦,开始出现“波纹”。这种现象在通带的边缘处尤其明显,这种效应称为吉布斯现象

    ​ 在频带较低的剪切频率f1和较高的剪切频率f2之间是共振频率,这里滤波器的增益最大,滤波器的带宽就是f2和f1之间的差值。

    接下来用Matlab程序设计带通滤波器生成m文件如下:

    function y=bandp(x,f1,f3,fsl,fsh,rp,rs,Fs)
    %带通滤波
    %使用注意事项:通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半
    %即,f1,f3,fs1,fsh,的值小于 Fs/2
    %x:需要带通滤波的序列
    % f 1:通带左边界
    % f 3:通带右边界
    % fs1:衰减截止左边界
    % fsh:衰变截止右边界
    %rp:边带区衰减DB数设置
    %rs:截止区衰减DB数设置
    %Fs:序列x的采样频率
    wp1=2*pi*f1/Fs;
    wp3=2*pi*f3/Fs;
    wsl=2*pi*fsl/Fs;
    wsh=2*pi*fsh/Fs;
    wp=[wp1 wp3];
    ws=[wsl wsh];
    % 设计切比雪夫滤波器;
    [n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs);
    [bz1,az1]=cheby1(n,rp,wp/pi);
    %查看设计滤波器的曲线
    [h,w]=freqz(bz1,az1,256,Fs);
    h=20*log10(abs(h));
    figure;plot(w,h);title('所设计滤波器的通带曲线');grid on;
    y=filter(bz1,az1,x);
    end
    
    带通滤波器案例设计:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %带通滤波器案例设计测试程序
    clear all; clc; close all;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % f1=300;f3=500;%通带截止频率上下限
    % fsl=200;fsh=600;%阻带截止频率上下限
    % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
    % Fs=2000;%采样率
    fs=2000;%采样频率
    t=(1:fs)/fs;%采样时间
    ff1=100;%信号频率100,400,700
    ff2=400;
    ff3=700;
    x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t)+sin(2*pi*ff3*t);%带测试的x序列
    figure;
    subplot(211);plot(t,x);%原序列图谱。
    subplot(212);plot_fft(x,fs,1);%原图像对应频谱
    % y=filter(bz1,az1,x);
    %带通测试
    y=bandp(x,300,500,200,600,0.1,30,fs);%带通滤波器函数测试。
    figure;
    subplot(211);plot(t,y);%带通滤波器输出序列图谱
    subplot(212);plot_fft(y,fs,1);%带通滤波器输出序列对应频谱
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xnyopwYK-1572535339794)(G:\研究生\项目小组任务\程序设计\第三周任务\bandp.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ZQvJWlPr-1572535339798)(G:\研究生\项目小组任务\程序设计\第三周任务\bandp_band.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-DSrXHmCQ-1572535339798)(G:\研究生\项目小组任务\程序设计\第三周任务\bandp_result.bmp)]

    总结:

    ​ 通过简单测试,以上三幅图分别是滤波前后的时频图,滤波器的滤波特性曲线图。通过图可以看出带通滤波器成功留下了400Hz的带通频带内成分而把不要的频带成分去除了。实现带通滤波器基本功能。

    4.带阻滤波器

    ​ 带阻滤波器(bandstop filters,简称BSF)是指能通过大多数频率分量、但将某些范围的频率分量衰减到极低水平的滤波器,与带通滤波器的概念相对。其中点阻滤波器(notch filter)是一种特殊的带阻滤波器,它的阻带范围极小,有着很高的Q值(Q Factor)。

    定义:

    ​ 在电路中将输入电压同时作用于低通滤波器和高通滤波器,再将两个电路的输出电压求和,就可以得到带阻滤波器。带阻滤波器一般分为腔体带阻滤波器和LC带阻滤波器。

    接下来用Matlab程序设计带阻滤波器生成m文件如下:

    function y=bands(x,f1,f3,fsl,fsh,rp,rs,Fs)
    %带阻滤波
    %使用注意事项:通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半
    %即,f1,f3,fs1,fsh,的值小于 Fs/2
    %x:需要带通滤波的序列
    % f 1:通带左边界
    % f 3:通带右边界
    % fs1:衰减截止左边界
    % fsh:衰变截止右边界
    %rp:边带区衰减DB数设置
    %rs:截止区衰减DB数设置
    %FS:序列x的采样频率
    wp1=2*pi*f1/Fs;
    wp3=2*pi*f3/Fs;
    wsl=2*pi*fsl/Fs;
    wsh=2*pi*fsh/Fs;
    wp=[wp1 wp3];
    ws=[wsl wsh];
    % 设计切比雪夫滤波器;
    [n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs);
    [bz1,az1]=cheby1(n,rp,wp/pi,'stop');
    %查看设计滤波器的曲线
    [h,w]=freqz(bz1,az1,256,Fs);
    h=20*log10(abs(h));
    figure;plot(w,h);title('所设计滤波器的通带曲线');grid on;
    y=filter(bz1,az1,x);
    end
    
    带阻滤波器案例设计:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %带阻滤波器案例设计测试程序
    clear all; clc; close all;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % f1=300;f3=500;%通带截止频率上下限
    % fsl=200;fsh=600;%阻带截止频率上下限
    % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
    % Fs=2000;%采样率
    fs=1000;%采样频率
    t=(1:fs)/fs;%采样时间
    ff1=100;%信号频率100,400,700
    ff2=150;
    ff3=200;
    x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t)+sin(2*pi*ff3*t);%带测试的x序列
    figure;
    subplot(211);plot(t,x);%原序列图谱。
    subplot(212);plot_fft(x,fs,1);%原图像对应频谱
    %带通测试
    y=bands(x,110,190,140,160,0.1,30,fs);%带通滤波器函数测试。
    figure;
    subplot(211);plot(t,y);%带通滤波器输出序列图谱
    subplot(212);plot_fft(y,fs,1);%带通滤波器输出序列对应频谱
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LBhilewO-1572535339802)(G:\研究生\项目小组任务\程序设计\第三周任务\bands.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1FiaUXpQ-1572535339806)(G:\研究生\项目小组任务\程序设计\第三周任务\bands_band.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LTmY1TFf-1572535339810)(G:\研究生\项目小组任务\程序设计\第三周任务\bands_result.bmp)]

    总结:

    ​ 通过简单测试,以上三幅图分别是滤波前后的时频图,滤波器的滤波特性曲线图。通过图可以看出带通滤波器成功去除了150Hz的带通频带内成分而把需要的频带成分留下来。实现带阻滤波器基本功能。

    展开全文
  • 要使用过滤器: 1)在变量img中创建一张图片2) 调用函数 filter 来创建与图像 'img' 大小相同的过滤器 在滤波器代码中指定截止频率 d0 的值。
  • MATLAB在频域实现理想滤波、Butterworth、Gaussian滤波器
  • 通过对理想滤波器的加hanning窗实现滤波器设计
  • 对于不同滤波器而言,每个频率的信号的强弱程度不同。当使用在音频应用时,它有时被称为高频剪切滤波器, 或高音消除滤波器。低通滤波器概念有许多不同的形式,其中包括电子线路(如音频设备中使用的hiss 滤波器)、...
  • 该程序为利用matlab进行图像处理的代码,主要针对傅里叶正反变换以及其频率域的呈现,以及理想(方形)低通滤波器、Butterworth滤波器的实现,代码费了不少心思,毕竟八百年没碰过matlab了,竟然写到一半开始拼...
  • 陷波滤波器matlab程序

    2014-11-06 15:54:31
    陷波器是一种特殊的带阻滤波器,其阻带在理想情况下只有一个频率点,因此也被称为点阻滤波器。这种滤波器主要用于消除某个特定频率的干扰,例如,在各种测量仪器和数据采集系统中用于消除电源干扰的工频陷波器。
  • 三、用频率抽样法设计一个理想低通滤波器,该低通滤波器的通带截止频率为3л/8,阻带截止频率为4л/8。 1)不加过渡点。 2)加一个过渡点。 3)加两个过渡点。 要求:掌握并了解在过渡带加过渡点对逼近滤波器的实际...
  • ------------------------- y = sinc_filter(x,Wn) y = sinc_filter(x,Wn,N) y = sinc_filter(x,Wn,N,dim) y = sinc_filter(x,Wn,[],dim) y = sinc_filter(x,Wn) 将接近理想的低通或带通砖墙滤波器应用于阵列 x,沿...
  • [Matlab]FIR滤波器设计:(基本窗函数FIR滤波器设计)

    万次阅读 多人点赞 2019-11-16 00:54:00
    [Matlab]FIR滤波器设计:(基本窗函数FIR滤波器设计) ​ IIR滤波器主要设计方法先设计一个模拟低通滤波器,然后把它转化为形式上的数字滤波器。但对于FIR滤波器来说,设计方法的关键要求之一就是保证线性相位条件。而...

    [Matlab]FIR滤波器设计:(基本窗函数FIR滤波器设计)

    ​ IIR滤波器主要设计方法先设计一个模拟低通滤波器,然后把它转化为形式上的数字滤波器。但对于FIR滤波器来说,设计方法的关键要求之一就是保证线性相位条件。而IIR滤波器的设计方法中只对幅值特性进行设计,因此无法保证相位。所以FIR滤波器的设计需要采用完全不同的方法。FIR滤波器的设计方法主要有窗函数法、频率采样法、切比雪夫逼近法等。

    ​ 由于理想滤波器在边界频率处不连续,故时域信号 h d ( n ) h_d(n) hd(n)一定是无限时宽的,也是非因果的。所以理想低通滤波器是无法实现的。如果实现一个具有理想线性相位特性的滤波器,其幅值特性只能采用逼近理想幅频特性的方法实现。如果对时域信号 h d ( n ) h_d(n) hd(n)进行截取,并保证截取过程中序列保持对称,而且截取长度为N,则对称点为 α = 1 2 ∗ ( N − 1 ) α=\frac{1}{2}*(N-1) α=21(N1)。截取后序列为 h ( n ) h(n) h(n),侧 h ( n ) h(n) h(n)可用下式子表示:
    h ( n ) = h d ( n ) ∗ w ( n ) h(n) = h_d(n)*w(n) h(n)=hd(n)w(n)
    式中, w ( n ) w(n) w(n)为截取函数,又称为穿函数。如果窗函数为矩形序列,则称之为矩形窗。窗函数有多种形式,为保证加窗后系统的线性相位特性,必须保证加窗后序列关于 α = 1 2 ∗ ( N − 1 ) α=\frac{1}{2}*(N-1) α=21(N1)点对称。常用的窗函数有矩形窗、汉宁窗、汉明窗、布莱克曼窗、凯塞窗。窗函数设计法的基本思想是用一个长度为N的序列 h ( n ) h(n) h(n)代替 h d ( n ) h_d(n) hd(n)作为实际设计的滤波器的单位脉冲响应。这种设计法成为窗函数设计法。显然在保证 h ( n ) h(n) h(n)对称性的前提下,窗函数长度N越长,则 h ( n ) h(n) h(n)越接近 h d ( n ) h_d(n) hd(n)。但是误差是肯定存在的,这种误差成为截断误差。

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-A0flZjTw-1573836815510)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\20170402112228645.jpg)]

    ​ 要确定如何设计一个FIR滤波器,首先得对加窗后的理想滤波器特性变化进行分析,并研究减少由截断引起的误差的途径,从而提出更好的滤波器设计方案。对于调整窗口长度可以有效地控制过渡带的宽度,但减少带内波动以及加大阻带衰减没有作用。所以必须挑选最为合适的窗函数对理想滤波器进行截取。下面简单介绍几种窗函数。一个实际的滤波器的单位脉冲响应可表示为:
    h ( n ) = h d ( n ) ∗ w ( n ) h(n) = h_d(n)*w(n) h(n)=hd(n)w(n)
    几种常见窗函数(只给出了窗函数的定义和幅度特性):
    W ( e j ∗ w ) = W ( w ) ( e − j α w ) W(e^{j*w})=W(w)(e^{-jαw}) W(ejw)=W(w)(ejαw)

    矩形窗FIR滤波器设计:

    矩形窗的窗函数为:
    w R ( n ) = R N ( n ) w_R(n)=R_N(n) wR(n)=RN(n)
    幅度函数为:
    R N ( w ) = s i n ( w N / 2 ) s i n ( w / 2 ) R_N(w) = \frac{sin(wN/2)}{sin(w/2)} RN(w)=sin(w/2)sin(wN/2)
    它的主瓣宽度为 4 π / N 4\pi/N 4π/N,第一瓣比主瓣地13dB.

    在Matlab中,实现矩形窗的函数为boxcar和recttwin ,其调用格式如下:

    w=boxcar(N);
    w=recttwin(N);
    %%%显示窗函数的GUI工具
    n = 60;
    w = rectwin(n);
    wvtool(w)%显示窗函数的GUI工具
    %还提供了显示窗函数的GUI工具,如wvtool可以显示用来显示窗的形状和频域图形,wintool可以打开窗设计和分析工具,如运行
    wvtool(hamming(64),hann(64),gausswin(64))
    %%可以对比汉明窗、汉宁窗和高斯窗
    

    其中,N是窗函数的长度,返回值w是一个N阶的向量,它的元素有窗函数的值组成。其中w=boxcar(N)等价于w=ones(N,1)。

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eknDhIxt-1573836815512)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\20170402112307695.png)]

    案例分析:

    利用矩形窗设计FIR带阻滤波器,主要参数如下:

    给定抽样频率为 Ω s = 2 ∗ p i ∗ 1.5 ∗ 1 0 4 ( r a d / s e c ) Ωs=2*pi*1.5*10^4(rad/sec) Ωs=2pi1.5104(rad/sec),

    通带截至频率为 Ω p 1 = 2 ∗ p i ∗ 0.75 ∗ 1 0 3 ( r a d / s e c ) , Ω p 2 = 2 ∗ p i ∗ 6 ∗ 1 0 3 ( r a d / s e c ) Ωp1=2*pi*0.75*10^3(rad/sec),Ωp2=2*pi*6*10^3(rad/sec) Ωp1=2pi0.75103(rad/sec),Ωp2=2pi6103(rad/sec)
    阻带截至频率为 Ω s t 1 = 2 ∗ p i ∗ 2.25 ∗ 1 0 3 ( r a d / s e c ) , Ω s t 2 = 2 ∗ p i ∗ 1.5 ∗ 1 0 3 ( r a d / s e c ) Ωst1=2*pi*2.25*10^3(rad/sec),Ωst2=2*pi*1.5*10^3(rad/sec) Ωst1=2pi2.25103(rad/sec),Ωst2=2pi1.5103(rad/sec)
    阻带衰减 δ 2 > = 50 d B {\delta}_2 >=50dB δ2>=50dB

    %%%%调用子程序1:
    function hd=ideal_bs(Wcl,Wch,m); 
    alpha=(m-1)/2; 
    n=[0:1:(m-1)];
    m=n-alpha+eps; 
    hd=[sin(m*pi)+sin(Wcl*m)-sin(Wch*m)]./(pi*m)
    %%%%调用子程序2:
    function[db,mag,pha,w]=freqz_m2(b,a)
    [H,w]=freqz(b,a,1000,'whole');
    H=(H(1:1:501))'; w=(w(1:1:501))';
    mag=abs(H);
    db=20*log10((mag+eps)/max(mag));
    pha=angle(H);
    %%%%运行MATLAB源代码如下:
    clear all;
    Wph=3*pi*6.25/15;%通带频率
    Wpl=3*pi/15;
    Wsl=3*pi*2.5/15;%阻带频率
    Wsh=3*pi*4.75/15;
    tr_width=min((Wsl-Wpl),(Wph-Wsh));%%过渡带带宽
    %过渡带宽度
    N=ceil(4*pi/tr_width);					%滤波器长度
    n=0:1:N-1;
    Wcl=(Wsl+Wpl)/2;						%理想滤波器的截止频率
    Wch=(Wsh+Wph)/2;
    hd=ideal_bs(Wcl,Wch,N);				%理想滤波器的单位冲击响应
    w_ham=(boxcar(N))';
    string=['矩形窗','N=',num2str(N)];
    h=hd.*w_ham;						%截取取得实际的单位脉冲响应
    [db,mag,pha,w]=freqz_m2(h,[1]);
    %计算实际滤波器的幅度响应
    delta_w=2*pi/1000;
    subplot(241);
    stem(n,hd);
    title('理想脉冲响应hd(n)')
    axis([-1,N,-0.5,0.8]);
    xlabel('n');ylabel('hd(n)');
    grid on
    subplot(242);
    stem(n,w_ham);
    axis([-1,N,0,1.1]);
    xlabel('n');ylabel('w(n)');
    text(1.5,1.3,string);
    grid on
    subplot(243);
    stem(n,h);title('实际脉冲响应h(n)');
    axis([0,N,-1.4,1.4]);
    xlabel('n');ylabel('h(n)');
    grid on
    subplot(244);
    plot(w,pha);title('相频特性');
    axis([0,3.15,-4,4]);
    xlabel('频率(rad)');ylabel('相位(Φ)');
    grid on
    subplot(245);
    plot(w/pi,db);title('幅度特性(dB)');
    axis([0,1,-80,10]);
    xlabel('频率(pi)');ylabel('分贝数');
    grid on
    subplot(246);
    plot(w,mag);title('频率特性')
    axis([0,3,0,2]);
    xlabel('频率(rad)');ylabel('幅值');
    grid on
    fs=15000;
    t=(0:100)/fs;
    x=cos(2*pi*t*750)+cos(2*pi*t*3000)+cos(2*pi*t*6100);
    q=filter(h,1,x);
    [a,f1]=freqz(x);
    f1=f1/pi*fs/2;
    [b,f2]=freqz(q);
    f2=f2/pi*fs/2;
    subplot(247);
    plot(f1,abs(a));
    title('输入波形频谱图');
    xlabel('频率');ylabel('幅度')
    grid on
    subplot(248);
    plot(f2,abs(b));
    title('输出波形频谱图');
    xlabel('频率');ylabel('幅度')
    grid on
    
    

    在这里插入图片描述

    汉宁窗FIR滤波器设计:

    汉宁窗(hanning window)又称升余弦窗,窗函数为:
    w H n ( n ) = 0.5 [ 1 − c o s ( 2 π ∗ n N − 1 ) ] ∗ R N ( n ) w_{Hn}(n) = 0.5[1-cos(\frac{2\pi*n}{N-1})]*R_N(n) wHn(n)=0.5[1cos(N12πn)]RN(n)
    幅值函数为:
    W H n ( w ) = 0.5 W R ( w ) + 0.25 [ W R ( w − 2 π N − 1 ) + W R ( w + 2 π N − 1 ) ] W_{Hn}(w) = 0.5W_R(w) + 0.25[W_R(w-\frac{2\pi}{N-1})+W_R(w+\frac{2\pi}{N-1})] WHn(w)=0.5WR(w)+0.25[WR(wN12π)+WR(w+N12π)]
    汉宁窗幅度函数由3部分相加而成,其结果是使主瓣集中了更多能量,而旁瓣3部分相加时相互抵消而变小,其代价是主瓣宽度增加到 8 π / N 8\pi/N 8π/N。第一瓣比主瓣低31dB,阻带衰减加大。

    在Matlab中,实现汉宁窗的函数为hanning和barthannwin ,其调用格式如下:

    w=hanning(N)
    w=barthannwin(N)
    

    案例1:绘制50个点的汉宁窗。

    N=49;n=1:N;
    wdhn=hanning(N);	
    figure(3);
    stem(n,wdhn,'.');
    grid on
    axis([0,N,0,1.1]);
    title('50点汉宁窗');
    ylabel('W(n)');
    xlabel('n');
    title('50点汉宁窗');
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-KhfUUgrt-1573836815515)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanning_50point.bmp)]

    案例2:已知连续信号为 x ( t ) = c o s ( 2 π f 1 t ) + 0.15 c o s ( 2 π f 2 t ) x(t)=cos(2{\pi}f_1t) +0.15cos(2{\pi}f_2t) x(t)=cos(2πf1t)+0.15cos(2πf2t),其中 f 1 = 100 H z , f 2 = 150 H z f_1=100Hz,f_2=150Hz f1=100Hz,f2=150Hz。若抽样频率 f s a m = 600 H z f_sam=600Hz fsam=600Hzd对信号进行抽样,利用不同宽度N的矩形截断该序列,N取40,观察不同的窗对普分析结果的影响。

    N=40; 
    L=512;
    f1=100;f2=150;fs=600;
    ws=2*pi*fs;
    t=(0:N-1)*(1/fs);
    x=cos(2*pi*f1*t)+0.25*sin (2*pi*f2*t);
    wh=boxcar(N)';
    x=x.*wh;
    subplot(221);stem(t,x);
    title('加矩形窗时域图');
    xlabel('n');ylabel('h(n)')
    grid on
    W=fft(x,L);
    f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi); 
    subplot(222);
    plot(f,abs(fftshift(W)))
    title('加矩形窗频域图');
    xlabel('频率');ylabel('幅度')
    grid on
    x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t); 
    wh=hanning(N)';
    x=x.*wh;
    subplot(223);stem(t,x);
    title('加汉宁窗时域图');
    xlabel('n');ylabel('h(n)')
    grid on
    W=fft(x,L);
    f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);
    subplot(224);
    plot(f,abs(fftshift(W)))
    title('加汉宁窗频域图');
    xlabel('频率');ylabel('幅度')
    grid on
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wc7oZ6mR-1573836815515)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\rectange_6.8.bmp)]

    用汉宁窗对谐波信号进行分析:

    clear; 
    % 原始数据:直流:0V; 基波:49.5Hz,100V,10deg; HR2:0.5V,40deg;
    hr0=0;f1=50.1; 
    hr(1)=25*sqrt(2);deg(1)=10; 
    hr(2)=0;deg(2)=0; 
    hr(3)=1.755*sqrt(2);deg(3)=40; 
    hr(4)=0;deg(4)=0; 
    hr(5)=0.885*sqrt(2);deg(5)=70; 
    hr(6)=0;deg(6)=0; 
    hr(7)=1.125;deg(7)=110; 
    M=7;f=[1:M]*f1;							%设定频率 
    % 采样 
    fs=10000; 
    N=2048;									% 约10个周期 
    T=1/fs; 
    n=[0:N-1];t=n*T; 
    x=zeros(size(t)); 
    for k=1:M 
        x=x+hr(k)*cos(2*pi*f(k)*t+deg(k)*pi/180); 
    end 
    %分析: 
    w=0.5-0.5*cos(2*pi*n/N);
    Xk=fft(x.*w); 
    amp=abs(Xk(1:N/2))/N*2;						%幅频 
    pha=angle(Xk(1:N/2))/pi*180;					%相频 
    for k=1:N/2 
        if(amp(k)<0.01) pha(k)=0;  %当谐波<10mV时,其相位=0 
        end 
        if(pha(k)<0) pha(k)=pha(k)+360;%调整到0-360度 
        end 
    end 
    fmin=fs/N; 
    xaxis=fmin*n(1:N/2);
    %横坐标为Hz 
    kx=round([1:M]*50/fmin);
    %各次谐波对应的下标(从0开始) 
    for m=1:M 
        km(m)=searchpeaks(amp,kx(m)+1);			%km为谱峰(从1开始) 
        if(amp(km(m)+1)<amp(km(m)-1)) 
            km(m)=km(m)-1; 
        end 
        beta(m)=amp(km(m)+1)./amp(km(m)); 
        delta(m)=(2*beta(m)-1)./(1+beta(m)); 
    end 
    fx=(km-1+delta)*fmin;							%估计频率 
    hrx=amp(km)*2.*pi.*delta.*(1-delta.*delta)./sin(pi*delta);
    degx=pha(km)-delta.*180/N*(N-1);				%估计相位 
    degx=mod(degx,360);							%调整到0-360度 
    efx=(fx-f)./f*100;								%频率误差 
    ehr=(hrx-hr)./hr*100;							%幅度误差 
    edeg=(degx-deg);								%相位误差 
    % 结果输出: 
    subplot(2,2,1);
    %画出采样序列 
    plot(t,x); 
    hold on; 
    plot(t,x.*w,'r');
    %加窗波形 
    hold off; 
    xlabel('x(k)'); 
    title('原信号和加窗信号 '); 
    subplot(2,2,2);
    %画出FFT分析结果 
    stem(xaxis,amp,'.r'); 
    xlabel('频率'); 
    title('幅频结果'); 
    subplot(2,2,4); 
    stem(xaxis,pha,'.r'); 
    xlabel('角频率'); 
    title('相频结果'); 
    subplot(2,2,3); 
    stem(ehr); 
    title('幅度误差(%)'); 
    %文本输出 
    fid=fopen('result.txt','w'); 
    fprintf(fid,'原始数据:f1=%6.1fHz, N=%.f,  fs=%.f \r\n\r\n',f1,N,fs); 
    fprintf(fid,'谐波次数      1      2      3      4      5      6     7\r\n'); 
    fprintf(fid,'设定频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',f); 
    fprintf(fid,'估计频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',fx); 
    fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',efx); 
    fprintf(fid,'设定幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hr); 
    fprintf(fid,'估计幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hrx); 
    fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',ehr); 
    fprintf(fid,'设定相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',deg); 
    fprintf(fid,'估计相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',degx); 
    fprintf(fid,'误差(度)  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n\r\n',edeg); 
    %其他数据 
    fprintf(fid,'谱峰位置理论值:\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',[1:M]*f1/fmin); 
    fprintf(fid,'谱峰位置估计值:\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',km-1+delta); 
    fprintf(fid,'误差(%%)\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',((km-1+delta)-[1:M]*f1/fmin)./([1:M]*f1/fmin)*100); 
    fprintf(fid,'delta     :\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',delta); 
    fclose(fid);
    
    %运行过程中的调用子程序:
    function index1=searchpeaks(x,index) 
    %在数组中寻找最大值对应的下标 
    %x为数组,index 为给定的下标(index不能取最前或最后2个下标),在前后2个数中(共5个数)查找最大值和紧邻的次最大值 
    % indexmax 返回两个谱峰位置中的前一个谱峰对应的下标 
    index1=index-2; 
    for k=-1:2 
        if(x(index+k)>x(index1)) 
            index1=index+k; 
        end 
    end 
    if x(index1-1)>x(index1+1) 
        index1=index1-1; 
    end
    
    
    
    #result.txt输出结果
    原始数据:f1=  50.1Hz, N=2048,  fs=10000 
    
    谐波次数      1      2      3      4      5      6     7
    设定频率 50.100 100.200 150.300 200.400 250.500 300.600 350.700
    估计频率 50.100 78.819 150.302 181.252 250.499 279.138 350.701
    误差(%)  -0.000 -21.338  0.001 -9.555 -0.000 -7.140  0.000
    
    设定幅值 35.355  0.000  2.482  0.000  1.252  0.000  1.125
    估计幅值 35.356  0.046  2.482  0.002  1.252  0.002  1.125
    误差(%)   0.001    Inf  0.009    Inf  0.004    Inf  0.003
    
    设定相位  10.00   0.00  40.00   0.00  70.00   0.00 110.00
    估计相位  10.03  31.67  39.97 338.35  70.06 329.86 110.05
    误差()    0.03  31.67  -0.03 338.35   0.06 329.86   0.05
    
    谱峰位置理论值:
     10.2605 20.5210 30.7814 41.0419 51.3024 61.5629 71.8234
    谱峰位置估计值:
     10.2605 16.1421 30.7818 37.1203 51.3022 57.1675 71.8235
    误差(%)
     -0.0002 -21.3385 0.0012 -9.5551 -0.0004 -7.1396 0.0002
    delta     :
     0.2605 0.1421 0.7818 0.1203 0.3022 0.1675 0.8235
    
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wL0WYF6F-1573836815517)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanning_han.bmp)]

    汉明窗FIR滤波器设计:

    汉宁窗(hanming window)又称改进升余弦窗,窗函数为:
    w H n ( n ) = [ 0.54 − 0.46 ∗ c o s ( 2 π ∗ n N − 1 ) ] ∗ R N ( n ) w_{Hn}(n) = [0.54-0.46*cos(\frac{2\pi*n}{N-1})]*R_N(n) wHn(n)=[0.540.46cos(N12πn)]RN(n)
    幅值函数为:
    W H n ( w ) = 0.54 W R ( w ) + 0.23 [ W R ( w − 2 π N − 1 ) + 0.23 W R ( w + 2 π N − 1 ) ] W_{Hn}(w) = 0.54W_R(w) + 0.23[W_R(w-\frac{2\pi}{N-1})+0.23W_R(w+\frac{2\pi}{N-1})] WHn(w)=0.54WR(w)+0.23[WR(wN12π)+0.23WR(w+N12π)]
    汉明窗主瓣宽度与汉宁窗相同, 8 π / N , 99.96 % 8{\pi}/N,99.96\% 8π/N99.96%的能量集中在主瓣,第一瓣比主瓣低41dB。

    在Matlab中,实现汉宁窗的函数为hanming ,其调用格式如下:

    w=hanming(N)
    

    案例分析1:设计一个汉明窗低通滤波器:

    %语音信号设计一个汉明窗低通滤波器:
    [x,FS,bits]=wavread('C:\Windows\Media\Windows Ringout');
    x=x(:,1);
    figure(1);
    subplot(211);plot(x);
    title('语音信号时域波形图')
    xlabel('n');ylabel('h(n)')
    grid on
    y=fft(x,1000);
    f=(FS/1000)*[1:1000];  
    subplot(212);
    plot(f(1:300),abs(y(1:300)));
    title('语音信号频谱图');
    xlabel('频率');ylabel('幅度')
    grid on
    %产生噪声信号并加到语音信号
    t=0:length(x)-1;
    zs0=0.05*cos(2*pi*10000*t/1024);
    zs=[zeros(0,20000),zs0];
    figure(2);
    subplot(211)
    plot(zs)
    title('噪声信号波形');
    xlabel('n');ylabel('h(n)')
    grid on
    zs1=fft(zs,1200);
    subplot(212)
    plot(f(1:600),abs(zs1(1:600)));
    title('噪声信号频谱');
    xlabel('频率');ylabel('幅度')
    grid on
    x1=x+zs';
    %sound(x1,FS,bits);
    y1=fft(x1,1200);
    figure(3);
    subplot(211);plot(x1);
    title('加入噪声后的信号波形');
    xlabel('n');ylabel('h(n)')
    grid on
    subplot(212);
    plot(f(1:600),abs(y1(1:600)));
    title('加入噪声后的信号频谱');
    xlabel('频率');ylabel('幅度')
    grid on
    %滤波
    fp=7500;
    fc=8500; 
    wp=2*pi*fp/FS;
    ws=2*pi*fc/FS;
    Bt=ws-wp; 
    N0=ceil(6.2*pi/Bt);     
    N=N0+mod(N0+1,2);
    wc=(wp+ws)/2/pi;         
    hn=fir1(N-1,wc,hamming(N)); 
    X=conv(hn,x);           
    X1=fft(X,1200);
    figure(4);
    subplot(211);
    plot(X);
    title('滤波后的信号波形');
    xlabel('n');ylabel('h(n)')
    grid on
    subplot(212);
    plot(f(1:600),abs(X1(1:600))); 
    title('滤波后的信号频谱')
    xlabel('频率');ylabel('幅度')
    grid on
    
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-IRdanKny-1573836815518)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_1.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zpDvFiHo-1573836815519)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_noise.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bBaOdODs-1573836815521)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_addnoise.bmp)]

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-iBooX8kL-1573836815523)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_hanfilter.bmp)]

    案例分析2:已知连续信号为 x a ( t ) = c o s ( 100 π t ) + s i n ( 100 π t ) + c o s ( 50 π t ) x_a(t)=cos(100{\pi}t) +sin(100{\pi}t)+cos(50{\pi}t) xa(t)=cos(100πt)+sin(100πt)+cos(50πt),用DFT分析其中 x a ( t ) x_a(t) xa(t)的频谱结构,选择不同的截取长度 T p T_p Tp。观察存在的截断效应,试用加窗的方法减少谱间干扰。

    clear;close all
    fs=400;T=1/fs;						%采样频率和采样间隔
    Tp=0.04;N=Tp*fs;					%采样点数N
    N1=[N,4*N,8*N];					%设定三种截取长度 
    for m=1:3
        n=1:N1(m);
        xn=cos(100*pi*n*T)+ sin(200*pi*n*T)+ cos(50*pi*n*T);
        Xk=fft(xn,4096);
    fk=[0:4095]/4096/T;
    subplot(3,2,2*m-1);plot(fk,abs(Xk)/max(abs(Xk)));
    if m==1 title('矩形窗截取');
    end
    end
    %hamming窗截断
    for m=1:3
        n=1:N1(m);
        wn=hamming(N1(m));
        xn=cos(200*pi*n*T)+ sin(100*pi*n*T)+ cos(50*pi*n*T).*wn';
        Xk=fft(xn,4096);
    fk=[0:4095]/4096/T;
    subplot(3,2,2*m);plot(fk,abs(Xk)/max(abs(Xk)));
    if m==1 title('hamming窗截取');
    end
    end
    
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-39eMxIRF-1573836815524)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_dft.bmp)]

    比较矩形窗与汉明窗的普分析结果可见,矩形窗比用汉明窗分辨率高(泄露小),但是谱间干扰大。因此汉明窗是以牺牲分辨率来换取谱间干扰降低。

    布莱克曼FIR滤波器设计:

    布莱克曼(Blackman window)的窗函数为:
    w B l ( n ) = [ 0.42 − 0.5 ∗ c o s ( 2 π ∗ n N − 1 ) + 0.08 ∗ c o s ( 4 π ∗ n N − 1 ) ] ∗ R N ( n ) w_{Bl}(n) = [0.42-0.5*cos(\frac{2\pi*n}{N-1})+0.08*cos(\frac{4\pi*n}{N-1})]*R_N(n) wBl(n)=[0.420.5cos(N12πn)+0.08cos(N14πn)]RN(n)
    幅值函数为:
    W H n ( w ) = 0.42 W R ( w ) + 0.25 [ W R ( w − 2 π N − 1 ) + W R ( w + 2 π N − 1 ) ] W_{Hn}(w) = 0.42W_R(w) + 0.25[W_R(w-\frac{2\pi}{N-1})+W_R(w+\frac{2\pi}{N-1})] WHn(w)=0.42WR(w)+0.25[WR(wN12π)+WR(w+N12π)]

    + 0.04 [ W R ( w − 4 π N − 1 + W R ( w + 4 π N − 1 ) ] +0.04[W_R(w-\frac{4\pi}{N-1}+W_R(w+\frac{4\pi}{N-1})] +0.04[WR(wN14π+WR(w+N14π)]

    布莱克曼窗幅度函数由5部分相加而成,5部分相加的结果使得旁瓣得到进一步抵消,阻带衰减加大而过渡带加大到 12 π / N 12{\pi}/N 12π/N

    在Matlab中,实现布莱克曼窗的函数为blackman ,其调用格式如下:

    w=blackman(N);
    

    :案例:用窗函数法设计数字带通滤波器。下阻带边缘: W s 1 = 0.3 π , A s = 65 d B W_{s1}=0.3{\pi},A_s=65dB Ws1=0.3πAs=65dB,下通带边缘: W p 1 = 0.4 π , R p = 1 d B W_{p1}=0.4{\pi},R_p=1dB Wp1=0.4πRp=1dB,上通带边缘: W p 2 = 0.6 π , R p = 1 d B W_{p2}=0.6{\pi},R_p=1dB Wp2=0.6πRp=1dB,上阻带边缘: W s 2 = 0.7 π , R p = 65 d B W_{s2}=0.7{\pi},R_p=65dB Ws2=0.7πRp=65dB。根据窗函数最小阻带衰减的特性,以及参照窗函数的基本参数表,选择布莱克曼窗可以达到75dB的最小阻带衰减,其过渡带为 11 π / N 11\pi/N 11π/N

    clear all;
    wp1=0.4*pi;
    wp2=0.6*pi;
    ws1=0.3*pi;
    ws2=0.7*pi;
    As=65;
    tr_width=min((wp1-ws1),(ws2-wp2)); 					%过渡带宽度 
    M=ceil(11*pi/tr_width)+1							%滤波器长度
    n=[0:1:M-1];
    wc1=(ws1+wp1)/2;									%理想带通滤波器的下截止频率
    wc2=(ws2+wp2)/2;									%理想带通滤波器的上截止频率
    hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);
    w_bla=(blackman(M))';								%布莱克曼窗
    h=hd.*w_bla;
    %截取得到实际的单位脉冲响应
    [db,mag,pha,grd,w]=freqz_m(h,[1]);
    %计算实际滤波器的幅度响应
    delta_w=2*pi/1000;
    Rp=-min(db(wp1/delta_w+1:1:wp2/delta_w))
    %实际通带纹波
    As=-round(max(db(ws2/delta_w+1:1:501)))
    As=75
    subplot(2,2,1);
    stem(n,hd);
    title('理想单位脉冲响应hd(n)')
    axis([0 M-1 -0.4 0.5]);
    xlabel('n');
    ylabel('hd(n)')
    grid on;
    subplot(2,2,2);
    stem(n,w_bla);
    title('布莱克曼窗w(n)')
    axis([0 M-1 0 1.1]);
    xlabel('n');
    ylabel('w(n)')
    grid on;
    subplot(2,2,3);
    stem(n,h);
    title('实际单位脉冲响应hd(n)')
    axis([0 M-1 -0.4 0.5]);
    xlabel('n');
    ylabel('h(n)')
    grid on;
    subplot(2,2,4);
    plot(w/pi,db);
    axis([0 1 -150 10]);
    title('幅度响应(dB)');
    grid on;
    xlabel('频率单位:pi');
    ylabel('分贝数')
    
    %调用小程序设计1:
    function [db,mag,pha,grd,w] = freqz_m(b,a);
    % Modified version of freqz subroutine
    % ------------------------------------
    % [db,mag,pha,grd,w] = freqz_m(b,a);
    %  db = Relative magnitude in dB computed over 0 to pi radians
    % mag = absolute magnitude computed over 0 to pi radians 
    % pha = Phase response in radians over 0 to pi radians
    % grd = Group delay over 0 to pi radians
    %   w = 501 frequency samples between 0 to pi radians
    %   b = numerator polynomial of H(z)   (for FIR: b=h)
    %   a = denominator polynomial of H(z) (for FIR: a=[1])
    %
    [H,w] = freqz(b,a,1000,'whole');
        H = (H(1:1:501))'; w = (w(1:1:501))';
      mag = abs(H);
       db = 20*log10((mag+eps)/max(mag));
      pha = angle(H);
    %  pha = unwrap(angle(H));
      grd = grpdelay(b,a,w);
    %  grd = diff(pha);
    %  grd = [grd(1) grd];
    %  grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];
    %  grd = median(grd)*500/pi;
    
    %调用小程序设计2:
    function hd=ideal_lp(wc,M);
    %计算理想低通滤波器的脉冲响应
    %[hd]=ideal_lp(wc,M)
    %hd=理想脉冲响应0到M-1
    %wc=截止频率
    % M=理想滤波器的长度
    alpha=(M-1)/2;
    n=[0:1:(M-1)];
    m=n-alpha+eps;
    %加上一个很小的值eps避免除以0的错误情况出现
    hd=sin(wc*m)./(pi*m);
    
    
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-52LsWSri-1573836815526)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\blackmen_bandp.bmp)]

    凯塞窗FIR滤波器设计:

    凯塞-贝塞窗(Kaiser-Basel window)的窗函数为:
    w k ( n ) = I 0 ( β ) I 0 ( α ) ∗ R N ( n ) w_{k}(n) = \frac{I_0( β )}{I_0( \alpha )}*R_N(n) wk(n)=I0(α)I0(β)RN(n)
    式中, β = α ( 1 − ( 2 n N − 1 ) − 1 ) 2 \beta= {\alpha}{\sqrt{(1-(\frac{2n}{N-1})-1)^2}} β=α(1(N12n)1)2 I 0 ( x ) I_0( x ) I0(x)是零阶第一类修正贝塞函数,可用下面级数计算:
    I 0 ( x ) = 1 + ∑ k = 1 + ∞ ( 1 k ! ( x 2 ) k ) 2 I_0( x ) = 1+\sum_{k=1}^{+ ∞}(\frac{1}{k!}({\frac{x}{2}})^{k})^{2} I0(x)=1+k=1+k!1(2x)k2
    I 0 ( x ) I_0( x ) I0(x)q取15-25项就可以满足精度要求。通常 α \alpha α用以控制窗的形状, α \alpha α加大,主瓣加宽,旁瓣减小,典型数据 4 < α < 9 4<\alpha<9 4<α<9。当 α = 5.44 \alpha=5.44 α=5.44s时,窗函数接近汉明窗;当 α = 7.865 \alpha=7.865 α=7.865s时,窗函数接近于布莱克曼窗。其幅值函数为:
    W k ( w ) = w k ( 0 ) + 2 ∑ n = 1 ( N − 1 ) 2 ( w k ( n ) c o s ( w n ) ) W_{k}(w) =w_k(0) + 2\sum_{n=1}^{\frac{(N-1)}{2}}(w_k(n)cos(wn)) Wk(w)=wk(0)+2n=12(N1)(wk(n)cos(wn))
    在Matlab中,实现汉宁窗的函数为kaiser,其调用格式如下:

    w=kaiser(N);

    在Matlab中设计标准响应FIR滤波器可使用fir1函数。fir1函数以经典方法实现加窗性相位FIR滤波器的设计,它可以设计出标准的低通、高通、带通、带阻滤波器。fir1函数用法为:

    b = fir1(n,Wn,‘ftype’,wimdow)

    各个参数的含义如下:

    • b -滤波器系数,n-滤波器阶数。
    • Wn -截止频率, 0 < = W n < = 1 0<=W_n<=1 0<=Wn<=1, W n = 1 W_n=1 Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时, W n = [ W 1 , W 2 ] , W 1 < w < W 2 W_n =[W_1,W_2],W_1<w<W_2 Wn=[W1,W2],W1<w<W2
    • ftype -当指定ftype时,可设计高通和带阻滤波器。ftype=hight时,设计高通FIR滤波器;ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。
    • window–窗函数。窗函数的长度应等于FIR滤波器系数的个数,即阶数n+1。

    案例分析:利用凯塞窗函数设计一个带通滤波器,上截止频率2500Hz,下截止频率1000Hz,过渡带宽200Hz,通带纹波允许差0.1,带阻纹波不大于允差0.02dB,通带幅值为1。

    Fs=8000;N=216;
    fcuts=[1000 1200 2300 2500];
    mags=[0 1 0];
    devs=[0.02 0.1 0.02];
    [n,Wn,beta,ftype]=kaiserord(fcuts,mags,devs,Fs);
    n=n+rem(n,2);
    hh=fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
    [H,f]=freqz(hh,1,N,Fs);
    plot(f,abs(H));
    xlabel('频率 (Hz)');
    ylabel('幅值|H(f)|');
    grid on;
    
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zSlKyrYh-1573836815527)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\kaser_bandp.bmp)]

    窗函数设计法:

    根据前面几节的分析:设计一个FIR低通滤波器通常按照下面的步骤执行:

    1. 根据滤波器设计要求,确定滤波器的过渡带宽和阻带衰减要求,选择合适窗函数的类型并进行估计窗函数的宽度N。
    2. 根据所求的理想滤波器求出单位脉冲响应 h d ( n ) h_d(n) hd(n)
    3. 根据求得的h(n)求出其频率响应。
    4. 根据频率响应验证是否满足技术指标。
    5. 若不满足指标要求,则应调整窗函数类型或者长度,然后重复(1),(2),(3)(4)步,直到满足要求为止。

    注意:matlab中数据通常是以列向量形式存储的,所以两个向量相乘必须进行转置。计算滤波器的单位脉冲响应h(n),根据窗函数设计理论 h ( n ) = h d ( n ) ∗ w ( n ) h(n)=h_d(n)*w(n) h(n)=hd(n)w(n),在matlab中用语句hn=hd*wd实现h(n).

    窗函数设计法程序设计如下:

    function [h]=usefir1(mode,n,fp,fs,window,r,sample)
    % mode:模式(1--高通; 2--低通; 3--带通; 4--带阻)
    % n:阶数, 加窗的点数为阶数加1
    % fp:高通和低通时指示截止频率, 带通和带阻时指示下限频率
    % fs:带通和带阻时指示上限频率
    % window:加窗(1--矩形窗; 2--三角窗; 3--巴特窗; 4--汉明窗; 
    %5--汉宁窗; 6--布莱克曼窗; 7--凯泽窗; 8--契比雪夫窗)
    % r代表加chebyshev窗的r值和加kaiser窗时的beta值
    % sample:采样率
    % h:返回设计好的FIR滤波器系数
    if window==1 w=boxcar(n+1);
    end
    if window==2 w=triang(n+1);end
    if window==3 w=bartlett(n+1);end
    if window==4 w=hamming(n+1);end
    if window==5 w=hanning(n+1);end
    if window==6 w=blackman(n+1);end
    if window==7 w=kaiser(n+1,r);end
    if window==8 w=chebwin(n+1,r);
    end
    wp=2*fp/sample;
    ws=2*fs/sample;
    if mode==1 h=fir1(n,wp,'high',w);
    end
    if mode==2 h=fir1(n,wp,'low',w);
    end
    if mode==3 h=fir1(n,[wp,ws],w);
    end
    if mode==4 h=fir1(n,[wp,ws],'stop',w);
    end
    m=0:n;
    subplot(131);
    plot(m,h);grid on;	
    title('冲激响应');
    axis([0 n 1.1*min(h) 1.1*max(h)]);
    ylabel('h(n)');xlabel('n');
    freq_response=freqz(h,1);
    magnitude=20*log10(abs(freq_response));
    m=0:511; f=m*sample/(2*511);
    subplot(132);
    plot(f,magnitude);grid on;
    title('幅频特性');
    axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)]);
    ylabel('f幅值');xlabel('频率');
    phase=angle(freq_response);
    subplot(133);plot(f,phase);grid on;
    title('相频特性');
    axis([0 sample/2 1.1*min(phase) 1.1*max(phase)]);
    ylabel('相位');xlabel('频率');
    
    

    案例分析:假设需要设计一个40阶的带通FIR滤波器,采用汉明窗,采样频率为10kHz,两个截止频率分别为2kHz和3kHz,则需要在Matlab的命令行窗口输入:

    h=usefir1(3,60,2000,3000,4,2,10000);
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8Mmw8ohg-1573836815529)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\usefir1.bmp)]

    展开全文
  • 通过对理想滤波器加海明窗实现滤波器的设计
  • %%%%%%%=================理想低通滤波器(Ideal Lowpass)==================== %定义圆的半径,也可以理解为截断频率 r=160; %沿图像中心取出(2n+1)*(2n+1)的矩形 filter_circle=J_shift(M0-r:M0+r,N0-r:N0+r); %...
    clc;
    clear all;
    close all;
    
    %读入图像,并转换为double型
    I=imread('D:\Gray Files\4-41.tif');
    I_D=im2double(I);
    %获得图像的高度和宽度
    [M,N]=size(I_D);
    %图像中心点
    M0=M/2;
    N0=N/2;
    J=fft2(I_D);
    J_shift=fftshift(J);
    
    %%%%%%%=================理想低通滤波器(Ideal Lowpass)====================
    %定义圆的半径,也可以理解为截断频率
    r=160;
    %沿图像中心取出(2n+1)*(2n+1)的矩形
    filter_circle=J_shift(M0-r:M0+r,N0-r:N0+r);
    %计算是否在半径为n的圆内,圆点为n+1,利用二象限的值,确定一、三、四象限
    n2=r^2;
    %获取矩阵的宽度
    w=length(filter_circle);
    for x=1:r+1
        x2=(x-r+1)^2;
        for y=1:r+1
            y2=(y-r+1)^2;
            if x2+y2>n2
               filter_circle(x,y)=0;
               filter_circle(w-x+1,y)=0;
               filter_circle(w-x+1,w-y+1)=0;
               filter_circle(x,w-y+1)=0;
            end
        end
    end
    J_shift=zeros(M,N);
    J_shift(M0-r:M0+r,N0-r:N0+r)=filter_circle;
    
    J=ifftshift(J_shift);
    I_D_rep=ifft2(J);
    imshow(I_D_rep,[]);

     

    展开全文
  • matlab程序 理想低通滤波器

    热门讨论 2009-11-28 22:50:07
    Ideal Low Pass Filter using Discrete Fourier Transform.
  • %%%%%%%=================理想高通滤波器(Ideal Lowpass)==================== %定义圆的半径,也可以理解为截断频率 r=30; %沿图像中心取出(2n+1)*(2n+1)的矩形 filter_circle=J_shift(M0-r:M0+r,N0-r:N0+r); %...
  • 理想低通滤波器 f = imread('Fig0441.tif'); f = im2double(f); % 计算填充图像大小 [M,N] = size(f); M2 = 2*M; N2 = 2*N; % 傅里叶变换 F = fftshift(fft2(f,M2,N2)); figure;imshow(mat2gray(log(1+abs(F))))...
  • [Matlab]维纳滤波器设计

    千次阅读 2019-11-30 16:42:26
    [Matlab]维纳滤波器设计 ​ 维纳滤波(wiener filtering) 一种基于最小均方误差准则、对平稳过程的最优估计器。这种滤波器的输出与期望输出之间的均方误差为最小,因此,它是一个最佳滤波系统。它可用于提取被平稳...
  • 理想的高通滤波器(IHPF)只需将所有低频分量切掉,这些低频分量距转换原点的距离为
  • 通过对理想滤波器的加blackman窗实现滤波器设计
  • 实验作业5: 1用理想低通滤波器在频率域实现低通滤波 程序代码 I二imread'd:/jpg; I=rgb2gray(l; figure(1,imshow(l; title'原图像) s=fftshift(fft2 (I; figure(2; imshow(abs(s,[]; title'图像傅里叶变换所得频谱; ...
  • MATLAB切比雪夫带通滤波器

    万次阅读 2019-02-22 21:09:53
    原始信号由5Hz,50Hz,110Hz三种频率的正弦信号构成...​使用通带为[10,100]Hz的切比雪夫滤波器,滤波后的信号时域曲线为: 图 2 滤波信号 ​对原始信号和滤波信号作傅里叶变换,观察频谱的变化(左图为原始信号...
  • 程序中以凯瑟窗作为窗口分别设计实现了各类滤波器(低通、带通、高通、带阻、双带通等),其中使用kaiserord函数设计fir滤波器的估算参数、kaiser函数生成凯瑟窗、fir1函数生成滤波器,并绘制出了各种滤波器的频率...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,778
精华内容 1,111
关键字:

matlab理想滤波器

matlab 订阅