zoukankan      html  css  js  c++  java
  • 《DSP using MATLAB》Problem 8.35

    代码:

    %% ------------------------------------------------------------------------
    %%            Output Info about this m-file
    fprintf('
    ***********************************************************
    ');
    fprintf('        <DSP using MATLAB> Problem 8.35 
    
    ');
    
    banner();
    %% ------------------------------------------------------------------------
    
    % Digital Highpass Filter Specifications:
    wphp = 0.6*pi;                 % digital passband freq in rad
    wshp = 0.4*pi;                 % digital stopband freq in rad
    Rp = 0.5;                      % passband ripple in dB
    As = 60;                       % stopband attenuation in dB
    
    Ripple = 10 ^ (-Rp/20)           % passband ripple in absolute
    Attn = 10 ^ (-As/20)             % stopband attenuation in absolute
    
    fprintf('
    *******Digital Highpass, Coefficients of DIRECT-form***********
    ');
    %[bhp, ahp] = butthpf(wphp, wshp, Rp, As)
    [bhp, ahp] = cheb1hpf(wphp, wshp, Rp, As)
    %[bhp, ahp] = cheb2hpf(wphp, wshp, Rp, As)
    %[bhp, ahp] = eliphpf(wphp, wshp, Rp, As)
    [C, B, A] = dir2cas(bhp, ahp);
    
    % Calculation of Frequency Response:
    %[dblp, maglp, phalp, grdlp, wwlp] = freqz_m(blp, alp);
    [dbhp, maghp, phahp, grdhp, wwhp] = freqz_m(bhp, ahp);
    
    % ---------------------------------------------------------------
    %    find Actual Passband Ripple and Min Stopband attenuation
    % ---------------------------------------------------------------
    delta_w = 2*pi/1000;
    Rp_hp = -(min(dbhp(ceil(wphp/delta_w+1):1:501)));                       % Actual Passband Ripple
    
    fprintf('
    Actual Passband Ripple is %.4f dB.
    ', Rp_hp);
    
    As_hp = -round(max(dbhp(1:1:ceil(wshp/delta_w)+1)));                % Min Stopband attenuation
    fprintf('
    Min Stopband attenuation is %.4f dB.
    
    ', As_hp);
    
    %% -----------------------------------------------------------------
    %%                             Plot
    %% -----------------------------------------------------------------  
    
    figure('NumberTitle', 'off', 'Name', 'Problem 8.35 Chebyshev-1 Highpass by dhpfd_bl function')
    set(gcf,'Color','white'); 
    M = 1;                          % Omega max
    
    subplot(2,2,1); plot(wwhp/pi, maghp); axis([0, M, 0, 1.2]); grid on;
    xlabel('Digital frequency in pi units'); ylabel('|H|'); title('Highpass Filter Magnitude Response');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.9441, 1]);
    
    subplot(2,2,2); plot(wwhp/pi, dbhp); axis([0, M, -100, 2]); grid on;
    xlabel('Digital frequency in pi units'); ylabel('Decibels'); title('Highpass Filter Magnitude in dB');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-70, -61, -60, -1, 0]);
    set(gca,'YTickLabelMode','manual','YTickLabel',['70'; '61'; '60';'1 ';' 0']);
    
    
    subplot(2,2,3); plot(wwhp/pi, phahp/pi); axis([0, M, -1.1, 1.1]); grid on;
    xlabel('Digital frequency in pi nuits'); ylabel('radians in pi units'); title('Highpass Filter Phase Response');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-1:1:1]);
    
    subplot(2,2,4); plot(wwhp/pi, grdhp); axis([0, M, 0, 25]); grid on;
    xlabel('Digital frequency in pi units'); ylabel('Samples'); title('Highpass Filter Group Delay');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [0:5:25]);
    
    
    
    % -----------------------------------------------------
    %              method 2
    % -----------------------------------------------------
    % Digital lowpass Filter Specifications:
    [wpLP, wsLP, alpha] = hp2lpfre(wphp, wshp);
    
    prompt = 'Please input the type of digital lp filter: 
    
      butter or cheby1 or cheby2 or ellip [butter]: ';
    type = input(prompt , 's');
    
    [bhp, ahp] = dhpfd_bl(type, wphp, wshp, Rp, As)
    
    [C, B, A] = dir2cas(bhp, ahp)
    
    
    
    % -----------------------------------------------------
    %              method 3  cheby1 function
    % -----------------------------------------------------
    % Calculation of Chebyshev-1 hp filter parameters:
    [N, wn] = cheb1ord(wphp/pi, wshp/pi, Rp, As);
    
    % Digital Chebyshev-1 Highpass Filter Design:
    [bhp, ahp] = cheby1(N, Rp, wn, 'high');
    
    [C, B, A] = dir2cas(bhp, ahp)
    
    % Calculation of Frequency Response:
    %[dblp, maglp, phalp, grdlp, wwlp] = freqz_m(blp, alp);
    [dbhp, maghp, phahp, grdhp, wwhp] = freqz_m(bhp, ahp);
    
    % ---------------------------------------------------------------
    %    find Actual Passband Ripple and Min Stopband attenuation
    % ---------------------------------------------------------------
    delta_w = 2*pi/1000;
    Rp_hp = -(min(dbhp(ceil(wphp/delta_w+1):1:501)));                       % Actual Passband Ripple
    
    fprintf('
    Actual Passband Ripple is %.4f dB.
    ', Rp_hp);
    
    As_hp = -round(max(dbhp(1:1:ceil(wshp/delta_w)+1)));                % Min Stopband attenuation
    fprintf('
    Min Stopband attenuation is %.4f dB.
    
    ', As_hp);
    
    %% -----------------------------------------------------------------
    %%                             Plot
    %% -----------------------------------------------------------------  
    
    figure('NumberTitle', 'off', 'Name', 'Problem 8.35 Chebyshev-1 Highpass by cheby1 function')
    set(gcf,'Color','white'); 
    M = 1;                          % Omega max
    
    subplot(2,2,1); plot(wwhp/pi, maghp); axis([0, M, 0, 1.2]); grid on;
    xlabel('Digital frequency in pi units'); ylabel('|H|'); title('Highpass Filter Magnitude Response');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.9441, 1]);
    
    subplot(2,2,2); plot(wwhp/pi, dbhp); axis([0, M, -100, 2]); grid on;
    xlabel('Digital frequency in pi units'); ylabel('Decibels'); title('Highpass Filter Magnitude in dB');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-70, -61, -60, -1, 0]);
    set(gca,'YTickLabelMode','manual','YTickLabel',['70'; '61'; '60';'1 ';' 0']);
    
    
    subplot(2,2,3); plot(wwhp/pi, phahp/pi); axis([0, M, -1.1, 1.1]); grid on;
    xlabel('Digital frequency in pi nuits'); ylabel('radians in pi units'); title('Highpass Filter Phase Response');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-1:1:1]);
    
    subplot(2,2,4); plot(wwhp/pi, grdhp); axis([0, M, 0, 25]); grid on;
    xlabel('Digital frequency in pi units'); ylabel('Samples'); title('Highpass Filter Group Delay');
    set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.4, 0.6, M]);
    set(gca, 'YTickMode', 'manual', 'YTick', [0:5:25]);
    

      运行结果:

            采用dhpfd_bl函数,获得Chebyshev-1型数字高通滤波器,系统函数直接形式的系数

            最小阻带衰减达到61dB,满足设计要求。

            采用dhpfd_bl函数和cheby1函数(MATLAB自带)设计数字高通滤波器,幅度谱、相位谱和群延迟对比,相差不大。

  • 相关阅读:
    【内网穿透】【natapp】web服务映射
    【javascript】日期转字符串
    【springcloud】Transaction rolled back because it has been marked as rollback-only
    MySQL 快速创建索引
    MySQL 快速导入大量数据 资料收集
    基于WinCE的JSON 类库 源码
    C# 模拟提交 Form表单的数据
    git恢复删除的分支及内容
    js数组push方法使用注意
    mint-ui的search组件如何在键盘显示搜索按钮
  • 原文地址:https://www.cnblogs.com/ky027wh-sx/p/11664291.html
Copyright © 2011-2022 走看看