Electric circuits 101 RC and RL circuits

RC and RL are one of the most basics examples of electric circuits and yet they are very rich in content. Let’s examine each one carefully

Here is the schematics of an RL circuit

The components used in an RL circuit are a DC voltage source (V1) a resistor (R1) and an inductor (L1). The inductor in the schematics is used to represent physical characteristic of the circuit. If I had to describe what the inductor represents I would say it represents the ‘electrical inertia’ of the circuit. As currents flows into the circuit it generates a magnetic field, that change in the magnetic field causes a change in the flux of the field concatenated to the circuit, this in turn, by the Faraday-Neumann-Lenz law generates a voltage in the circuit that is opposite to the voltage that is generating the magnetic field (V1). This is the reason because the current in the circuit will not jump immediately to its full value of V/R given by Ohm’s law. The current will eventually reach the value V/R after an infinite amount of time, but for practical purposes we can consider it at that level after some time has passed (namely 4 times tau = L/R should be sufficient).

It so happens that, given the materials and the geometry, the inductance coefficients of the circuit I am modelling is 0.0229 H. Note that inductance depends only on these two parameters (materials and geometry). I calculated the coefficient using FEMM since the analytical calculation was not that easy to do. You can download FEMM here: http://www.femm.info/wiki/HomePage . This software is useful for many electromagnetic simulations and fairly user-friendly.

What we would expect is that the current will obey the following differential equation given by Ohm’s law at each point in time:

The equation above yields the following solution

Now, the same phenomenon (reversed) happens when you open the circuit, the inductor opposes the change in flux by generating a voltage that tries to keep up the shrinking flux of the magnetic field through the circuit. That means that, if you were to switch on and off the current at particular intervals you would obtain this behaviour:

As you can see, at time 0 the circuit is closed and the current is reaching its theoretical value, after 4 times tau (time constant of the circuit) the circuit is opened again and the current exponentially decays to zero. By knowing tau, one could find the frequency at which the switch should click in order to maintain this behaviour. This process yields the on and off cycle (in green). Where the green function hits the x axis the switch clicks and opens/closes the circuit.

Here below you can find the script I used to produce this graph and model the circuit.

If you are not satisfied with this crude simulation, perhaps you could use LTspiceIV which is a great free simulation software for electrical and electronic circuits. I used it to produce the schematics above and to simulate the behaviour of the current in the first 40ms after startup. Below is the LTspice simulation graph

As you can see by plotting the current through the resistor in LTspice, it seems our calculation were right!

Now let’s move on to the RC circuit, here it is:

Aside from the voltage source, the RC circuit is composed of a capacitor and a resistor, we therefore assume that self-inductance is negligible. When we close the circuit, there is no inductance here so the current will just jump to the V/R value and since the capacitor is charging up and building a voltage, we can expect the current at the resistor to drop as times goes by. The differential equation we have now is the following:

By solving it, one finds:

Therefore

The current through the resistor drops exponentially, depending on a time constant tau = 1/RC. This behaviour is indeed what we find by running the simulation in Python and LTspice

Python code used for the RC circuit:

PDEs time again: the Transport equation

Apparently, it looks like I’m all focused on PDEs at this time of the year. There’s a good reason though! I find PDEs extremely fascinating in that they allow you to model complex behaviours of objects such as waves propagating and strings vibrating and so on. If not for their beautiful mathematical content you must love them for their applications in physics!
Take the transport equation, imagine you’d like to model the behaviour of an ocean wave. Well, in an ideal world, once generated, our wave would go on forever at a constant speed, say k. In order to find the function that represents our wave at each point at every time t, one should solve the following problem



It turns out solutions are all in the form of:



For those of you already familiar with the wave equation, this should look extremely familiar. It is indeed the equation of a wave travelling towards the positive x at speed k
If we choose a Gaussian function for instance, then we would obtain the following:



and by running the same code we used for the heat equation, here is what we get:

Since the code is pretty much the same as the heat equation I will not post it here. However, it is downloadable via dropbox by clicking here if you’d like to check it out.
Now a question should arise. Our wave is an ideal wave, it does not take into account friction and if you think of an ocean wave, it dissolves itself quickly as it meets the beach. Can we model this with the transport equation? Sure! We can use the transport equation with exponential decay. The problem we need to solve now looks a little different:



where gamma is a constant which relates to the decay rate. The general solution now looks something like this



Note that it decays and tends to zero as time tends to infinity and that it is still a wave travelling ‘to the right’. If we then use the same Gaussian as before, we obtain a much more realistic wave. Here is a short clip generated using the code at the bottom of the page (the animation is slower since I set dt=0.01 for a better resolution).

You can play with the parameters I set and see how the behaviour of the wave changes!

Heat Equation part 2 a slight modification

While writing the code for the previous post I slightly modified the code in order to add 2 ‘peaks of heat’. Now, by looking at what I wrote I am not quite sure that the solution is consistent (physically speaking). At first glance we can see that there are three ‘heat spots’. If you imagine of heating the rod in two spots and cooling it in between, the situation is the following:



By setting L = 3*pi and making some other small changes we obtain:

Physically we could interpret the result as follows:
the heat flows spontaneously towards areas where the temperature is lower, therefore the temperature should go up in the middle of the rod (since it has been cooled down and has a lower temperature compared its surroundings) and go down in the sides. Assuming that the left and right boundaries of the rod are kept at constant temperature through an ideal thermostat.


Why not  checking if this solution is mathematically and physically consistent?

The Heat Equation: a Python implementation

By making some assumptions, I am going to simulate the flow of heat through an ideal rod.
Suppose you have a cylindrical rod whose ends are maintained at a fixed temperature and is heated at a certain x for a certain interval of time. Suppose that the temperature in each section with infinitesimal width dx is uniform so that we can describe the temperature in the rod using a function of only x and t.


Mathematically speaking, problem we are now facing is the following:



where k is a constant called thermal diffusivity and is different according to the different materials. By using the method of separation of variables, we can find the solution we need and by applying the initial conditions we find a particular solution for f(x) = sin(x) and L = pi
Our solution looks something like this:

Now we only need to evaluate our function at each x and t. Remember that if u(x,y) is differentiable, then:



holds. We can throw out the last term and approximate our function using the above relation since partial derivatives of u must exists and we can easily get them. Thinking about it, the second term is useless too, since we are not moving along the x axis, therefore we are left with the following:



By using matplotlib and the animation function I generated this short clip:

Here is the python implementation of the solution and the code used to graph the solution evolving with time.


Hope this was useful!

Estimating data parameters using R

Say we have some data and we are pretty confident that it comes from a random variable which follows a Normal distribution, now we would like to estimate the parameters of that distribution. Since the best estimator for the population mean is the sample mean and the best estimator for the variance is the corrected variance estimator, we could use those two estimators to compute a point estimate of the parameters we need. But, what if we would like to have a rough idea of what could be the range of those parameters within a certain level of confidence? Well, then we would have to find an interval that contains the parameters at a, say, 5% confidence level.

In order to do this, since the variance is unknown and needs to be estimated, we use the Student-t distribution and the following formula for the two sided interval:

In the same way, we could find unilateral boundaries for the value of the population mean, simply by adjusting the formula above:

Furthermore one could visualise the outcome by plotting the selected regions on a standard t-Student with n-1 degrees of freedom:

1) Confidence interval2) Upper bound3)Lower bound

Here below is he implementation in R of the process described above:

By clicking here you can download the code of the plotting function I used to generate the plots above.

How to make a rough check to see if your data is normally distributed

Now then, in the previous article I wrote about hypothesis testing with data that is normally distributed, in this article I’m going to post some quick test you can do to check if it is fairly safe to assume your data is normal. I would like to highlight the fact that you can never be 100% sure that your data follows a particular distribution, however, in order to be able to say something about the phenomenon you are trying to understand, it is often practical to approximate its distribution with the best fit according to the measurements (the samples). At the very least this is a starting point.

First, I should point out that if you have a small sample (< 10 samples) then it is really hard to draw some conclusions, but if the sample is large enough (as a broad and very generous rule of thumb > 30) then we can consider some indices and plots.

Skewness:
First of all, if your data is normally distributed, then it should be approximately symmetric. In order to asses this characteristic you can either use a boxplot or the skewness index g4. However I suggest to use both. Suppose you have two datasets as below:

If we plot the boxplots we obtain:

Remembering that the black thick line is the median and that the coloured part of the boxplot contains 50% of the data, we can already see that the distribution on the left is approximately symmetric, while the one on the right has a right tail (right skew).
Then, by running the skewness function we check the numbers in order to double check our guess

Since the skewness of the first dataset is close to zero it seems reasonable to assume the first dataset is normally distributed. Note that the beta simulated dataset has a positive skewness significantly higher than 0. A final check at the histograms can help a bit more

Ok, the histogram may cast some doubt on whether the first dataset is normal, however it certainly helps us ruling out normality for the second dataset

Finally, we run the kurtosis test, which gives us some information on the shape of the bell curve

We are now in a pretty safe position to assume that the first dataset is normally distributed, since the kurtosis is pretty close to 3, and we can see that the second seems to have a “peak” which is indeed what the histogram tells us. Having so many samples in the second dataset makes us pretty confident in ruling out normality, however it could also not be the case, for instance, for some reason, say the instruments with which we made the measurements were biased, we could have ended up with biased data.

Edit: Last day, while I was writing this article I forgot a very important plot: the qqplot! The function qqnorm lets you compare  the quantiles of your data with those of the normal distribution, giving you another bit of information that could point you in the right direction.
By calling the function qqnorm(data_beta) and qqnorm(data_norm) we obtain the following plots

Now, ideally, what you would expect if your data is approximately normally distributed, is that the quantiles of the sample match the theoretical samples, therefore you should get a 45 degree line. Can you guess which one of the two plots is the one generated by qqnorm(data_beta)? ;)

Hypothesis testing on normally distributed data in R

Hypothesis testing is a useful statistical tool that can be used to draw a conclusion about the population from a sample. Say for instance that you are interested in knowing if the average value of a certain parameter differs significantly from a given value within a well defined confidence level: you would probably set up your test like this:

Note that the two hypothesis above are mutually exclusive.

Now  what can you do to check which one is more probable?
First, we should check what kind of distribution our data follows. Roughly speaking, if the number of samples is > 30 we can plot a histogram to get a visual grasp of the distribution and then run a few simple function, to assess the skewness and kurtosis of the distribution. If these parameters are close to those of a Normal distribution, then we could assume that the data comes from a Normal distribution. This assumption is not 100% sure but it is reasonable if the above parameters are close to those of a Normal distribution. If you are not sure, the other option is to run some normality tests or gather more data. Read how to make a rough check to see if your data is normally distributed.

Now that we are confident that our data is Normally distributed and hopefully the samples are independent and identically distributed, we can estimates of the sample mean and as far as the variance is concerned:
- If the population variance is known, we can use the Normal distribution in R
- If the population variance is unknown, then we should use a t-distribution with n-1 degrees of freedom (n is the number of observations). As n tends to infinity t-distribution tends to a Normal, therefore if n is sufficiently large (say > 60) you could approximate the t-distribution with a Normal one.

Now, by fixing a certain alpha (confidence level), we decide to reject the null Hypothesis (H0) if:

Where is the sample mean, the standard deviation and z the percentile of the Normal (or Student) distribution.

Again, roughly speaking, the idea behind the test is that we should reject the null Hypothesis if the data shows that it is highly unlikely H0 to be true. For the right tail test, we should reject H0 if the statistics above (without the absolute value operator) is to the right of the red line, in the cyan shaded region:

Aside from this test you can also make tests to check whether the population mean is higher or lower than a certain value using the same process. Here below is the implementation of the test in R assuming that the variance is unknown (and therefore using the t distribution):

In this case we can say that at 5% confidence level we do not reject the null Hypothesis in all three cases.

In case you’d like to run a z-test because you know the population variance, by clicking here you can download the script which is using the normal distribution. Essentially you only have to replace the functions related to the student distribution in the script above with those related to the Normal. Finally, here is the code used to make the plot above.

Biot-Savart law: magnetic field of a straight wire

Magnetism and magnetism related phenomena are fascinating almost for everyone, in fact, I remember being a intrigued by the interaction between magnets since I was a little kid. On the other side, Python is great for plotting vector fields, although now that I am using also Matlab I found out most functions in numpy and matlab are similar to those used in Matlab (I don’t know who got inspired by who though).

Here is a small sketch of a vector field. Remember the Biot-Savart law for a straight infinitely long wire, which (assuming the current is flowing upwards) looks something like this:

Below is a 3d representation of the vector field and the Python code used to generate it. The tricky part (if I’m allowed to say so) is that you need to convert the value of the vector field from cylindrical coordinates into cartesian coordinates in order to plot it. Perhaps there’s a way to plot a vector field using also spherical and cylindrical coordinates however I still do not know how to do that.

Hope this was interesting! :)

Spectrogram representation of the B note from a guitar

A spectrogram is a visual representation of the spectrum of frequencies as they change with respect of another variable (in our case time). Matlab offers a great function to plot the frequencies against time, therefore I decided to try it on the recordings I used in the previous post. Here is the result I got

Now, some comments. This representation of the frequencies tells us that the stronger components of our signal are below 1.5 KHz (1500 Hz) which is what we expected, however it looks like there are no other significant frequencies above 7 KHz and this is weird since in our FT of the signal we found higher frequencies (although with lower magnitude). Notice how the the higher the magnitude in the FT, the longer the frequency appear in the spectrogram.

Now, by using a simple tool I generated another spectrogram which seems to contain also the higher frequencies the one in Matlab did not have.

This spectrogram is better in tune with our FT of the signal, since it shows that our signal has frequencies up to 15 KHz as the FT showed

The program I used to generate the second spectrogram is Spek, it is free, and fast, however it has very few options. Below are some useful links

http://spek.cc/

http://en.wikipedia.org/wiki/Spectrogram

Applying Fourier to a real signal: guitar musical note B

I recorded a sound from the second string of my electric guitar (musical note B) and made a .wave file out of it. By using Matlab I tried to separate the sound of the note from the background noise and made some experiments.

At first I wanted to use Python for doing this, since I’m more knowledgeable about that than Matlab, however it looks like the most common ways to open a .wave file (scipy.audiolab and similar) are not available for Python 3 yet (either that or I could not find anything). Matlab provides simple functions for handling .wave files and great graphing tools but I could not find something to convert bins into frequencies as the function in numpy does, therefore I’m not sure the frequencies in the graph are correct.

Anyway, you can download the recording at this link, and here is the time domain signal representation:

By taking the FFT of the signal and filtering out the frequencies we don’t need, we filter out the noise and/or sounds we don’t need. My filtering rule lets through all the frequencies above 200 and below 1600 Hz. Of course this is a very naive way of filtering, however that’s a start.

You can play around and block out one, two, three peaks a time and hear the difference between the original recording and the edited one. Now, ideally B of the guitar should have a frequency of about 249 Hz (more or less), however it looks like the frequencies between 650 and 1020 Hz in our frequency domain representation are the main components of the note. Perhaps that’s because I messed up the bins. Anyway I’m glad that by applying a so simple filtering rule I’ve been able to remove the pitch sound at the beginning and smoothing out the musical note. Finally I also found out that Matlab has great publishing tools that enable you to print out an HTML file of your code, the results and graphs all together. Here is the code I used:

Calculate Fourier transform of a gaussian using Matlab

While wondering around in the Matlab documentation I found out there is a simple way to calculate the Fourier transform of any function using Matlab. Furthermore by using the function pretty() you can print out the result in a more friendly manner.

Let’s calculate for example the Fourier transform of the function u(x) = e^(-x^2). This calculation could also be done on paper: by applying the definition I wrote in the earlier article one finds out that:

Now, by using some properties of the transform (you can check them at http://en.wikipedia.org/wiki/Fourier_transform), we can play around and plot some functions. Here is the original function and its FT:

Let’s take the FT of e^(-(3x)^2)

It’s worth noting that when the function gets thinner, its FT gets more spread and viceversa. This phenomena is useful to understand the uncertainty principle.

This you tube video explains the link between this mathematical tools and physics.
https://www.youtube.com/watch?v=V7UNvDN_EZo

If you’d like to know more about the FT, check the wikipedia page:
http://en.wikipedia.org/wiki/Fourier_transform

Here is the Matlab code I used

Filtering an ideal signal with FFT

Hi everyone! Remember that in the last article I wrote that you can use the FFT to clean a signal from background noise? Well here is an example of signal filtering.

We start by generating a signal and then add some random noise using the random number generator in numpy. Here is the noisy signal:

Now, our signal is made up of two main frequencies: 20 and 30 Hz while the rest is mainly background noise. You can check this in the FFT of the signal:

Suppose we’d like to isolate the 20 Hz bit, that means we’d have to set to 0 each value in the Spectrum which is greater or lower than 20 Hz. Let’s block out everything outside 18 and 22 Hz and then apply the inverse FTT. Here is the signal we found (blue) compared to the exact signal we wanted to find (red). At the bottom of the graph there is the error between the two signals, notice that it is not that high, (at worst we missed 0.08, however on average we missed less).

Here is the python code I used to make this.

Have fun!

Fourier Transform: first trials

The Fourier transform is surely one of the most amazing pieces of maths applied to the sciences. For starters it is defined as follows:

The definition is fine provided that you can take the integral of u(x) over R. (Technically you should say that u(x) belongs to L1(R). But aside from these technicalities, how awesome is that definition? Calculating the FT of a function with pen and paper isn’t that straight forward at all, except in some cases. It involves mainly using complex analysis techniques. If you are familiar with calculating “standard integrals” it is rare that you can apply the same techniques here.

But how can we apply this to something? Well it turns out that if you have a signal u(x) in the time domain, its Fourier transform is the signal representation in the frequency domain. Essentially the FT enables you to decompose your time domain signal into all the frequencies that make up that signal. A first simple application is deleting noise from a signal. By checking the signal spectrum (its FT) you can choose which frequencies you need to delete from the signal in order to remove the noise. Then you just apply the inverse FT and get back your time domain signal cleaned up (more or less). Now obviously the math behind this is hairy, however this should not scare you at all, and if you are interested in these subjects I suggest you to check out some resources to get a better grasp of how all of this works.

In order to compute the FT of a signal with Python we need to use the ftt function built in into Scipy. FTT stands for Fast Fourier Transform and it is a brilliant algorithm to speed up the otherwise very long calculation of the FT of a real signal.

The signal I’ve made is a linear combination of the sine and cosine function, pretty straightforward here is its representation in the time domain:

A closer look to our signal:

Now by taking the FFT we get a list of complex numbers. Then we take the absolute value of those numbers and plot them to represent the spectrum:

Eventually, we just need to fix the lower axis and switch from bins to Hz. Since the plot is symmetrical we can take out the right half without losing information. As you can see below, we have three spikes in the frequency domain. Those three peaks corresponds to the three frequencies we set earlier in the signal (85,50,75 Hz).

 

Hope you enjoyed this, soon I’ll make a post about applying this to a real signal, such as the recording of guitar notes. Below is the code I used to plot the data above. Have fun!