unexpected discontinuity in graphic

Hello everyone,

I am writing a program to compare two Chebyshev Type 1 filters of different orders. On the 16th order of the filter, there is an unexpected discontinuity in the graphic. Does anyone know why this happens?

Here is my code:

function LAB3_ASK2
close all
clearvars

% Example usage:
order1 = 2;
order2 = 16;
Pripple = 3;
Ts = 0.2;

% Compare filters and plot responses
chebyshevType1Comparison(order1, order2, Pripple, Ts);
end

function chebyshevType1Comparison(order1, order2, Pripple, Ts)
    % Design and compare Chebyshev Type I filters of different orders

    wc = 2;         % Cutoff angular frequency
    fc = wc/(2*pi); % Cutoff natural frequency
    fs = 1/Ts;      % Sample frequency
    N = 256;        % Number of samples

    % Order 1 (2nd order)
    [num,den] = cheby1(order1, Pripple, 2*fc/fs, 'high');
    [H2, f2] = freqz(num, den, N, fs);

    % Order 2 (16th order)
    [num,den] = cheby1(order2, Pripple, 2*fc/fs, 'high');
    [H16, f16] = freqz(num, den, N, fs);

    % Plotting
    plotFilterResponse(2*f2/fs, H2, 2*f16/fs, H16, order1, order2);
end

function plotFilterResponse(f2, H2, f16, H16, order1, order2)
    % Plot the frequency response of Chebyshev Type I filters

    % Plot individual responses
    figure;
    plot(f2, 20*log10(abs(H2)), 'b');
    hold on;
    plot(f16, 20*log10(abs(H16)), 'r');
    grid on;
    title('Frequency response');
    ylabel('Magnitude (dB)');
    xlabel('Normalized frequency (π * radians/sample)');
    legend([num2str(order1) 'nd order'], [num2str(order2) 'th order'], 'Location', 'southwest');

    % Save the plot
    saveas(gcf, ['Lab3/Chebyshev_Type1_Comparison_' num2str(order1) '_' num2str(order2) '.png']);
end

Answer 1

[Star Strider]

My guess is that the filter is unstable. This occurs relatively frequently with transfer function realiations of filters.

I would instead use zero-pole-gain output and use second-order-section implementation:

[z,p,k] = cheby1(order1, Pripple, 2*fc/fs, 'high');

[sos1,g1] = zp2sos(z,p,k);

[H2, f2] = freqz(sos1, N, fs);

That should be stable.

Do the same for the second filter.

LAB3_ASK2

function LAB3_ASK2
close all
clearvars

% Example usage:
order1 = 2;
order2 = 16;
Pripple = 3;
Ts = 0.2;

% Compare filters and plot responses
chebyshevType1Comparison(order1, order2, Pripple, Ts);
end

function chebyshevType1Comparison(order1, order2, Pripple, Ts)
    % Design and compare Chebyshev Type I filters of different orders

    wc = 2;         % Cutoff angular frequency
    fc = wc/(2*pi); % Cutoff natural frequency
    fs = 1/Ts;      % Sample frequency
    N = 256;        % Number of samples

    % Order 1 (2nd order)
    [z,p,k] = cheby1(order1, Pripple, 2*fc/fs, 'high');
    [sos1,g1] = zp2sos(z,p,k);
    
    [H2, f2] = freqz(sos1, N, fs);

    % Order 2 (16th order)
    [z,p,k] = cheby1(order2, Pripple, 2*fc/fs, 'high');
    [sos2,g2] = zp2sos(z,p,k);
    
    [H16, f16] = freqz(sos2, N, fs);

    % Plotting
    plotFilterResponse(2*f2/fs, H2, 2*f16/fs, H16, order1, order2);
end

function plotFilterResponse(f2, H2, f16, H16, order1, order2)
    % Plot the frequency response of Chebyshev Type I filters

    % Plot individual responses
    figure;
    plot(f2, 20*log10(abs(H2)), 'b');
    hold on;
    plot(f16, 20*log10(abs(H16)), 'r');
    grid on;
    title('Frequency response');
    ylabel('Magnitude (dB)');
    xlabel('Normalized frequency (π * radians/sample)');
    legend([num2str(order1) 'nd order'], [num2str(order2) 'th order'], 'Location', 'southwest');

    % Save the plot
    % saveas(gcf, ['Lab3/Chebyshev_Type1_Comparison_' num2str(order1) '_' num2str(order2) '.png']);
end

EDIT — (22 Nov 2023 at 14:28)

Tested code with my suggestions.

.

[Paul]

I don't think the plot is correct for sos1. freqz requires at least two sections when using the sos form. If the sos input has one row it is interpreted as a numerator in the b,a syntax.

[Star Strider]

The 16-order filter was the issue. That got solved.

Answer 2

[Paul]

Hi Konstantinos,

Your plot, the plot I got on my local machine, and the plot generated here aren't quite the same,presumbably due to machine and/or library differences in computing very small numbers. But the discontinuity you see looks like a result of H16 evaluating to exactly 0 at a few points. Here, they are the 1st, 2nd, and 6th points, and also the 1st point of H2.

LAB3_ASK2

Indices where H2 is identically zero

ans = 1

Indices where H16 is identically zero

ans = 3×1

1 2 6

xlim([0 0.2])

When 0 is converted to dB we get

20*log10(0)

ans = -Inf

-Inf values aren't included in line plots, hence the plot looks choppy.

Having said that, the tf form for a high order filter is not recommended, as discussed at cheby1 (Limitations section).

function LAB3_ASK2
close all
clearvars
% Example usage:
order1 = 2;
order2 = 16;
Pripple = 3;
Ts = 0.2;
% Compare filters and plot responses
chebyshevType1Comparison(order1, order2, Pripple, Ts);
end
function chebyshevType1Comparison(order1, order2, Pripple, Ts)
    % Design and compare Chebyshev Type I filters of different orders
    wc = 2;         % Cutoff angular frequency
    fc = wc/(2*pi); % Cutoff natural frequency
    fs = 1/Ts;      % Sample frequency
    N = 256;        % Number of samples
    % Order 1 (2nd order)
    [num,den] = cheby1(order1, Pripple, 2*fc/fs, 'high');
    [H2, f2] = freqz(num, den, N, fs);
    disp('Indices where H2 is identically zero')
    find(H2 == 0)

    % Order 2 (16th order)
    [num,den] = cheby1(order2, Pripple, 2*fc/fs, 'high');
    [H16, f16] = freqz(num, den, N, fs);
    disp('Indices where H16 is identically zero')
    find(H16 == 0)
    % Plotting
    plotFilterResponse(2*f2/fs, H2, 2*f16/fs, H16, order1, order2);
end
function plotFilterResponse(f2, H2, f16, H16, order1, order2)
    % Plot the frequency response of Chebyshev Type I filters
    % Plot individual responses
    figure;
    plot(f2, 20*log10(abs(H2)), 'b');
    hold on;
    plot(f16, 20*log10(abs(H16)), 'r');
    grid on;
    title('Frequency response');
    ylabel('Magnitude (dB)');
    xlabel('Normalized frequency (π * radians/sample)');
    legend([num2str(order1) 'nd order'], [num2str(order2) 'th order'], 'Location', 'southwest');
    % Save the plot
%    saveas(gcf, ['Lab3/Chebyshev_Type1_Comparison_' num2str(order1) '_' num2str(order2) '.png']);
end