LoopFilterMatrix.m

%here we have a standard loop filter as:
%         V2         V3         V4
%   I1->------+--R2--------R3-----
%        |   |       |          |
%       R1  C2      C3         C4
%        |
%     V1 |
%       C1
%This progam is a test of the solution for the loop filter
%it uses sine wave where we can check the phase and amplitude, the circuit
%is meant to be dropped into a PLL time domain analysis program, later.
%The presentation, here, looks like standard matrix equations so we can
%understand what is going on.

clear
NumberSamples=5000;
SampleFreq=10e4;
dT=1/SampleFreq;
Bin=SampleFreq/NumberSamples;
%pre-allocate V for speed, V could be a 2D array, who cares
V1=zeros(1,NumberSamples);
V2=zeros(1,NumberSamples);
V3=zeros(1,NumberSamples);
V4=zeros(1,NumberSamples);
R1=1e3;
R2=1e3;
R3=1e3;
C1=1e-6;
C2=100e-12;
C3=1e-9;
C4=.5e-9;

freq=100;
I1=1e-4*cos(2*pi*freq*(1:NumberSamples)*dT);

%note the symmetry in the following matrix
A=[  1/R1+C1/dT  -1/R1            0                0;...
    -1/R1        1/R1+C2/dT+1/R2  -1/R2            0; ...
    0            -1/R2            1/R2+C3/dT+1/R3  -1/R3;...
    0            0                -1/R3            1/R3+C4/dT];
    
%we can do the following in MatLab much faster as an indexed approach,
%but the following is easy to impliment in C if we make our own matrix
%inverter, besides it is easy to see the relationship to a standard matrix
%statement.
    for SN=2:NumberSamples
        B=[        C1/dT*V1(SN-1) ;...
            I1(SN)+C2/dT*V2(SN-1) ;...
                   C3/dT*V3(SN-1) ;...
                   C4/dT*V4(SN-1)];
               
        x=A\B;%Simple step, but it is the heart of the process
        
        V1(SN)=x(1); 
        V2(SN)=x(2); 
        V3(SN)=x(3); 
        V4(SN)=x(4); 
    end
 plot(1:NumberSamples,V4,1:NumberSamples,1000*I1,'r')
%plot((1:NumberSamples)*Bin,20*log10(abs(fft(V3))))
grid on