Using Matlab to elegantly calculate frequency response

[Steven O.]

Another question related to my low-pass filter. Basically, I can
calculate the frequency response of the circuit, one frequency at a
time, using Matlab; what I can't figure out is how to use this method
to calculate the response for a whole bunch of frequencies at once.
Maybe someone can help me out. Here are the calculations for a single
frequency at a time:

EDU>> % The circuit analysis calculations are:
EDU>> % (Vin - V2)/R2 + (Vout - V2)/c2 = (V2 - Vout)/R4
EDU>> % (V2 - Vout)/R4 = Vout/c1

EDU>> % Re-arranging....
EDU>> % V2*(-1/R2 - 1/c2 - 1/R4) + Vout*(1/c2 + 1/R4) = -Vin/R2
EDU>> % V2*(1/R4) + Vout*(-1/R4 - 1/c1) = 0

So the goal is to set up a matrix: Vm * V = Vr, and then calculate:
V = inv(Vm) * Vr. Again, for a single frequency:

EDU>> f=2500;
EDU>> w=2*pi*f;
EDU>> R2=330;
EDU>> R4=1000;
EDU>> c1=1./(j.*w.*470e-9);
EDU>> c2=c1;
EDU>> Vin=1;
EDU>> Vm = [(-1/R2 - 1/c2 - 1/R4) (1/c2 + 1/R4) ; (1/R4) (-1/R4 -
1/c1)];
EDU>> Vr = [ -Vin/R2 ; 0 ];
EDU>> V=inv(Vm)*Vr;
EDU>> abs(V)
ans =
0.37972
0.050967

The second value for "ans" correctly matches the measured output from
the circuit at f = 2500 Hz; it also matches the result if I use the
following equation to calculate the output voltage:
EDU>> % Some pencil and paper work then yields
EDU>> % Vout = (Vin/R2) / [(R4/c1) * (1/R2 + 1/c2) + 1/c1 + 1/R2]

So far, so good. And, with the formula for Vout, I can use an input
vector f = [1 10 100 1000 etc.] for the input frequencies, and get a
vector for the output voltages. But the fact that I had to work out a
pencil and paper solution for Vout is not elegant. The elegant
solution was the one shown above -- set up the matrix for the circuit
behavor, and use V=inv(Vm)*Vr to let Matlab worry about the algebra.

The problem is when I try to combine the two approaches. That is, I
use the matrix algebra above, culminating with V=inv(Vm)*Vr; and I
also try to use an input vector for the frequencies, f = [1 10 100
1000 etc.], rather than a single value such as f = 1000. Matlab
chokes on the Matrix algebra plus input vector. Specifically, what I
get is:

EDU>> f=[1 100 1000 2500];
EDU>> w=2*pi*f;
EDU>> R2=330; R4=1000; c1=1./(j.*w.*470e-9); c2=c1; Vin=1;
EDU>> Vm = [(-1./R2 - 1./c2 - 1./R4) (1./c2 + 1./R4) ; (1./R4) (-1./R4
- 1./c1)];
??? Error using ==> vertcat
All rows in the bracketed expression must have the same
number of columns.

So right here, I'm stuck. I'll played with setting up the matrix in
all kinds of ways -- this posting would run about 500 lines if I
showed all of them -- but could not get this to work. Can someone
kindly show me how to use the input vector f, while retaining the
elegant matrix calculation approach?

In case anyone wants to see the circuit at hand, you can view it at:
http://www.oppenheimercommunications.com/OpAmpProblem.htm

Thanks in advance for all replies.

Steve O.


"Spying On The College Of Your Choice" -- How to pick the college that is the Best Match for a high school student's needs.
www.SpyingOnTheCollegeOfYourChoice.com

 

[Rune Allnor]

Steven O. skrev:

Another question related to my low-pass filter. Basically, I can
calculate the frequency response of the circuit, one frequency at a
time, using Matlab; what I can't figure out is how to use this method
to calculate the response for a whole bunch of frequencies at once.

Usually, nothing is done "at once". One gets an expression H(w) that
varies with angular frequency w, and computes the parameters
one at the time.

Maybe someone can help me out. Here are the calculations for a single
frequency at a time: ....
The second value for "ans" correctly matches the measured output from
the circuit at f = 2500 Hz; it also matches the result if I use the
following equation to calculate the output voltage:
EDU>> % Some pencil and paper work then yields
EDU>> % Vout = (Vin/R2) / [(R4/c1) * (1/R2 + 1/c2) + 1/c1 + 1/R2]

So you get a computed result that matches the measuerements?

So far, so good. And, with the formula for Vout, I can use an input
vector f = [1 10 100 1000 etc.] for the input frequencies, and get a
vector for the output voltages. But the fact that I had to work out a
pencil and paper solution for Vout is not elegant.

Hmm... depends on what you mean by "elegant." There are
program packages that do what you want. (P)SPICE used to
be one, for all I know it still is.

However, doing things by pen and paper is a way better
method of teaching than just plugging some numbers
into a simulator. For some reason I get the impression
you take some intro course in electronics.

The elegant
solution was the one shown above -- set up the matrix for the circuit
behavor, and use V=inv(Vm)*Vr to let Matlab worry about the algebra.

Well, I am sure somebody do things that way. I am not so
sure if that's a very good way of learning the trade.

Rune

 

[Ted Edwards]

Another question related to my low-pass filter. Basically, I can
calculate the frequency response of the circuit, one frequency at a
time, using Matlab; what I can't figure out is how to use this method
to calculate the response for a whole bunch of frequencies at once.
Maybe someone can help me out. Here are the calculations for a single
frequency at a time:

EDU>> % The circuit analysis calculations are:
EDU>> % (Vin - V2)/R2 + (Vout - V2)/c2 = (V2 - Vout)/R4
EDU>> % (V2 - Vout)/R4 = Vout/c1
...

Ever heard of APL? You could get a demo copy from IBM's site. APL2 is
probably the most powerful array processing language and it supports
complex numbers so no problem there. There's a bit of a learning curve
but once you have these expressions expressed for f{is}2500, you simply
let f be a bunch of frequencies (a vector) e.g. for f from about 12Hz to
300KHz with 1/4 octave spacing, set
f{is}10{times}2*.25{times}{iota}60
and run your function again.

Ted

 

[Ban]

Steven O. wrote:

Forget about Matlab, learn first to get the Transfer function. You started
right, the currents into each node are zero.
Wrong. Do you remember that the impedance of a capacitor is 1/jwC? So in the
second term the current would be: (V3-V2)/(1/jwC2). The third current can
not flow into the opamp, so it flows first through R4 and then C1 to gnd.
Now a series of R and C are?
V2/(R4 + 1/jwC1)
Do not forget the j being the imaginary operator sqrt(-1).We have now to make a loop where te sum of the voltages is zero:
V1 + (V2-V1) + (V3-V2) + V3 = 0

You can now rearrange those two equations, substitute V2 and finally
calculate V3/V1. This is called the transfer function and gives you values
for all frequencies. This is the goal.

Everything below is unnecessary, and in fact you have wasted your time with
it. When you think you got the transfer function right, post it here and we
can do the Laplace transform to get an easier way of writing it.

 

[Ban]

Steven O. wrote:

Forget about Matlab, learn first to get the Transfer function. You
started right, the currents into each node are zero.

Wrong. Do you remember that the impedance of a capacitor is 1/jwC? So
in the second term the current would be: (V3-V2)/(1/jwC2). The third
current can not flow into the opamp, so it flows first through R4 and
then C1 to gnd. Now a series of R and C are?
V2/(R4 + 1/jwC1)
Do not forget the j being the imaginary operator sqrt(-1).
We have now to make a loop where te sum of the voltages is zero:
V1 + (V2-V1) + (V3-V2)`+ V3 = 0

Here I made a mistake, the voltages have to be in the same direction as the
loop, so if V1 is positive, then the last term(V3) must be negative

 

[Fred Bloggs]

Another question related to my low-pass filter. Basically, I can
calculate the frequency response of the circuit, one frequency at a
time, using Matlab; what I can't figure out is how to use this method
to calculate the response for a whole bunch of frequencies at once.
Maybe someone can help me out. Here are the calculations for a single
frequency at a time:

EDU>> % The circuit analysis calculations are:
EDU>> % (Vin - V2)/R2 + (Vout - V2)/c2 = (V2 - Vout)/R4
EDU>> % (V2 - Vout)/R4 = Vout/c1

EDU>> % Re-arranging....
EDU>> % V2*(-1/R2 - 1/c2 - 1/R4) + Vout*(1/c2 + 1/R4) = -Vin/R2
EDU>> % V2*(1/R4) + Vout*(-1/R4 - 1/c1) = 0

So the goal is to set up a matrix: Vm * V = Vr, and then calculate:
V = inv(Vm) * Vr. Again, for a single frequency:

EDU>> f=2500;
EDU>> w=2*pi*f;
EDU>> R2=330;
EDU>> R4=1000;
EDU>> c1=1./(j.*w.*470e-9);
EDU>> c2=c1;
EDU>> Vin=1;
EDU>> Vm = [(-1/R2 - 1/c2 - 1/R4) (1/c2 + 1/R4) ; (1/R4) (-1/R4 -
1/c1)];
EDU>> Vr = [ -Vin/R2 ; 0 ];
EDU>> V=inv(Vm)*Vr;
EDU>> abs(V)
ans =
0.37972
0.050967

The second value for "ans" correctly matches the measured output from
the circuit at f = 2500 Hz; it also matches the result if I use the
following equation to calculate the output voltage:
EDU>> % Some pencil and paper work then yields
EDU>> % Vout = (Vin/R2) / [(R4/c1) * (1/R2 + 1/c2) + 1/c1 + 1/R2]

So far, so good. And, with the formula for Vout, I can use an input
vector f = [1 10 100 1000 etc.] for the input frequencies, and get a
vector for the output voltages. But the fact that I had to work out a
pencil and paper solution for Vout is not elegant. The elegant
solution was the one shown above -- set up the matrix for the circuit
behavor, and use V=inv(Vm)*Vr to let Matlab worry about the algebra.

The problem is when I try to combine the two approaches. That is, I
use the matrix algebra above, culminating with V=inv(Vm)*Vr; and I
also try to use an input vector for the frequencies, f = [1 10 100
1000 etc.], rather than a single value such as f = 1000. Matlab
chokes on the Matrix algebra plus input vector. Specifically, what I
get is:

EDU>> f=[1 100 1000 2500];
EDU>> w=2*pi*f;

This is your problem- line above tells MATLAB to make w a row vector.

EDU>> R2=330; R4=1000; c1=1./(j.*w.*470e-9); c2=c1; Vin=1;
EDU>> Vm = [(-1./R2 - 1./c2 - 1./R4) (1./c2 + 1./R4) ; (1./R4) (-1./R4
- 1./c1)];
??? Error using ==> vertcat

I would say the same thing if asked to do that crap with w a row vector.

All rows in the bracketed expression must have the same
number of columns.

So right here, I'm stuck. I'll played with setting up the matrix in
all kinds of ways -- this posting would run about 500 lines if I
showed all of them -- but could not get this to work. Can someone
kindly show me how to use the input vector f, while retaining the
elegant matrix calculation approach?

Set up a loop with index ranging from 1 to Len(f), make w a scalar and
build Vm index by index, where Vm=blahblah( some algebra with f).

[snip- nothing of interest]

 

[Steven O.]

Ban,

Thanks for taking the time to reply. However, you read a little too
hastily. In my full Matlab calculations, included in the original
posting, c2 is clearly *defined* as 1/(j*w*<value_of_cap2>). I'm
quite sure I got the transfer function correct, both because I checked
my algebra carefully, and because the resulting plot is identical to
the plot generated by a Multisim simulation of the circuit.

The problem -- again, I refer you to the original posting -- is that
the transfer function is complicated enough that extracting w as a
function of all the circuit parameters is not even remotely trivial.

Steve O.

Steven O. wrote:

Forget about Matlab, learn first to get the Transfer function. You started
right, the currents into each node are zero.

Wrong. Do you remember that the impedance of a capacitor is 1/jwC? So in the
second term the current would be: (V3-V2)/(1/jwC2). The third current can
not flow into the opamp, so it flows first through R4 and then C1 to gnd.
Now a series of R and C are?
V2/(R4 + 1/jwC1)
Do not forget the j being the imaginary operator sqrt(-1).
We have now to make a loop where te sum of the voltages is zero:
V1 + (V2-V1) + (V3-V2) + V3 = 0

You can now rearrange those two equations, substitute V2 and finally
calculate V3/V1. This is called the transfer function and gives you values
for all frequencies. This is the goal.

Everything below is unnecessary, and in fact you have wasted your time with
it. When you think you got the transfer function right, post it here and we
can do the Laplace transform to get an easier way of writing it.

"Spying On The College Of Your Choice" -- How to pick the college that is the Best Match for a high school student's needs.
www.SpyingOnTheCollegeOfYourChoice.com

 

[Rune Allnor]

Steven O. skrev:

The problem -- again, I refer you to the original posting -- is that
the transfer function is complicated enough that extracting w as a
function of all the circuit parameters is not even remotely trivial.

Eh... might this be part of the problem? There is no need to
"extract" w from the transfer function. It is the variable of the
tranfer function H(w). All the components make up the
parameters.

Rune

 

[Fred Bloggs]

Steven O. skrev:




Eh... might this be part of the problem? There is no need to
"extract" w from the transfer function. It is the variable of the
tranfer function H(w). All the components make up the
parameters.

Rune

Obviously you're missing the point that realization of elegance can only
be obtained through a distant non-trivialization of this extraction.

 

[Steven O.]

Steven O. skrev:

Eh... might this be part of the problem? There is no need to
"extract" w from the transfer function. It is the variable of the
tranfer function H(w). All the components make up the
parameters.

The goal of the exercise is to say, "In order for H(w) to equal 0.707,
w must be <value>." So what I want is w = f (Vin, Vout, component
values).

Steve O.

"Spying On The College Of Your Choice" -- How to pick the college that is the Best Match for a high school student's needs.
www.SpyingOnTheCollegeOfYourChoice.com

 

[Greg Heath]

Another question related to my low-pass filter. Basically, I can
calculate the frequency response of the circuit, one frequency at a
time, using Matlab; what I can't figure out is how to use this method
to calculate the response for a whole bunch of frequencies at once.
Maybe someone can help me out. Here are the calculations for a single
frequency at a time:

EDU>> % The circuit analysis calculations are:
EDU>> % (Vin - V2)/R2 + (Vout - V2)/c2 = (V2 - Vout)/R4
EDU>> % (V2 - Vout)/R4 = Vout/c1

Clarification:

1. The DC biases on the amp (V2 = 12V on pin 7 and V3 = -12V on
pin 4) have the same labels (V2,V3) as other nodes in the circuit.
Since this is an AC problem, lets call the DC biases V20 and V30
to eliminate that confusion.
2. The OP Amp is a voltage amplifier with a large negative gain.
Looking at the voltage gain between the OP amp output and inputs
yields

Vout = -A*(V3-Vout)
= A*V3 / (A-1)
~ V3 (for A >> 1).

3. If you assume Vout = V3, then R4 and C2 are in parallel and
the same current flows through R2, R4&C2, and C1. Therefore

(V2-Vin)/R2 = (Vout-V2)/(R4 + 1/jwC2) = jwC1 * Vout

Consequently

( 1/(R4+1/jwC2) ) * Vout - ( (1/R2)+ 1/(R4 + 1/jwC2)) * V2 = -Vin/R2
( 1/(R4+1/jwC2) -jwC1) * Vout - ( 1/(R4 + 1/jwC2)) * V2 = 0

EDU>> % Re-arranging....
EDU>> % V2*(-1/R2 - 1/c2 - 1/R4) + Vout*(1/c2 + 1/R4) = -Vin/R2
EDU>> % V2*(1/R4) + Vout*(-1/R4 - 1/c1) = 0

Which is not the same as mine.

So the goal is to set up a matrix: Vm * V = Vr, and then calculate:
V = inv(Vm) * Vr.

No!

Never again try to calculate a matrix inverse or ratio of determinants
to solve linear equations. Enter the command

help slash

then use

V = V\Vr.

See if this solves your problem.

Hope this helps.

Greg

Again, for a single frequency:

EDU>> f=2500;
EDU>> w=2*pi*f;
EDU>> R2=330;
EDU>> R4=1000;
EDU>> c1=1./(j.*w.*470e-9);
EDU>> c2=c1;
EDU>> Vin=1;
EDU>> Vm = [(-1/R2 - 1/c2 - 1/R4) (1/c2 + 1/R4) ; (1/R4) (-1/R4 -
1/c1)];
EDU>> Vr = [ -Vin/R2 ; 0 ];
EDU>> V=inv(Vm)*Vr;
EDU>> abs(V)
ans =
0.37972
0.050967

The second value for "ans" correctly matches the measured output from
the circuit at f = 2500 Hz; it also matches the result if I use the
following equation to calculate the output voltage:
EDU>> % Some pencil and paper work then yields
EDU>> % Vout = (Vin/R2) / [(R4/c1) * (1/R2 + 1/c2) + 1/c1 + 1/R2]

So far, so good. And, with the formula for Vout, I can use an input
vector f = [1 10 100 1000 etc.] for the input frequencies, and get a
vector for the output voltages. But the fact that I had to work out a
pencil and paper solution for Vout is not elegant. The elegant
solution was the one shown above -- set up the matrix for the circuit
behavor, and use V=inv(Vm)*Vr to let Matlab worry about the algebra.

The problem is when I try to combine the two approaches. That is, I
use the matrix algebra above, culminating with V=inv(Vm)*Vr; and I
also try to use an input vector for the frequencies, f = [1 10 100
1000 etc.], rather than a single value such as f = 1000. Matlab
chokes on the Matrix algebra plus input vector. Specifically, what I
get is:

EDU>> f=[1 100 1000 2500];
EDU>> w=2*pi*f;
EDU>> R2=330; R4=1000; c1=1./(j.*w.*470e-9); c2=c1; Vin=1;
EDU>> Vm = [(-1./R2 - 1./c2 - 1./R4) (1./c2 + 1./R4) ; (1./R4) (-1./R4
- 1./c1)];
??? Error using ==> vertcat
All rows in the bracketed expression must have the same
number of columns.

So right here, I'm stuck. I'll played with setting up the matrix in
all kinds of ways -- this posting would run about 500 lines if I
showed all of them -- but could not get this to work. Can someone
kindly show me how to use the input vector f, while retaining the
elegant matrix calculation approach?

In case anyone wants to see the circuit at hand, you can view it at:
http://www.oppenheimercommunications.com/OpAmpProblem.htm

Thanks in advance for all replies.

Steve O.


"Spying On The College Of Your Choice" -- How to pick the college that is the Best Match for a high school student's needs.
www.SpyingOnTheCollegeOfYourChoice.com

 

[Rune Allnor]

Steven O. skrev:

The goal of the exercise is to say, "In order for H(w) to equal 0.707,
w must be <value>." So what I want is w = f (Vin, Vout, component
values).

I don't understand what you mean. If you design the circuit from
stratch,
there is a spec somewher that says "H(w0) == 0.707" for some specified
frequency w0. If you are given a cirquit, with components in place,
what you want to do is to plot H(w) over a range of w.

If you want to take a cirquit diagram and solve for the w where
H(w)==0.0707,
I think you are in for a big job.

Rune

 

[Ban]

Ban,

Thanks for taking the time to reply. However, you read a little too
hastily. In my full Matlab calculations, included in the original
posting, c2 is clearly *defined* as 1/(j*w*<value_of_cap2>). I'm
quite sure I got the transfer function correct, both because I checked
my algebra carefully, and because the resulting plot is identical to
the plot generated by a Multisim simulation of the circuit.

The problem -- again, I refer you to the original posting -- is that
the transfer function is complicated enough that extracting w as a
function of all the circuit parameters is not even remotely trivial.

Steve O.

Steve, do not quarrel with me, look what is happening. You have made up one
loop with your equations, and you are able to express V2 as a function of
V3, but when you enter this into your first equation , you get what I wrote
for the loop 0 = 0 which is true, but doesn't help you to solve the transfer
function. :-(

In order to solve the thing you have to set up another equation to solve the
dependency of V2 from V1. We need another loop to determine the current
through C2 and then we can find the current through R2 as the sum of those
other two currents.

And please do not write c2 = 1/jwC2 which is misleading and sloppy, write
Xc2 or X2 instead, so your professor can understand you.

Tell us what your professor answered to you, maybe it was not so far away
from what I said?


You are still wasting your time. Matlab is only able to solve those
equations numerically, it doesn't have intelligence. This requires *your*
understanding.

 

[Ban]

Since you still think it is the problem with Matlab, I looked a little
closer. What you have is only of first order and that is clearly wrong.
So go back to your circuit and draw it on a paper. Look at the buffer. This
is a dependent voltage controlled voltage source. Make two loops where you
can superimpose the two voltage sources, treat the dependent source as you
learned it. Then derive the transfer function V3/V1(w) and finally replace
s=jw/wg (the sinosoidal steady state).

 

[Greg Heath]

It can be done by solving a quartic equation. See below.

The goal of the exercise is to say, "In order for H(w) to equal 0.707,
w must be <value>." So what I want is w = f (Vin, Vout, component
values).

Define absH2 = H*conj(H) and h2 = 1/2. Then you are looking for the
frequency at which absH2 = h2.

H(jw) = Vout/Vin will be complex and the ratio of two polynomials in
jw. Since Vout is the voltage across a capacitor, H(0) = 1 and
H(j*inf) = 0. Since there are two independent storage devices in
the circuit, the order of the polynomial in the denominator will be
2. Therefore, absH2 will be real and the ratio of two polynomials
in w. The order of the polynomial in the denominator will be 4.

Method1: Convert the equation absH2 = h2 into a quartic equation
that can be solved for w. This can be solved
analytically or numerically.

Method2: Use one of MATLABs minimization routines to minimize
(absH-h0)^2

Method3: Define a vector of frequencies. Next find the indices
for which max(find(absH2 >= h2)) and
min(find(absH2 <= h2)). Then use linear interpolation.

Using the voltage divider formula,

H(jw) = (1/jwC1) / ( (1/jwC1)+ R2 + 1/(G4+jwC2) )
= (G4 + jwC2) / ( G4+jwC2 + jwC1(G4+jwC2)R2+jwC1 )
--> G4/G4 as w --> 0
--> jwC2/(jwC1)(jwC2)R2 as w --> inf

Hope this helps.

Greg