Bode plot of a 157 order transfer function

I am trying to evaluate a cascade of 480 stages of second-order low-pass filters (to mimic human cochleas) in s-domain. The resulted filter is the 157 order (has 158 non-zero coeffiecients in denominator). The problem is that I got nothing from functions bode(myfilter) or lsim(myfilter, input_signal, t). But I can plot some kind of frequency response with freqresp(). Is it because bode() and lsim() have order limitations? I also attached the code below:

% A second order low pass filter
% Cascade of 480 stages

Q = 0.79;
numerator = 1;

mResult = 1; % Accumulated transfer function


for k = 1:480
    tau(k) = 1.014^(k-1)/(40000*pi);
    myFilters(k) = tf(1, [tau(k)^2, 1/Q*tau(k), 1]); 
end

% Place selection
for t = 1:480
    mResult = mResult*myFilters(t);
end

f = 1:1:30000; % 1-30kHz
w = 2*pi*f;

filterUndertest = mResult;


out = freqresp(filterUndertest, w);
for k = 1:30000
    AbsOut (k) = abs(out(:, :, k));
end

semilogx(f,AbsOut);  % This thing works

bode(mResult);       % But this not

Answer 1

[Paul]

Hi Dongxu Guo,

That's quite a filter!

Using a tf for a for such a high order filter is not likely to work due to large rounding errors. BTW, isn't the filter order 480*2 = 960?

For Bode plots, my quick experiments indicate the frd representation seems to be the best approach, though I have no idea if the final result looks like it should.

% A second order low pass filter
% Cascade of 480 stages
Q = 0.79;
numerator = 1;

f = 1:1:30000; % 1-30kHz
f = logspace(0,5,1000); % use less points for run time 
f(f>30000) = [];
w = 2*pi*f;

mResult = frd(tf(1),w); % Accumulated transfer function

N = 480;
for k = 1:N %1:480
    tau(k) = 1.014^(k-1)/(40000*pi);
    myFilters{k} = tf(1, [tau(k)^2, 1/Q*tau(k), 1]);
end

% Place selection
for t = 1:N  %480
    mResult = mResult*frd(myFilters{t},w);
end

%{
f = 1:1:30000; % 1-30kHz
w = 2*pi*f;
filterUndertest = mResult;
out = freqresp(filterUndertest, w);
for k = 1:30000
    AbsOut (k) = abs(out(:, :, k));
end
%}

[m,p,w]=bode(mResult,w);
bode(mResult)

The Bode plot doesn't go all the way to f(end) because eventually the magnitudes are numerically 0 (I assume the phase plot stops to be consistent with the magnitude plot, though the calculated phase wil likely be inncacurate anyway).

m(end),p(end)

ans = 0

ans = -41040

figure
semilogx(f,abs(m(:)));

For lsim, you might want to try using a loop to lsim each one with the output of the previous stage being the input of the current stage in the loop.

[Sam Chak]

The 2nd-order low-pass filter stages are cascaded 480 times in the for-loop, making it a massively ultra-high-order filter. The order of the filter can be checked using order() command.

NS = order(myFilters)

NS = 960

[Dongxu Guo]

Thanks! It looks nice to me, and I'll try looping lsim to get a output in time domain. Sorry for mistaking the order - it should be 960.

[Paul]

@Dongxu Guo

Here is the code from the question

Q = 0.79;
numerator = 1;
mResult = 1; % Accumulated transfer function
for k = 1:480
    tau(k) = 1.014^(k-1)/(40000*pi);
    myFilters(k) = tf(1, [tau(k)^2, 1/Q*tau(k), 1]); 
end

I just want to point out that myFilters constructed this way is not really an arrary of transfer functions, but rather a transfer function with 1 output and 480 inputs. Assuming we really want to just store separate transfer functions, consider using a cell array (thought the current code seems to work fine), which I've now updated in the Answer.

whos myFilters

  Name           Size              Bytes  Class    Attributes

  myFilters      1x480            131585  tf                 

size(myFilters)

Transfer function with 1 outputs and 480 inputs.