Bode plot of a 157 order transfer function
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Dongxu Guo
el 25 de Oct. de 2023
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 comentarios
Respuesta aceptada
Paul
el 25 de Oct. de 2023
Editada: Paul
el 26 de Oct. de 2023
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)
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 comentarios
Paul
el 26 de Oct. de 2023
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
size(myFilters)
Más respuestas (0)
Ver también
Categorías
Más información sobre Digital Filter Analysis en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!