RC circuit and its frequency response

I am trying to model a RC circuit's response to the complex exponential input. But the resulting graph is not really good. I have tried to change the modeling period, the transfer function, and the expression for exponential gain. But it only get worse. Here is my code. So basically, I am trying to build the Bode plot manually.

%% 2
f_s = 44.1*10^3; % sampling frequency
tau = 0.005; % tau
dt = 1/f_s; % delta time

i = 1;
M_1 = zeros(1,9); % the magnitude of low pass
M_2 = zeros(1,9); % magnitude of high pass
P_1 = zeros(1,9); % the phase of low pass
P_2 = zeros(1,9); % phase of high pass
H1 = zeros(1,9); % initializing the complex gain for low pass
H2 = zeros(1,9); % initializing the complex gain for high pass

for f = logspace(1,4,48)

    w = 2*pi*f; % angular frequency
    T = 1/f;
    tau = 0.005;

    t_end = max(10*T,8*tau); % find the steady state function
    t = 0:dt:t_end;

    vin = exp(1j*w*t); % This is input as complex exponential;

    a = [1, 1/tau]; % low pass filter
    b = [1/tau];
    sys_low = tf(b,a);

    a = [1, 1/tau]; % high pass filter
    b = [1, 0];
    sys_high = tf(b,a);

    [y1,t_out_1] = lsim(sys_low,vin,t); % low pass for f using lsim
    [y2,t_out_2] = lsim(sys_high,vin,t); % high pass for f using lsim

    y1 = y1';
    y2 = y2'; % make same dimension

    idx = round(length(t)*0.8); % use the last 20% of the data

    t_segment = t(idx:end)'; % make the segment of the last period
    y1_segment = y1(idx:end)';
    y2_segment = y2(idx:end)';
    x_segment = vin(idx:end)';

    H1(i) = (y1_segment'*x_segment)/(x_segment'*x_segment); % the complex gain
    H2(i) = (y2_segment'*x_segment)/(x_segment'*x_segment);

    i = i+1;

end
f = logspace(1,4,48);

M_1 = 20*log10(abs(H1)); % the magnitude of steady state output of low pass
M_2 = 20*log10(abs(H2)); % the magnitude of steady state output of high pass

P_1 = (angle(H1)/pi); % the phase of the steady state output of low pass
P_2 = (angle(H2)/pi); % the phase of the steady state output of high pass

figure;
subplot(2,2,1);
semilogx(f,M_1); % plot the magnitude over frequency
% xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Figure 1');
grid on;

subplot(2,2,2);
semilogx(f,M_2); % plot the magnitude over frequency for high pass filter
% xlabel('Frequency (Hz)');
% ylabel('Magnitude (dB)');
title('Figure 2');
grid on;

subplot(2,2,3);
semilogx(f,P_1); % plot the phase over frequency for low pass filter
xlabel('Frequency (Hz)');
ylabel('Phase (π radians)');
title('Figure 3');
grid on;

subplot(2,2,4);
semilogx(f,P_2); % plot the phase over frequency for high pass filter
xlabel('Frequency (Hz)');
% ylabel('Phase (π radians)');
title('Figure 4');
grid on;

% sgtitle('Phase and Magnitude Response of Steady State for Low and High Pass Filters');
sgtitle('Bode Plots');

●

% Fig. 1: Magnitude of Steady State Response Over Frequencies for Low Pass Filter
% Fig. 2: Magnitude of Steady State Response Over Frequencies for High Pass Filter
% Fig. 3: Phase of Steady State Response Over Frequencies for Low Pass Filter
% Fig. 4: Phase of Steady State Response Over Frequencies for High Pass Filter

[Mathieu NOE]

I am not familiar with this way of computing a complex gain (directly from complex time signals ??)

H1(i) = (y1_segment'*x_segment)/(x_segment'*x_segment); % the complex gain

smells similar to a cross spectrum divided by the input autospectrum but without any fft operation

who invented this ? sources ?

I would have done with a more traditional way (send a real sine wave and project the output on cos / sin carriers (basically a DFT ))

[Paul]

Seems like the gain of the transfer function would just be the steady-state amplitude of the complex output (because the input is unit amplitude) and the phase of the transfer function would be the steady steady state phase of the ratio of the output to the input.

Answer 1

[Mathieu NOE]

still don't know if your method is valid , so I suggest mine (simply uses the same equations as a classical DFT)

it's not 100% accurate as a correct dft would require the computaion of the real and imaginary part of H with sin/cos carrier signals starting and ending at zero crossing points. this is not easy to do with your sine sweep and fixed sampling rate

a better approach would be to have a sampling rate which is proportionnal to your sine wave frequency (integer ratio) - then the H estimation would be more accurate

still even with those inacurracies the result looks quite good in my eyes - a few other modifications coulb de done (some operations did'nt need to be inside the for loop .. i counter is unnecessary...)

new code :

%% 2
f_s = 44.1*10^3; % sampling frequency
tau = 0.005; % tau
dt = 1/f_s; % delta time

a = [1, 1/tau]; % low pass filter
b = [1/tau];
sys_low = tf(b,a);
a = [1, 1/tau]; % high pass filter
b = [1, 0];
sys_high = tf(b,a);

nf = 48;
f = logspace(0,4,nf);

M_1 = zeros(1,nf); % the magnitude of low pass
M_2 = zeros(1,nf); % magnitude of high pass
P_1 = zeros(1,nf); % the phase of low pass
P_2 = zeros(1,nf); % phase of high pass
H1 = zeros(1,nf); % initializing the complex gain for low pass
H2 = zeros(1,nf); % initializing the complex gain for high pass

for k = 1:nf 
    w = 2*pi*f(k); % angular frequency
    T = 1/f(k);
    t_end = max(10*T,8*tau); % find the steady state function
    t = 0:dt:t_end;
    vin_sin = sin(w*t); % 
    vin_cos = cos(w*t); % 

    [y1,t_out_1] = lsim(sys_low,vin_sin,t); % low pass for f using lsim
    [y2,t_out_2] = lsim(sys_high,vin_sin,t); % high pass for f using lsim
    y1 = y1';
    y2 = y2'; % make same dimension
    idx = round(length(t)*0.8); % use the last 20% of the data
    t_segment = t(idx:end)'; % make the segment of the last period
    y1_segment = y1(idx:end)';
    y2_segment = y2(idx:end)';
    sin_segment = vin_sin(idx:end)';
    cos_segment = vin_cos(idx:end)';    
    samples = numel(t_segment);
    H1(k) = 2/samples*(sum(y1_segment.*sin_segment) + 1i*sum(y1_segment.*cos_segment)); % the complex gain
    H2(k) = 2/samples*(sum(y2_segment.*sin_segment) + 1i*sum(y2_segment.*cos_segment)); % the complex gain
end
M_1 = 20*log10(abs(H1)); % the magnitude of steady state output of low pass
M_2 = 20*log10(abs(H2)); % the magnitude of steady state output of high pass
P_1 = (angle(H1)/pi); % the phase of the steady state output of low pass
P_2 = (angle(H2)/pi); % the phase of the steady state output of high pass
figure;
subplot(2,2,1);
semilogx(f,M_1); % plot the magnitude over frequency
% xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Figure 1');
grid on;
subplot(2,2,2);
semilogx(f,M_2); % plot the magnitude over frequency for high pass filter
% xlabel('Frequency (Hz)');
% ylabel('Magnitude (dB)');
title('Figure 2');
grid on;
subplot(2,2,3);
semilogx(f,P_1); % plot the phase over frequency for low pass filter
xlabel('Frequency (Hz)');
ylabel('Phase (π radians)');
title('Figure 3');
grid on;
subplot(2,2,4);
semilogx(f,P_2); % plot the phase over frequency for high pass filter
xlabel('Frequency (Hz)');
% ylabel('Phase (π radians)');
title('Figure 4');
grid on;
% sgtitle('Phase and Magnitude Response of Steady State for Low and High Pass Filters');
sgtitle('Bode Plots');

[Paul]

Hi Mathieu,

If simulating with a unit-amplitude, complex sinusoid input we only need to ensure the output is in steady state to compute the gain and phase directly.

f_s = 44.1*10^3; % sampling frequency
tau = 0.005; % tau
dt = 1/f_s; % delta time
a = [1, 1/tau]; % low pass filter
b = [1/tau];
sys_low = tf(b,a);
t = (0:dt:10)';  % assume 10-seconds is long enough to get to steady state
f = 10;
h = freqresp(sys_low,2*pi*f); % truth, for comparison
u = exp(1j*2*pi*f*t);
y = lsim(sys_low,u,t);
[abs(y(end)),abs(h)]  % gain

ans = 1×2

0.9540 0.9540

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[angle(y(end)./u(end)),angle(h(end))] % phase

ans = 1×2

-0.3044 -0.3044

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Mathieu NOE]

hi @Paul

ok yes that makes sense

the only benefit of my proposal is that is something you can do in real world (and it averages the computation of the TF over multiple periods so the noise effects are reduced (again in real life))

Answer 2

[Andrew Ouellette]

Hi @Zeyuan,

Have you tried using the bodeplot function?

f_s = 44.1*10^3; % sampling frequency
Ts = 1/f_s; % sample time
tau = 0.005; % tau
a = [1, 1/tau]; % low pass filter
b = [1/tau];
sys_low = tf(b,a);
a = [1, 1/tau]; % high pass filter
b = [1, 0];
sys_high = tf(b,a);
sys = c2d([sys_low sys_high],Ts); % I assume your transfer functions are in continuous time from your notation
h = bodeplot(sys);
h.PhaseUnit = "rad";
grid on;