精华内容
参与话题
问答
  • 下两个滤波器都是切比雪夫I型数字滤波器,不是巴特沃尔滤波器,请使用者注意!1.低通滤波器使用说明:将下列代码幅值然后以m文件保存,文件名要与函数名相同,这里函数名:lowp。function y=lowp(x,f1,f3,rp,rs,Fs) ...

    本文为转载内容,原文地址为点击打开链接

    下两个滤波器都是切比雪夫I型数字滤波器,不是巴特沃尔滤波器,请使用者注意!

    1.低通滤波器

    使用说明:将下列代码幅值然后以m文件保存,文件名要与函数名相同,这里函数名:lowp。

    function y=lowp(x,f1,f3,rp,rs,Fs)
    %低通滤波
    %使用注意事项:通带或阻带的截止频率的选取范围是不能超过采样率的一半
    %即,f1,f3的值都要小于 Fs/2
    %x:需要带通滤波的序列
    % f 1:通带截止频率
    % f 3:阻带截止频率
    %rp:边带区衰减DB数设置
    %rs:截止区衰减DB数设置
    %FS:序列x的采样频率
    % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
    % Fs=2000;%采样率
    %
    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
    --------------------------------------

    低通滤波器使用例子的代码

    fs=2000;
    t=(1:fs)/fs;
    ff1=100;
    ff2=400;
    x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t);
    figure;
    subplot(211);plot(t,x);
    subplot(212);hua_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);hua_fft(y,fs,1);%hua_fft()函数是画频谱图的函数,代码在下面给出,要保存为m文件调用
    %这段例子还调用了我自己写的专门画频谱图的函数,也给出,不然得不出我的结果;

    %画信号的幅频谱和功率谱
    %频谱使用matlab例子表示
    function hua_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
        ifnargin==3
           plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));%matlab的帮助里画FFT的方法
           %ylabel('幅值');xlabel('频率');title('信号幅值谱');
           %plot(y_f,abs(y_ft(1:nfft/2)));%论坛上画FFT的方法
        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));
        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

    下面三幅图分别是滤波前的时频图,滤波器的滤波特性曲线图和滤波后的时频图,通过图可以看出成功留下了100Hz的低频成分而把不要的高频成分去除了。

    2.高通滤波器

    function y=highp(x,f1,f3,rp,rs,Fs)
    %高通滤波
    %使用注意事项:通带或阻带的截止频率的选取范围是不能超过采样率的一半
    %即,f1,f3的值都要小于 Fs/2
    %x:需要带通滤波的序列
    % f 1:通带截止频率
    % f 2:阻带截止频率
    %rp:边带区衰减DB数设置
    %rs:截止区衰减DB数设置
    %FS:序列x的采样频率
    % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
    % Fs=2000;%采样率
    %
    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

    下面是高通滤波器的例子

    fs=2000;
    t=(1:fs)/fs;
    ff1=100;
    ff2=400;
    x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t);
    figure;
    subplot(211);plot(t,x);
    subplot(212);hua_fft(x,fs,1);
    
    %------高通测试
    z=highp(x,350,300,0.1,20,fs);
    figure;
    subplot(211);plot(t,z);
    subplot(212);hua_fft(z,fs,1);


    下面三幅图分别是滤波前的时频图,滤波器的滤波特性曲线图和滤波后的时频图,通过图可以看出成功留下了400Hz的高频成分而把不要的低频成分100Hz去除了。


    3.带通滤波器

    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的采样频率
    % f1=300;f3=500;%通带截止频率上下限
    % fsl=200;fsh=600;%阻带截止频率上下限
    % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
    % Fs=2000;%采样率
    %
    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

    带通滤波器使用例子

    %--------------
    %带通滤波器测试程序
    fs=2000;
    t=(1:fs)/fs;
    ff1=100;
    ff2=400;
    ff3=700;
    x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t)+sin(2*pi*ff3*t);
    figure;
    subplot(211);plot(t,x);
    subplot(212);hua_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);hua_fft(y,fs,1);


    %调用到的hua_fft()函数代码如下

    function hua_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
        ifnargin==3
           plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));%matlab的帮助里画FFT的方法
           %ylabel('幅值');xlabel('频率');title('信号幅值谱');
           %plot(y_f,abs(y_ft(1:nfft/2)));%论坛上画FFT的方法
        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));
        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
    运行结果如下图,第一幅是滤波前测试信号的时频图,第二幅是滤波器的滤波曲线图,第三幅是经滤波后的测试信号时频图。

    4.带阻滤波器

    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的采样频率
    % f1=300;f3=500;%通带截止频率上下限
    % fsl=200;fsh=600;%阻带截止频率上下限
    % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
    % Fs=2000;%采样率
    %
    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

    使用例子

    %带阻滤波器测试
    fs=1000;
    t=(1:fs)/fs;
    y=sin(2*pi*100*t)+sin(2*pi*150*t)+sin(2*pi*200*t);
    figure;hua_fft(y,fs,1);
    z=bands(y,110,190,140,160,0.1,30,fs);
    figure;hua_fft(z,fs,1);

    运行结果如下图,第一幅是滤波前测试信号的频谱图,第二幅是滤波器的滤波曲线图,第三幅是经滤波后的测试信号频谱图。




    展开全文
  • 带阻滤波器分两类,一类是窄带抑制带阻滤波器(简称窄带阻滤波器),另一类是宽带抑制带阻滤波器(简称宽带阻滤波器)。窄带阻滤波器一般用带通滤波器和减法器电路组合起来实现。窄带阻滤波器通常用作单一频率的陷波...
  • 带通滤波器所使用的归一化变换是将s+1/s代入低通传递函数公式得到高通滤波器的传递函数,如果将这一关系代入高通滤波器则获得一个带阻滤波器。图1表示了高通滤波器频率响应与变换后的带阻滤波器频率响应之间的等效...
  • 带阻滤波器设计指标的归一化过程和频率响应曲线参数的定义。像带通滤波器一样,带阻滤波器网络也可以由归一化低通滤波器通过适当的变换获得。  讨论了用级联低通和高通滤波器的方法设计宽带带通滤波器。用类似的...
  • 采用与全极点滤波器同样的方法,把这些电路先变换成高通滤波器,然后变换成带阻滤波器。  由于每个归一化低通滤波器可以有两种形式实现,带阻滤波器结果也有两种不同的结构,如图1所示。  图1 椭圆函数...
  • 复数低通极点变换得到一组带阻参量,其中五和∴并... 图1 带阻滤波器的vcvs实现  vcvs结构有一些不希望的特性。虽然当K>1时,调节可变的R6或R7可以改变Q值,但是由于输出的3dB带宽受传输零点影响,Q值不可能被独立
  • 原始音频信号受到一正弦噪声(幅度0.1、频率1000Hz)的污染,请使用MATLAB设计一个带阻滤波器,滤除正弦噪声信号。 根据题目的要求,思路非常清晰,使用MATLAB完成设计任务的总体思路可按照以下四步进行: 1、...
  • 二阶有源带阻滤波器

    2012-08-07 15:20:41
    二阶有源带阻滤波器multism仿真, 二阶有源带阻滤波器multism仿真。
  • 本文采用支线式结构设计了一种中等阻带带宽的带阻滤波器,使用悬置带线结构获得较高的Q值,陡峭的边带过渡特性,该滤波器还具有低成本、易于加工的优点,满足工程应用的需要。
  • Multisim 带阻滤波器

    2015-10-03 12:12:36
    求大神指导如何设计有源带阻滤波器≥﹏≤ 基于Multisim仿真的课程设计题(ಥ_ಥ)
  • 我们将u=0和v=0上式,我们可以得到如下式子。 根据上式,可以到F(0,0)的值是非常大的。这里,我们将F(0,0)称为直流分量,直流分量比其他的成分要大好几个数量级。所以,这也就是傅里叶谱为什么需要使用对数...

           1.高通滤波器

           首先,对一副图像进行如下二维傅里叶变换。

    我们将u=0和v=0带上式,我们可以得到如下式子。

    根据上式,可以到F(0,0)的值是非常大的。这里,我们将F(0,0)称为直流分量,直流分量比其他的成分要大好几个数量级。所以,这也就是傅里叶谱为什么需要使用对数变换才能看清楚的原因。
           这里,对于高通滤波器而言,由于直流分量被衰减,所以,所得到的图像的动态范围是非常狭窄的,也就造成了图像偏灰。进一步而言,保持直流(DC)分量,对别的部分进行增幅,可以增强图像的细节。这样的滤波器称为锐化滤波器。这一小节主要介绍高通滤波器与锐化滤波器。

            1.1理想高通滤波器

            
            这里的D0是滤波器的阻带半径,而D(u,v)是点到滤波器中央的距离。理想高通的滤波器的振幅特性如下所示。

    用这个滤波器对图像进行处理,可得到如下所示的结果。我们可以看到,与理想的低通滤波器一样,所得到的图像有很明显的振铃现象。结果图像从视觉上来看,有些偏暗,这是因为图像的直流分量被滤掉的原因。


           1.2巴特沃斯高通滤波器

           同样的,巴特沃斯高通滤波器也可以通过改变次数n,对过度特性进行调整。过大的n会造成振铃现象。


           1.3高斯高通滤波器

           高斯滤波器的过度特性很好,所以不会发生振铃现象。


           1.4 上述三种滤波器的实现Matlab代码

    close all;
    clear all;
    
    %% ---------Ideal Highpass Filters (Fre. Domain)------------
    f = imread('characters_test_pattern.tif');
    f = mat2gray(f,[0 255]);
    
    [M,N] = size(f);
    P = 2*M;
    Q = 2*N;
    fc = zeros(M,N);
    
    for x = 1:1:M
        for y = 1:1:N
            fc(x,y) = f(x,y) * (-1)^(x+y);
        end
    end
    
    F = fft2(fc,P,Q);
    
    H_1 = ones(P,Q);
    H_2 = ones(P,Q);
    
    for x = (-P/2):1:(P/2)-1
         for y = (-Q/2):1:(Q/2)-1
            D = (x^2 + y^2)^(0.5);
            if(D <= 60)  H_1(x+(P/2)+1,y+(Q/2)+1) = 0; end    
            if(D <= 160) H_2(x+(P/2)+1,y+(Q/2)+1) = 0; end
         end
    end
    
    G_1 = H_1 .* F;
    G_2 = H_2 .* F;
    
    g_1 = real(ifft2(G_1));
    g_1 = g_1(1:1:M,1:1:N);
    
    g_2 = real(ifft2(G_2));
    g_2 = g_2(1:1:M,1:1:N);         
    
    for x = 1:1:M
        for y = 1:1:N
            g_1(x,y) = g_1(x,y) * (-1)^(x+y);
            g_2(x,y) = g_2(x,y) * (-1)^(x+y);
        end
    end
    %% -----show-------
    figure();
    subplot(1,2,1);
    imshow(f,[0 1]);
    xlabel('a).Original Image');
    
    subplot(1,2,2);
    imshow(log(1 + abs(F)),[ ]);
    xlabel('b).Fourier spectrum of a');
    
    figure();
    subplot(1,2,1);
    imshow(H_1,[0 1]);
    xlabel('c).Ideal Highpass filter(D=60)');
    
    subplot(1,2,2);
    h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
    set(h,'EdgeColor','k');
    axis([0 P 0 Q 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    figure();
    subplot(1,2,1);
    imshow(log(1 + abs(G_1)),[ ]);
    xlabel('d).Result of filtering using c');
    
    subplot(1,2,2);
    imshow(g_1,[0 1]);
    xlabel('e).Result image');
    
    figure();
    subplot(1,2,1);
    imshow(H_2,[0 1]);
    xlabel('f).Ideal Highpass filter(D=160)');
    
    subplot(1,2,2);
    h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
    set(h,'EdgeColor','k');
    axis([0 P 0 Q 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    figure();
    subplot(1,2,1);
    imshow(log(1 + abs(G_2)),[ ]);
    xlabel('g).Result of filtering using e');
    
    subplot(1,2,2);
    imshow(g_2,[0 1]);
    xlabel('h).Result image');
    close all;
    clear all;
    
    %% ---------Butterworth Highpass Filters (Fre. Domain)------------
    f = imread('characters_test_pattern.tif');
    f = mat2gray(f,[0 255]);
    
    [M,N] = size(f);
    P = 2*M;
    Q = 2*N;
    fc = zeros(M,N);
    
    for x = 1:1:M
        for y = 1:1:N
            fc(x,y) = f(x,y) * (-1)^(x+y);
        end
    end
    
    F = fft2(fc,P,Q);
    
    H_1 = zeros(P,Q);
    H_2 = zeros(P,Q);
    
    for x = (-P/2):1:(P/2)-1
         for y = (-Q/2):1:(Q/2)-1
            D = (x^2 + y^2)^(0.5);
            D_0 = 100;
            H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^2);   
            H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^6);
         end
    end
    
    
    G_1 = H_1 .* F;
    G_2 = H_2 .* F;
    
    g_1 = real(ifft2(G_1));
    g_1 = g_1(1:1:M,1:1:N);
    
    g_2 = real(ifft2(G_2));
    g_2 = g_2(1:1:M,1:1:N);         
    
    for x = 1:1:M
        for y = 1:1:N
            g_1(x,y) = g_1(x,y) * (-1)^(x+y);
            g_2(x,y) = g_2(x,y) * (-1)^(x+y);
        end
    end
    
    %% -----show-------
    figure();
    subplot(1,2,1);
    imshow(f,[0 1]);
    xlabel('a).Original Image');
    
    subplot(1,2,2);
    imshow(log(1 + abs(F)),[ ]);
    xlabel('b).Fourier spectrum of a');
    
    figure();
    subplot(1,2,1);
    imshow(H_1,[0 1]);
    xlabel('c)Butterworth Lowpass (D_{0}=100,n=1)');
    
    subplot(1,2,2);
    h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
    set(h,'EdgeColor','k');
    axis([0 P 0 Q 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    figure();
    subplot(1,2,1);
    imshow(log(1 + abs(G_1)),[ ]);
    xlabel('d).Result of filtering using c');
    
    subplot(1,2,2);
    imshow(g_1,[0 1]);
    xlabel('e).Result image');
    
    figure();
    subplot(1,2,1);
    imshow(H_2,[0 1]);
    xlabel('f).Butterworth Lowpass (D_{0}=100,n=3)');
    
    subplot(1,2,2);
    h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
    set(h,'EdgeColor','k');
    axis([0 P 0 Q 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    figure();
    subplot(1,2,1);
    imshow(log(1 + abs(G_2)),[ ]);
    xlabel('g).Result of filtering using e');
    
    subplot(1,2,2);
    imshow(g_2,[0 1]);
    xlabel('h).Result image');
    close all;
    clear all;
    
    %% ---------Gaussian Highpass Filters (Fre. Domain)------------
    f = imread('characters_test_pattern.tif');
    f = mat2gray(f,[0 255]);
    
    [M,N] = size(f);
    P = 2*M;
    Q = 2*N;
    fc = zeros(M,N);
    
    for x = 1:1:M
        for y = 1:1:N
            fc(x,y) = f(x,y) * (-1)^(x+y);
        end
    end
    
    F = fft2(fc,P,Q);
    
    H_1 = zeros(P,Q);
    H_2 = zeros(P,Q);
    
    for x = (-P/2):1:(P/2)-1
         for y = (-Q/2):1:(Q/2)-1
            D = (x^2 + y^2)^(0.5);
            D_0 = 60;
            H_1(x+(P/2)+1,y+(Q/2)+1) = 1 - exp(-(D*D)/(2*D_0*D_0));   
            D_0 = 160;
            H_2(x+(P/2)+1,y+(Q/2)+1) = 1 - exp(-(D*D)/(2*D_0*D_0));
         end
    end
    
    G_1 = H_1 .* F;
    G_2 = H_2 .* F;
    
    g_1 = real(ifft2(G_1));
    g_1 = g_1(1:1:M,1:1:N);
    
    g_2 = real(ifft2(G_2));
    g_2 = g_2(1:1:M,1:1:N);         
    
    for x = 1:1:M
        for y = 1:1:N
            g_1(x,y) = g_1(x,y) * (-1)^(x+y);
            g_2(x,y) = g_2(x,y) * (-1)^(x+y);
        end
    end
    
    
    %% -----show-------
    figure();
    subplot(1,2,1);
    imshow(f,[0 1]);
    xlabel('a).Original Image');
    
    subplot(1,2,2);
    imshow(log(1 + abs(F)),[ ]);
    xlabel('b).Fourier spectrum of a');
    
    figure();
    subplot(1,2,1);
    imshow(H_1,[0 1]);
    xlabel('c)Gaussian Highpass (D_{0}=60)');
    
    subplot(1,2,2);
    h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
    set(h,'EdgeColor','k');
    axis([0 P 0 Q 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    figure();
    subplot(1,2,1);
    imshow(log(1 + abs(G_1)),[ ]);
    xlabel('d).Result of filtering using c');
    
    subplot(1,2,2);
    imshow(g_1,[0 1]);
    xlabel('e).Result image');
    
    figure();
    subplot(1,2,1);
    imshow(H_2,[0 1]);
    xlabel('f).Gaussian Highpass (D_{0}=160)');
    
    subplot(1,2,2);
    h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
    set(h,'EdgeColor','k');
    axis([0 P 0 Q 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    figure();
    subplot(1,2,1);
    imshow(log(1 + abs(G_2)),[ ]);
    xlabel('g).Result of filtering using e');
    
    subplot(1,2,2);
    imshow(g_2,[0 1]);
    xlabel('h).Result image');
    

           1.5 锐化滤波器

           按照之前所说的,锐化滤波器是将傅里叶谱的直流分量保留,然后将其余的成分增幅。使用锐化滤波器,可以对图像的细节进行增强,使得细节凸显出来。锐化滤波器的表达式如下所示。

          其实上式的目的很明显,就是先将原图的傅里叶保留下来,然后叠加上高通滤波器的结果,所得到的图像就是锐化后的图像了。这里为了调整锐化程度,引入了两个变量可以调整直流分量的衰减程度,可以调整高频分量的增幅程度。

           下面是代码。
    close all;
    clear all;
    clc;
    
    %% ---------The High-Fre-Emphasis Filters (Fre. Domain)------------
    f = imread('blurry_moon.tif');
    f = mat2gray(f,[0 255]);
    
    [M,N] = size(f); 
    P = 2*M;
    Q = 2*N;
    fc = zeros(M,N);
    
    for x = 1:1:M
        for y = 1:1:N
            fc(x,y) = f(x,y) * (-1)^(x+y);
        end
    end
    
    F = fft2(fc,P,Q);
    
    H_HP = zeros(P,Q);
    
    for x = (-P/2):1:(P/2)-1
         for y = (-Q/2):1:(Q/2)-1
            D = (x^2 + y^2)^(0.5);
            D_0 = 80;
            H_HP(x+(P/2)+1,y+(Q/2)+1) = 1 - exp(-(D*D)/(2*D_0*D_0));   
         end
    end
    
    G_HP = H_HP .* F;
    
    G_HFE = (1 + 1.1 * H_HP) .* F;
    
    g_1 = real(ifft2(G_HP));
    g_1 = g_1(1:1:M,1:1:N);
    
    g_2 = real(ifft2(G_HFE));
    g_2 = g_2(1:1:M,1:1:N);
    
    for x = 1:1:M
        for y = 1:1:N
            g_1(x,y) = g_1(x,y) * (-1)^(x+y);
            g_2(x,y) = g_2(x,y) * (-1)^(x+y);
        end
    end
    
    g = histeq(g_2);
    
    %g_1 = mat2gray(g_1);
    %% -----show-------
    figure();
    subplot(1,2,1);
    imshow(f,[0 1]);
    xlabel('a).Original Image');
    
    subplot(1,2,2);
    imshow(log(1 + abs(F)),[ ]);
    xlabel('b).Fourier spectrum of a');
    
    figure();
    subplot(1,2,1);
    imshow(g_1,[0 1]);
    xlabel('c).Result image of High-pass Filter');
    
    subplot(1,2,2);
    imshow(log(1 + abs(G_HP)),[ ]);
    xlabel('d).Result of filtering using e');
    
    figure();
    subplot(1,2,1);
    imshow(g_2,[0 1]);
    xlabel('e).Result image of High-Fre-Emphasis Filter');
    
    subplot(1,2,2);
    imshow(log(1 + abs(G_HFE)),[ ]);
    xlabel('f).Result of filtering using e');
    

           2.带阻滤波器

           同样的,带阻滤波器也有三种特性。高斯、巴特沃斯和理想,三种类型,其数学表达式如下所示。

           其带通滤波器可以使用上面的表格转化而得。

           带阻滤波器可以用于去除周期性噪声,为了体现带阻滤波器的特性,我们先对一幅图像增加很严重的噪声。


           在原图的傅里叶谱上添加了几个很明显的亮点。在对其做IDFT,可以看到,原图被严重的周期噪声污染了。此时,我们可以使用带阻滤波器,可以有很好的去噪效果。为了避免振铃现象,选择使用如下所示巴特沃斯带阻滤波器,所用滤波器的次数为2次。使用空间域的操作,要去除这种噪声基本是不可能的,这也是频域内的操作的优点。

      
    close all;
    clear all;
    clc;
    %% ---------------------Add Noise-------------------------
    f = imread('left.tif');
    f = mat2gray(f,[0 255]);
    
    [M,N] = size(f);
    P = 2*M;
    Q = 2*N;
    fc = zeros(M,N);
    
    for x = 1:1:M
        for y = 1:1:N
            fc(x,y) = f(x,y) * (-1)^(x+y);
        end
    end
    
    F = fft2(fc,P,Q);
    
    H_NP = ones(P,Q);
    
    for x = (-P/2):1:(P/2)-1
         for y = (-Q/2):1:(Q/2)-1
            D = 2;
            
            v_k = -200; u_k = 150;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            
            v_k = 200; u_k = 150;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            
            v_k = 0; u_k = 250;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            
            v_k = 250; u_k = 0;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            
            
            H_NP(x+(P/2)+1,y+(Q/2)+1) = 1 + 700*(1 - H_NP(x+(P/2)+1,y+(Q/2)+1));
         end
    end
    
    G_Noise = F .* H_NP;
    
    g_noise = real(ifft2(G_Noise));
    g_noise = g_noise(1:1:M,1:1:N);     
    
    for x = 1:1:M
        for y = 1:1:N
            g_noise(x,y) = g_noise(x,y) * (-1)^(x+y);
        end
    end
    
    
    %% ---------Bondpass Filters (Fre. Domain)------------
    H_1 = ones(P,Q);
    
    for x = (-P/2):1:(P/2)-1
         for y = (-Q/2):1:(Q/2)-1
            D = (x^2 + y^2)^(0.5);
            D_0 = 250;
            W = 30;
            H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+((D*W)/((D*D) - (D_0*D_0)))^6);   
         end
    end
    
    G_1 = H_1 .* G_Noise;
    
    g_1 = real(ifft2(G_1));
    g_1 = g_1(1:1:M,1:1:N);     
    
    for x = 1:1:M
        for y = 1:1:N
            g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        end
    end
    
    %% -----show-------
    close all;
    
    figure();
    subplot(1,2,1);
    imshow(f,[0 1]);
    xlabel('a).Original Image');
    
    subplot(1,2,2);
    imshow(log(1 + abs(F)),[ ]);
    xlabel('b).Fourier spectrum of a');
    
    
    figure();
    subplot(1,2,2);
    imshow(log(1 + abs(G_Noise)),[ ]);
    xlabel('c).Fourier spectrum of b');
    
    subplot(1,2,1);
    imshow(g_noise,[0 1]);
    xlabel('b).Result of add noise');
    
    
    figure();
    subplot(1,2,1);
    imshow(H_1,[0 1]);
    xlabel('d).Butterworth Bandpass filter(D=355 W=40 n =2)');
    
    subplot(1,2,2);
    h = mesh(1:20:Q,1:20:P,H_1(1:20:P,1:20:Q));
    set(h,'EdgeColor','k');
    axis([0 Q 0 P 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    
    figure();
    subplot(1,2,2);
    imshow(log(1 + abs(G_1)),[ ]);
    xlabel('e).Fourier spectrum of f');
    
    subplot(1,2,1);
    imshow(g_1,[0 1]);
    xlabel('f).Result of denoise');
    

           3.陷波滤波器(Notch Filter)

           陷波滤波器也用于去除周期噪声,虽然带阻滤波器也能可以去除周期噪声,但是带阻滤波器对噪声以外的成分也有衰减。而陷波滤波器主要对,某个点进行衰减,对其余的成分不损失。使用下图将更容易理解。

           从空间域内看,图像存在着周期性噪声。其傅里叶频谱中,也能看到噪声的所在之处(这里我用红圈标注出来了,他们不是数据的一部分)。我们如果使用带阻滤波器的话,是非常麻烦的,也会使得图像有损失。陷波滤波器,能够直接对噪声处进行衰减,可以有很好的去噪效果。
           其表达式如下所示,陷波滤波器可以通过对高通滤波器的中心,进行位移就可以得到了。

    这里,由于傅里叶的周期性,傅里叶频谱上不可能单独存在一个点的噪声,必定是关于远点对称的一个噪声对。这里的需要去除的噪声点,其求取的方式如下所示。
           针对于上图,我们设计如下滤波器,去进行去噪。
    (图片下标错了,后续更新改过来!)
           所得到的结果,如下所示。噪声已经被去除了,画质得到了很大的改善。

                                                                  (图片下标错了,后续更新改过来!)
    close all;
    clear all;
    clc;
    
    %% ---------Butterworth Notch filter (Fre. Domain)------------
    f = imread('car_75DPI_Moire.tif');
    f = mat2gray(f,[0 255]);
    
    [M,N] = size(f);
    P = 2*M;
    Q = 2*N;
    fc = zeros(M,N);
    
    for x = 1:1:M
        for y = 1:1:N
            fc(x,y) = f(x,y) * (-1)^(x+y);
        end
    end
    
    F = fft2(fc,P,Q);
    
    H_NF = ones(P,Q);
    
    for x = (-P/2):1:(P/2)-1
         for y = (-Q/2):1:(Q/2)-1
            D = 30;
            
            v_k = 59; u_k = 77;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            
            v_k = 59; u_k = 159;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            
            v_k = -54; u_k = 84;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            
            v_k = -54; u_k = 167;
            D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
            D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5);
            H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);
         end
    end
    
    G_1 = H_NF .* F;
    
    g_1 = real(ifft2(G_1));
    g_1 = g_1(1:1:M,1:1:N);     
    
    for x = 1:1:M
        for y = 1:1:N
            g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        end
    end
    
    %% -----show-------
    close all;
    
    figure();
    subplot(1,2,1);
    imshow(f,[0 1]);
    xlabel('a).Original Image');
    
    subplot(1,2,2);
    imshow(log(1 + abs(F)),[ ]);
    xlabel('b).Fourier spectrum of a');
    
    figure();
    subplot(1,2,1);
    imshow(H_NF,[0 1]);
    xlabel('c).Butterworth Notch filter(D=30 n=2)');
    
    subplot(1,2,2);
    h = mesh(1:10:Q,1:10:P,H_NF(1:10:P,1:10:Q));
    set(h,'EdgeColor','k');
    axis([0 Q 0 P 0 1]);
    xlabel('u');ylabel('v');
    zlabel('|H(u,v)|');
    
    figure();
    subplot(1,2,2);
    imshow(log(1 + abs(G_1)),[ ]);
    xlabel('e).Fourier spectrum of d');
    
    subplot(1,2,1);
    imshow(g_1,[0 1]);
    xlabel('d).Result image');
    
    原文发于博客:http://blog.csdn.net/thnh169/ 

    展开全文
  • 然后,低通和高通滤波器单独设计,并使其输人端并联,把两个输出相加,就形成了带阻滤波器。  对于全极点滤波器,宽带设计方法适用于截止频率之间的间隔在一倍频程以上的情况,这样在输出端相加时阻带内的相互影响...
  • 先给出一张图像,来说明带阻滤波器
    先给出一张图像,来说明带阻滤波器和带通滤波器的关系。
          
    由上图可以看出这张图片的中低成份占的比较多,而高频成份比较少。因为白色的像素都集中在中点和离中心的附近。对于频谱图,由白色代表某一频率点有响应,也就是原图含有该频率的成份。可以举个例子,假如我对一张图片加入椒盐噪声,这个椒盐噪声属于全频谱的。可以得到下图示
      
        很明显图像被污染了,频谱图出现了很多白色点,就代表这张图片在低频成份很多的情况下,又额外的增加了高频中频低频成份,如上第三图就出现了凹凸不平山峰一样。
        解释完频率后来看带通滤波器和带阻滤波器。

    再分别给出它们傅里叶反变换图

    图4为带阻滤波器的三维视图,图5只含中频信号的图像,图6只含高频和低频信号。读者可以自己举例子来进行高频和低频的分析,一样的原理。

    只给出带阻滤波器的相关代码
    I=imread('C:\Users\hlx\Desktop\1.jpg');   %读入原图像文件
    I2=rgb2gray(I);
    I1=double(I2);
    fftI=fft2(I1);         % 二维离散傅立叶变换
    sfftI=fftshift(fftI);    % 直流分量移到频谱中心
    [N1,N2]=size(sfftI);
    n=2;
    d0=10;
    d1=200;
    n1=floor(N1/2);
    n2=floor(N2/2);
    for i=1:N1
       for j=1:N2
        d=sqrt((i-n1)^2+(j-n2)^2);
        if d<=d0 || d>=d1
                 h=0;
        else
            h=1;
        end
        result(i,j)=h*sfftI(i,j);
       end
    end
    sfftI=result;
    RR=real(sfftI);      % 取傅立叶变换的实部
    II=imag(sfftI);      % 取傅立叶变换的虚部
    A=sqrt(RR.^2+II.^2); % 计算频谱幅值
    A=(A-min(min(A)))/(max(max(A))-min(min(A)))*255;    %归一化
    imshow(A);





    展开全文
  • 本文所提出的带阻滤波器单元大大改善了带阻滤波器的阻带衰减,扩展了带宽,具有良好的带外响应。由于时间有限没有对该结构的进行深层研究,我们在下阶段将对该结构产生高阻带衰减以及为何能将带宽展宽进行更深入的...
  • 设计FIR带阻滤波器,有源程序,欢迎下载
  • 这是我课程设计的一部分,自己编写的程序,呵呵.
  • 理想带阻滤波器

    千次阅读 2016-07-01 18:44:37
    理想带阻滤波器是频域滤波的一种,主要是抑制距离频域中心D0、一个圆环区域的频域成分,因此可以使用理想带阻滤波器来消除频率分布在圆环上的周期噪声。 公式如下   产生的3D理想带阻滤波器的图形如下    ...

         

         理想带阻滤波器是频域滤波的一种,主要是抑制距离频域中心D0、一个圆环区域的频域成分,因此可以使用理想带阻滤波器来消除频率分布在圆环上的周期噪声。

    公式如下


                      


    产生的3D理想带阻滤波器的图形如下


                            

     

     

    现在看一下对一个图像加点白噪声,然后滤波一下,下面是原图和滤波后的图


                            

    展开全文
  • 我们将u=0和v=0上式,我们可以得到如下式子。 根据上式,可以到F(0,0)的值是非常大的。这里,我们将 F(0,0)称为直流分量,直流分量比其他的成分要大好几个数量级。所以,这也就是傅里叶谱为什么需要使用...
  • 要求 将以下带阻滤波器的设计指标归一化。  带阻滤波器,在200Hz和800Hz处下降3dB,在300Hz和5O0Hz处最小衰减为40dB。  解 ①确定上限截止频率与下限截止频率之比:  所以,此滤波器是带宽型。  ②把设计...
  • 低通:低频通过 高通:高频通过 带通:一定频率范围通过 带阻:阻止一定频率范围的信号  带通滤波器工作原理:  一个理想的带通滤波器应该有一个完全平坦的... 带阻滤波器工作原理:  能通过大多数频率分量...
  • 要求 按例1的指标设计一个有源带阻滤波器,其增益为+6dB。  解 ①例1的带阻变换可得表1所示的3节滤波器的技术指标。  表1 3阶滤波器的技术指标  图1 fr=f∞的带阻电路  欢迎转载,信息来自维库电子...
  • 这个文档涉及了带通和带阻滤波器实现过程中用到的基础知识,实现公式,其中有几个实例,对于初学者能很快上手
  • matlab带阻滤波器设计

    千次阅读 2019-12-19 11:50:33
    MATLAB带阻滤波器设计## ** %程序设计 %任务书中给出的要求为中心频率200Hz,带宽150Hz。 %故设上通带截止频率为110Hz,下通带截止频率290Hz,阻带上限频率140Hz,阻带下限频率260Hz。 %此处仅以Boxcar窗为示例,...
  • matlab设计模拟带阻滤波器

    万次阅读 热门讨论 2018-10-24 22:13:11
    简单记录下在matlab上如何设计出模拟的带阻滤波器,包括:巴特沃斯滤波器、切比雪夫I型滤波器、切比雪夫II型滤波器、椭圆型滤波器。 %设计带阻滤波器 %巴特沃斯、切比雪夫I型、切比雪夫II型、椭圆型滤波器 clear ...
  • 基于共面波导的可调带阻滤波器设计,舒昶,邓中亮,本文设计了一种基于共面波导结构的具有MEMS开关的可调带阻滤波器,该可调带阻滤波器通过改变电压调节悬空梁的高度,使得中心频率��
  • 1、低通:(Low-pass filter)是容许低于截止频率的信号通过,但高于截止频率的信号不能通过的电子...3、带通:是指能通过某一频率范围内的频率分量、但将其他范围的频率分量衰减到极低水平的滤波器,与带阻滤波器的...
  • FIR高通/低通/带通/带阻滤波器设计,包括M文件和Simulink两种方法的原码
  • 数电数电实验五 RC有源低通与带阻滤波器

空空如也

1 2 3 4 5 ... 20
收藏数 23,410
精华内容 9,364
关键字:

带阻滤波器