I have a differential equation that has been proven to be correct.
The transfer function obtained by Laplace domain analysis and Matlab freqs match up and all is well.
The problem is somewhere in my attempt to digitize the differential equation. I am required to do this using the trapezium rule to approximate the analog signal.
The goal is to get a digitized equation in the form of..
Vo^n+1 = alpha1(Vi^n+1) + alpha2(Vi^n) + alpha3(Vo^n)
Where each alpha value is made up of a different combination of T, R1, R2 and C all numerical parameters.
I am following some guidelines set by my lecturer but if the approach seems wrong please let me know.
Current incomplete solution..
Then I am trying to convert that equation to the Z domain and therefore use the freqz tool in Matlab. Please see the working out below.
Here is the Matlab plot that shows what has gone wrong. The reason for the multiple lines are varying values of capacitance for C.
The blue lines represent the validated analogue freqs(laplace domain transfer function) frequency response.
The green lines show my attempt at using the freqz(z domain discrete transfer function) frequency response.
If it is of any further interest, the code to generate those plots can be found below as well.
Any help would be sincerely appreciated, clutching at straws and the occasional hair at the moment!
%% Comparing the impulse amlitude response using two different methods % % Basically checking if freqs matches freqz clc; clear all; % Circuit Parameters R1 = 2.2e3; C = 0.01e-6; R2 = 10e3; fs = 10000; % Sampling frequency (Hz) T = 1/fs; % Timestep % Coefficients den = (T*R1 + T*R2 + 2*R1*R2*C); N = 1000; % Digital parameters y = zeros(N,1); y(1) = 1; result = zeros(N,1); % Analogue parameters fa = 0:fs/N:fs-fs/N; w = 2*pi*fa; j =sqrt(-1); lastValue = 0; clf; figure(1); hold on; for i=1:5 % Analogue Recalculations b = 1/(R1*C); a = (R1 + R2)/(R1*R2*C); Ha = freqs(b,[1 a], w); % Digital Recalculations alpha1 = (T*R2)/den; alpha2 = alpha1; alpha3 = (2*R1*R2*C - T*R1 - T*R2)/den; [Hd,fd] = freqz([alpha1 alpha1],[1 alpha3],N,fs); for n=1:N-1 currentValue = alpha1*y(n+1) + alpha2*y(n) - alpha3*lastValue; result(n) = currentValue; lastValue = currentValue; end C = 1.5*C; plot(fa,abs(Ha),fd,abs(Hd)),xlabel('Frequency'),ylabel('H(s)'); end hold off;