Trouble simulating an RC low and high pass filter with trapoziodal aproximation

Hi,

I have been working to recreate the example code found in Andrew Simper's 2020 lecture at Cytomic but have run into a roadblock. I can't seem to get a reasonable capacitor approximation and am not sure why. You can see the relevant code on page 45 of the PDF here (lecture). I've attached my code below as well.

If you have any suggestions, please let me know. Thank you for your time.

%% circuit variables
Is=1e-12;   % Saturation current (in Amperes)
n = 2;         % Ideality factor (typically between 1 and 2)
temp = 300;       % Temperature in Kelvin
k = 1.380649e-23;  % Boltzmann constant (J/K)
q = 1.602176634e-19; % Electron charge (C)

% Calculate thermal voltage
vt=(k*temp)/q;  % Thermal voltage (in Volts)

R1=2.2e3;
R2=6.8e3;
C1=0.47e-6;
C2=0.01e-6;

%% simulation loop
fs=48e3; %sample rate
Fd=1000; %input frequnecy
T=(1/1000)*10;
time=0:1/fs:T; %create time vector
vs=2.*sawtooth(Fd*2*pi*time); %create sawtooth input

%figure,plot(time,vs), title('input wavefrom (voltage v time')

gr1=1/R1;%deinfe gr1 as the reciprocal of r1 for easier eq formation
gr2=1/R2;
m=1/2; %define division integral aproximation
h=1/fs; %define hight of the traps
gc1=C1/(h*m);
gc2=C2/(h*m);
ic1eq=0;
ic2eq=0;


%initalize vectors
v1=zeros(1,length(vs));
v2=zeros(1,length(vs));
v3=zeros(1,length(vs));

va1=0;
va2=0;
va3=0;

for k=1:length(vs)
    va1=vs(k);
    v1(k)=va1;
    va2=(gr1*va1+gc1*va3-ic1eq)/(gr1+gc1);
    v2(k)=va2;
    va3=(gc1*va2+ic1eq+ic2eq)/(gc2+gc1+gr2);
    %va3=(gc1*ic2eq+gr1*(ic1eq+ic2eq+gc1*va1))/(gr1*(gc2+gr2)+gc1*(gc2+gr1+gr2));
    v3(k)=va3;
    ic1eq = ic1eq + (1/m)*(gc1*(va3-va2)-ic1eq);
    ic2eq = ic2eq + (1/m)*(gc2*(va3)-ic2eq);
end
figure, plot(time,v3)
title('input wave form vs output'), xlabel('time'), ylabel('voltage')
hold on
plot(time,v1), legend('v3','v1')

[Umar]

Hi Dayan,

Could you be more specific explaining, “I can't seem to get a reasonable capacitor approximation and am not sure why.” It would help to guide you better and provide suggestions regarding this issue.

[Dayan]

Yes. Sorry for the lack of specificity, and thank you for your help. I have been staring at this problem for a while.

My output signal, which you can see plotted in the post as v3, is very different from the waveform Andrew obtained in his presentation and what I know a high/low-pass filter should look like. It has a very low amplitude and does not filter or smooth the sawtooth how I would normally expect a filter to do so.

Please let me know if this clarifies any of your questions.

[Umar]

Hi Dayan,

Thanks for clarification. Based on your current setup, potential adjustments might involve refining the equations governing `v2` and `v3` updates within the simulation loop. Here's a hypothetical adjustment focusing on refining the update equation for `v3`:

for k = 1:length(vs)

    va1 = vs(k);

    v1(k) = va1;

    va2 = (gr1*va1 + gc1*va3 - ic1eq) / (gr1 + gc1);

    v2(k) = va2;

    va3 = (gc1*va2 + ic1eq + ic2eq) / (gc2 + gc1 + gr2);

    % Ensure this equation reflects the intended filter behavior

    v3(k) = va3; 

    ic1eq = ic1eq + (1/m) * (gc1 * (va3 - va2) - ic1eq);

    ic2eq = ic2eq + (1/m) * (gc2 * va3 - ic2eq);

end

Make sure that your simulation accurately reflects the desired filtering behavior by reviewing the integration method, component values, and simulation loop structure. By systematically verifying and adjusting these aspects, you can align the simulated output `v3` with the expected characteristics of a high/low-pass filter as per your requirements. Please let me know if you have any further questions.

Answer 1

[Divyajyoti Nayak]

Hi @Dayan,

The reason your plot is not matching to the one in the lecture is because you are calculating ‘va2’ with old value of ‘va3’. Also, you have commented out the right equation for ‘va3’. You should calculate both ‘va2’ and ‘va3’ with respect to ‘va1’ only, or you can calculate one of them with ‘va1’ and then use this new value to calculate the other.

To fix your plot you can just switch the order of calculating ‘va2’ and ‘va3’.

for k=1:length(vs)
    va1=vs(k);
    v1(k)=va1;
    
    %va3=(gc1*va2+ic1eq+ic2eq)/(gc2+gc1+gr2);
    va3=(gc1*ic2eq+gr1*(ic1eq+ic2eq+gc1*va1))/(gr1*(gc2+gr2)+gc1*(gc2+gr1+gr2));
    v3(k)=va3;
    %va3 is being calculated before va2 as it is dependent only on va1.
    va2=(gr1*va1+gc1*va3-ic1eq)/(gr1+gc1);
    v2(k)=va2;
    ic1eq = ic1eq + (1/m)*(gc1*(va3-va2)-ic1eq);
    ic2eq = ic2eq + (1/m)*(gc2*(va3)-ic2eq);
end

Here is the result after making the change, it is quite like the plot in the lecture.

Hope this helps!

[Dayan]

Thanks for the advice!