Bode plot of a 157 order transfer function

Question

Dongxu Guo 2023-10-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2038556-bode-plot-of-a-157-order-transfer-function

编辑： Paul 2023-10-26

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2023-10-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2038556-bode-plot-of-a-157-order-transfer-function#answer_1340751

编辑：Paul 2023-10-26

在 MATLAB Online 中打开

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Dongxu Guo 2023-10-26

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 2023-10-26

在 MATLAB Online 中打开

@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.

请先登录，再进行评论。