chebyshev II BPF magnitude response is not equal to 1 in the passband - incorrect magnitude response when using freqz(b,a)

I need to exctract frequency bands of eeg signal using iir filterbanks , i got the optimum order of chebyshev type II BPF , I designed the filter but I get magntitude response not equal to 1 in the passbands , also I get incorrect magntuide response when using freqz(b,a) , code given below and magnitude response figure when using freqz(sos,[],Fs)and freqz(b,a) , kindly check .

    clear 
    close all
    clc
    % Delta (0.5 - 4) Hz
    % theta (4 - 7.9) Hz
    % Lower Alpha (7.9 - 10) Hz
    % Upper Alpha (10.1 - 12.9) Hz
    % LOW BETA (12-15) Hz
    % MID BETA (15-18) Hz
    Fs=173.61;% Sampling Frequency
    Fn=Fs/2; % Nyquist Frequency
    Wpm=[0.5 4;4 7.9;7.9 10;10.1 12.9;12 15;15 18];% Passbands
    Rp=3;% Passband Ripple
    Rs=40;% Stopband Ripple
    for m=1:size(Wpm,1)
        Wp=Wpm(m,:)/Fn;
        Ws=Wp.*[0.9 1.05];
        %% Design a bandpass filters
        [n(m,:),Ws]=cheb2ord(Wp,Ws,Rp,Rs); % Determine filter Order
        [b{m},a{m}]=cheby2(ceil(n(m,:)/2),Rs,Ws,"bandpass");
        [z,p,k]=cheby2(ceil(n(m,:)/2),Rs,Ws,"bandpass"); % Transfer Function Coefficients
        [sos{m},g{m}]=zp2sos(z,p,k);
        [H{m},f]=freqz(sos{m},[],Fs);
        [H1{m},f]=freqz(b{m},a{m},[],Fs);
        subplot(size(Wpm,1),2,2*m-1)
        plot(f,abs(H{m}),"LineWidth",1.5)
        title(sprintf('Frequency Response of digital Chebyshev II order n = %d',n(m,:)))
        ylim([0 1.2*max(abs(H{m}))]),xlabel('Hz') , ylabel('Magnitude')
      grid on , xticks([Wp*Fn Fn])
        subplot(size(Wpm,1),2,2*m)
        plot(f,abs(H1{m}),"LineWidth",1.5)
        title(sprintf('Frequency Response of digital Chebyshev II order n = %d',n(m,:)))
        ylim([0 1.2*max(abs(H1{m}))]),xlabel('Hz'), ylabel('Magnitude'),
      grid on ,    xticks([Wp*Fn Fn])
        grid on , xticks([Wp*Fn Fn])
    end

Answer 1

[Star Strider]

In general, using the polynomial numerator and denominator implementation of digital filters is not optimal. Using the second-order-section representation is preferred, largely for stability reasons (although stability might not be an issue with your filters).

If you have a MATLAB version that supports the bandpass function, I suggest using it instead. It will design a very efficient elliptic filter for each set of frequencies, and present it in the second output of the function. You can use that digital-filter object directly with filtfilt. You can also design elliptic filters similarly to the way you designed your Chebychev filters (using command-line functions), however the bandpass function makes the process much easier.

I looked at your code and ran it. I agree with oyur interpretation of the results. Using elliptic filters instead might be preferable. This has been my experience with signal processing of physiologic signals.

[Mohamad]

Ok i will try elliptic filters design , thanks a lot

[Star Strider]

My pleasure!

To force bandpass and related functions to produce an elliptic filter, include:

ImpulseResponse='iir'

as an argument in the original call.

.

Answer 2

[Paul]

Hi Mohamad,

When using the two-output form of zp2sos the gain, g, must be applied to the output of freqz to get the correct frequency response. See code below where g is included in the plot and ylim commands for the sos case. Alternatively, use the one-output form of zp2sos.

I wasn't sure in the original code why the filter order from cheb2ord was halved, so I also removed that from the code below.

The sos freqz results are now as expected for all cases. The freqz results from the transfer function form are not good for the first two cases, as can sometimes happen (documentation discussion). Stick with the sos or ctf forms (the latter introduced in 2024b).

Fs=173.61;% Sampling Frequency
Fn=Fs/2; % Nyquist Frequency
Wpm=[0.5 4;4 7.9;7.9 10;10.1 12.9;12 15;15 18];% Passbands
Rp=3;% Passband Ripple
Rs=40;% Stopband Ripple
for m = 1:size(Wpm,1)
    Wp=Wpm(m,:)/Fn;
    Ws=Wp.*[0.9 1.05];
    %% Design a bandpass filters
    [n(m,:),Ws]=cheb2ord(Wp,Ws,Rp,Rs); % Determine filter Order
    %[b{m},a{m}]=cheby2(ceil(n(m,:)/2),Rs,Ws,"bandpass");
    [b{m},a{m}]=cheby2(n(m),Rs,Ws,"bandpass");
    %[z,p,k]=cheby2(ceil(n(m,:)/2),Rs,Ws,"bandpass"); % Transfer Function Coefficients
    [z,p,k]=cheby2(n(m),Rs,Ws,"bandpass"); % Transfer Function Coefficients
    [sos{m},g{m}]=zp2sos(z,p,k);
    [H{m},f]=freqz(sos{m},[],Fs);
    [H1{m},f]=freqz(b{m},a{m},[],Fs);
    subplot(size(Wpm,1),2,2*m-1)
    
    plot(f,g{m}*abs(H{m}),"LineWidth",1.5)  % correction here!
    
    title(sprintf('Frequency Response of digital Chebyshev II order n = %d',n(m,:)))
    
    % changed scaling here
    ylim([0 1.2*max(abs(g{m}*H{m}))]),xlabel('Hz') , ylabel('Magnitude'), grid on, xticks([Wp*Fn Fn])
    
    subplot(size(Wpm,1),2,2*m)
    plot(f,abs(H1{m}),"LineWidth",1.5)
    title(sprintf('Frequency Response of digital Chebyshev II order n = %d',n(m,:)))
    ylim([0 1.2*max(abs(H1{m}))]),xlabel('Hz'), ylabel('Magnitude'), grid on, xticks([Wp*Fn Fn])
end

[Mohamad]

Dear Paul, thanks a lot for explaining how to get the correct freqz results . I divide the filter order by 2 , because documentation mentioned for BPF & BSF the resulting order is 2n so i divide by 2 to get the optimum filter order computed by cheby2ord .

[Paul]

The documentation of cheb2ord shows an example where the n output from cheb2ord is used as the input to cheby2 (not n/2) to design a filter that meets the Rp and Rs specifications. In other words, cheb2ord and cheby2 are "coordinated" as to the meaning of n for bandpass and bandstop filters.

[Mohamad]

[b,a] = cheby2(3,50,[0.2 0.6],'stop'); freqz(b,a) But here it returns 6th order BSF , as shown in the documentation .

[Paul]

Yes, the order of the banstop filter is 2*n. And cheb2ord returns a value of n such that the filter returned from cheby2 with n as the input will meet the specification inputs to cheb2ord. Check the doc link I provided in the previous comment.

[Mohamad]

Ok , Paul , now I understand in order to meet the passband & stopband ripples , I must use the computed n from cheby2ord , then I have to enter it to cheby2 , thanks and best regards