
LC LowPass Filter Analysis using the opensource computing language R .


If you are an Electronics Engineer working on circuit analysis and simulation, R may not be the first tool that comes to mind. In fact, there are excellent (and free) simulation tools available like the LTSpice package from Linear technologies for example. So why bother using something like R for this purpose? Well, this started as one of those “I'll do it just because I can" kind of project, but by the time I was done, I realized there are actually some very useful features in R that make this type of software an attractive option.
Figure 1 shows a typical LC filter topology. A circuit like this is commonly used as a noise filter to prevent switching power supply noise from reaching sensitive analog circuits. The inductors used in these applications are not ideal, so the R component models the small Inductor series resistance (sometimes referred to as “DCR” or DC resistance).
Figure 1  LC Filter Topology
A primary consideration in the design of this type of filter is obviously the 3dB corner frequency. One usually wants this corner frequency to be well below the switching power supply switching frequency so that the noise can be filtered effectively. A second consideration is the influence of the DCR series resistance in the passband peaking (amplification) this type of topology produces. While one wishes to keep the DCR at minimum to reduce voltage drops in series with the circuit, these lower values of DCR actually increase peaking.
let's see how this type of analysis can be done in R. I started by generating the input frequency range vector. For this purpose, I constructed the fseries() function that yields a vector of frequency values geometrically scaled. You just input the minimum/maximum frequencies and the number of points per decade and the vector is produced. We will evaluate the transfer function at each of these frequency values:
# fseries(min,max,perdec)
# This fuction returns a series of frequencies, geometrically distributed. Inputs:
# min: the first value
# max: the last value
# perdec: number of points per decade
fseries < function(min, max, perdec) {
expmin = floor(log10(min))
expmax = floor(log10(max))
step = 1/perdec
return(10^(seq(expmin,expmax,by= step)))
}
Since R handles complex numbers natively, writing a complex transfer function is a breeze. The following function computes the complex transfer function. Notice that I computed the series and parallel impedances first and then computed the transfer function as an impedance divider. I could have jumped straight into the wellknow Laplace H(s) transfer function for this type of filter, but I think this formulation is more generic and can be easily adapted to other, possibly more complex, filter topologies:
# Hf(f, R, L ,C)
# Transfer function for LC lowpass filter with series R.
# Returns (complex) transfer function value. Inputs
# f = Frequency in Hz. Can be a vector!
# R = Inductor series resistance in Ohm
# L = Series Inductance in Henry
# C = Parallel Capacitance in Farad
Hf < function(f, R, L ,C) {
s = as.complex(2i*pi*f)
Zser = as.complex( R + s*L)
Zpar = as.complex(1/(s*C))
return(Zpar/(Zpar+Zser))
}
So now we just have to create the series of frequencies we want to analyze and plot the results.
I'll be using the excellent ggplot2 package for plotting as the results are more to my liking than with the 'base' plotting system. Note that R also makes it easy to find things like the maximum peaking. For example:
library(ggplot2) library(scales) f = fseries(100,1e6,100) #100Hz to 1MHz, 100pts/decade R = 0.06; L = 1E6; C = 47E6 amplitude = 20*log10(Mod(Hf(f,R,L,C))) #amplitude in dB = 2olog(H(f)) df < data.frame(x = f, y = amplitude) p1 < ggplot(df,aes(x,y)) + geom_line(lwd = 1, col = 'blue') + scale_y_continuous(breaks=seq(65,25,by = 5)) + labs(x = "f [Hz]", y = "H(f) [dB]") + scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+ annotation_logticks(sides = "bt") p1
paste("Maximum Peaking [dB] :", max(amplitude))
## [1] "Maximum Peaking [dB] : 7.89803469994997"
let's find the frequency at which this maximum occurs:
indx = which(max(amplitude) == amplitude)
paste("Peaking Frequency [Hz] :", f[indx])
## [1] "Peaking Frequency [Hz] : 22387.2113856834"
We can also easily find the 3dB corner frequency:
indx = which(round(amplitude) == 3)
# After rounding, several values can match 3dB. Use the mean
paste("3dB Corner Frequency [Hz] :", mean(f[indx]))
## [1] "3dB Corner Frequency [Hz] : 35077.5119843054"
It's well know that a LC filter of this type should yield a 40dB/decade attenuation beyond the cutoff frequency. Let's confirm this:
f0 = 100e3
slope = 20*log10(Mod(Hf(f0,R,L,C)))  20*log10(Mod(Hf(10*f0,R,L,C)))
paste("Slope in the stopband [dB/decade] :", slope)
## [1] "Slope in the stopband [dB/decade] : 40.4329008350589"
One may be interested in knowing how the (undesirable) peaking changes with different values of R. We can use a 'for' loop for this:
L = 1E6; C = 47E6
df = data.frame(x = numeric(), y = numeric(), R = numeric())
for (R in seq(0.01,0.11, by = 0.02)){
amplitude = 20*log10(Mod(Hf(f,R,L,C)))
dfi < data.frame(x = f, y = amplitude, R = R)
df = rbind(df,dfi)
}
df$R = as.factor(df$R)
p < ggplot(df,aes(x,y, colour = R)) + geom_line(lwd = 1) +
labs(x = "f [Hz]", y = "H(f) [dB]") +
scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+
scale_y_continuous(breaks=seq(65,25,by = 5)) +
annotation_logticks(sides = "bt")
p
Zooming on the peaking area:
df = subset(df, x < 1e5 & x > 1000)
p < ggplot(df,aes(x,y, colour = R)) + geom_line(lwd = 1) +
labs(x = "f [Hz]", y = "H(f) [dB]") +
scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+
scale_y_continuous(breaks=seq(25,25,by = 5)) +
annotation_logticks(sides = "bt")
p
One strategy to dampen the peaking is to add a resistor in parallel with the inductor. See Figure 2 below. Let's modify the transfer function to accommodate this extra resistor 'Rp':
Figure 2  Modified Topology
# Hfp(f, R, Rp, L ,C)
# Transfer function for LC lowpass filter with series R.
# Returns (complex) transfer function value. Inputs
# f = Frequency in Hz. Can be a vector!
# R = Inductor series resistance in Ohm
# L = Series Inductance in Henry
# C = Parallel Capacitance in Farad
# Rp = Parallel resistor (parallel with inductor)
Hfp < function(f, R, Rp, L ,C) {
s = as.complex(2i*pi*f)
Zser = as.complex( ((R + s*L) * Rp)/( (R + s*L) + Rp))
Zpar = as.complex(1/(s*C))
return(Zpar/(Zpar+Zser))
}
And now we print the resulting transfer function:
R = 0.06; L = 1E6; C = 47E6; Rp = 0.1
amplitude = 20*log10(Mod(Hf(f,R,L,C)))
df1 < data.frame(x = f, y = amplitude, Rp = "None")
amplitudep = 20*log10(Mod(Hfp(f,R, Rp, L,C)))
df2 < data.frame(x = f, y = amplitudep, Rp = "0.1 Ohm")
df = rbind(df1,df2)
df$Rp = as.factor(df$Rp)
p2 < ggplot(df,aes(x,y, colour = Rp)) + geom_line(lwd = 1) +
labs(x = "f [Hz]", y = "H(f) [dB]") +
scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+
scale_y_continuous(breaks=seq(65,25,by = 5)) +
annotation_logticks(sides = "bt")
p2
As we can see, the peaking can be reduced significantly, but at a great penalty in stopband attenuation, which is why this method isn't always used.
Iterating over Rp:
R = 0.06; L = 1E6; C = 47E6
df = data.frame(f = numeric(), y = numeric(), Rp = numeric())
for (Rp in seq(0.05,0.5, by = 0.05)){
amplitude = 20*log10(Mod(Hfp(f,R,Rp,L,C)))
dfi < data.frame(f = f, y = amplitude, Rp = Rp)
df = rbind(df,dfi)
}
df$Rp = as.factor(df$Rp)
p < ggplot(df,aes(x = f,y = y , colour = Rp)) + geom_line(lwd = 1) +
labs(x = "f [Hz]", y = "H(f) [dB]") +
scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+
scale_y_continuous(breaks=seq(65,25,by = 5)) +
annotation_logticks(sides = "bt")
p
and zooming in:
df$Rp = as.factor(df$Rp)
df = subset(df, f < 1e5 & f > 1000)
p < ggplot(df,aes(x = f,y = y , colour = Rp)) + geom_line(lwd = 1) +
labs(x = "f [Hz]", y = "H(f) [dB]") +
scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+
scale_y_continuous(breaks=seq(65,25,by = 5)) +
annotation_logticks(sides = "bt")
p
I was pleasantly surprised by how easy it was to model this type of filter in R. Simulations ran very fast and it's relatively easy to iterate over component values. I could also imagine using R's powerful optimization tools to derive optimal component values in a circuit given a set of constraints. This flexibility, combined with ggplot's nice graphics outputs make R a tool to seriously consider for filter analysis.
Comments, questions, suggestions? You can reach me at: contact (at sign) paulorenato (dot) com