SciLearn

Sunday, March 8, 2009

peak detection

I have been recently working on a Python program that detects peaks in an oscilloscope trace. With help from this post I found, I managed to find the best algorithm to deal with my situation.

Mathematically, finding the peaks (local maxima) of any graph is straightforward, as long as the function is simple enough. All one has to do is to find the x-values where the first derivative is zero, and then apply the second derivative test to differentiate the maxima (peaks) from the minima (troughs).

However, real data tend to be rather coarse, so peak detection can be rather difficult. From the article I mentioned earlier, this issue can be resolved by "derivative smoothing". The typical way of estimating the derivative of a function represented by a set of discrete points is by taking every pair of points and calculating the gradient. Derivative smoothing involves taking more than two points to calculate the derivative. The article says this can be achieved by taking a group of points and calculating the gradient using just the first and the last points in this group, but I decided to try something different: making a linear fit on this group of points and find the gradient of this fit.

I don't know if this is more accurate or appropriate, but I tried it out anyway, and so far it's working well. The only downside is that the calculations are much slower since fitting a set of data points is computationally intensive. In any of these peak detection methods, the precision of the peak is always limited by the discreteness of the values; to improve this precision, one can attempt to fit the peak itself with a fit curve, say a Gaussian (if symmetric), and use the fitted parameters to find a more precise value for the peak.

Friday, January 2, 2009

commutative power

The power operator ("a raised to the power of b") is certainly not commutative, i.e. in general a^b ≠ b^a. This should be no surprise, e.g. 2³ ≠ 3².

What I am interested is find pairs of natural numbers, in which a^b does equal b^a. I already know one pair: 2⁴ = 4² (in addition to the trivial pairs, e.g. 3³ = 3³), so I need to find out what's so special about this pair of numbers.

To do this, I break the equation into a series of steps:
2⁴ = 2^2×2 = (2²)² = 4².
Now the relationship becomes clear; I want a pair of numbers such that:
a^b = a^c×a = (a^c)^a = b^a
where c is another natural number that satisfies:
b = c×a and b = a^c
I can combine these equations into one:
c a = a^c
So how many (nontrivial) pairs of natural numbers satisfy this equation? By nontrivial, I exclude cases where c = 1.

Now, to solve this equation:
c a = a^c
c a - a^c = 0
(c - a^{c - 1}) a = 0
Either a = 0 (which is trivial) or c - a^{c - 1} = 0. Consider the latter nontrivial case:
c - a^{c - 1} = 0
c = a^{c - 1}
c^{1/(c - 1)} = a
Therefore, if a value of c is known, a can be found by taking the (c - 1)^th root of c. If c = 2, then a = 2^{1/(2 - 1)} = 2, and therefore b = c a = 4; this shows why 2⁴ = 4². For any other c > 2, (c - 1)^th root of c will not be an integer, so there are no other nontrivial natural-number combinations of a, b and c.

Well, looks like that answered the question. It seems 2⁴ = 4² is pretty special case, after all.

Friday, December 19, 2008

iron's odor

Well, just a short post on an interesting finding. I was curious why iron tends to produce strange metallic smells (it occurs to some other metals as well as blood), and it looks like a quick Googling gives me a neat explanation of this phenomenon. Go check it out.

It's funny that some common phenomena can be so easily ignored, and we never get to appreciate the science behind it.

Saturday, December 6, 2008

numerical pendulum

There was this "flash of idea" that I could solve "any" differential equations numerically, which implies I could "simulate" Maxwell's differential equations using a computer program. Sadly, that was too tall an order, so I could only manage to simulate the most basic of differential equations.

So far I managed to derive a simple method to simulate pendulums directly using the rigid pendulum equation (c is a negative constant)m
$\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} = c \sin \theta$
I do know Euler's method, which is really simple. But how do I apply it to second derivatives? That's the trick here.

I found that if I broke this equation into two by defining a variable ω,
$\begin{align*} \omega &= \frac{\mathrm{d}\theta}{\mathrm{d}t} \\ \frac{\mathrm{d}\omega}{\mathrm{d}t} &= c \sin \theta \end{align*}$
I can easily apply Euler's method using the finite differences approximation,
$\begin{align*} \frac{\Delta \theta}{\Delta t} &= \omega \\ \frac{\Delta \omega}{\Delta t} &= c \sin \theta \end{align*}$
Therefore, I can convert this into a form useful for programming,
$\begin{align*} \Delta \omega_i &= c \sin \theta_i \Delta t \\ \Delta \theta_i &= \omega_i \Delta t \\ \omega_{i + 1} &= \omega_i + \Delta \omega_i \\ \theta_{i + 1} &= \theta_i + \Delta \theta_i \end{align*}$
This is easily done by iteration; simply keep track of the index i and the values of ω and θ. By choosing the Δt to be small, one can approximate the pendulum with some precision. I have tested the method and it seems to be quite stable, though I'm not sure how good the precision.

Monday, November 24, 2008

101 posts?

Hey, I just realized that the last post was the 100^th one.

Happy Thanksgiving!
(Hopefully, I can write a full SciLearn post on Thursday.)

Saturday, November 15, 2008

wave interference

This post deals with the interference of light.

Light is a strange thing. It's made of photons, yet most of the time, it behaves as a wave. Now I won't go into the quantum mechanical details of light; I'm just going to talk about one particular nature of light: interference, which occurs whenever light behaves as a wave. I'll provide a systematic summary of wave interference and slit interference patterns.

The first step to interference is to begin with the basic equations of wave motion. Let's pick a simple sine wave in 1 dimensions:

d(x, t) = a sin ϕ = a sin[kx - ωt + ϕ₀]

where:
d = displacement;
a = amplitude;
k = angular wavenumber;
x = position;
ω = angular frequency;
t = time;
ϕ = phase;
ϕ₀ = initial phase.

When two waves interfere, their displacement superpose on each other, so they can simply be added to produce the resultant wave (for simplicity, amplitude of both waves are equal):

$\begin{align*} D (x, t) &= a \sin \phi_1 + a \sin \phi_2 \\ &= 2 a \cos \frac{\phi_1 - \phi_2}{2} \sin \frac{\phi_1 + \phi_2}{2} \\ &= 2 a \cos \frac{\Delta \phi}{2} \sin \frac{\sum \phi}{2} \\ &= 2 a \cos \frac{\Delta \phi}{2} \sin \langle\phi\rangle \end{align*}$

where
ϕ₁ = k₁x₁ - ω₁t + ϕ_0,1;
ϕ₂ = k₂x₂ - ω₂t + ϕ_0,2;
Δϕ = ϕ₁ - ϕ₂;
Σϕ = ϕ₁ + ϕ₂.
⟨ϕ⟩ = Σϕ / 2
The above equation involved the use of a "sum-to-product" trigonometric identity.

There are special scenarios under which the above equation simplifies:

Standing wave:
x₁ = x = x₂
k₁ = k = k₂
-ω₁ = ω = ω₂
ϕ_0,1 = 0 = ϕ_0,2
The displacement equation is:

$\begin{align*} & \quad \: D (x, t) \\ &= 2 a \cos \frac{\Delta \phi}{2} \sin \frac{\sum \phi}{2} \\ &= 2 a \cos \frac{\left(k x + \omega t\right) - \left(k x - \omega t\right)}{2} \\ & \quad \: \times \sin \frac{\left(k x + \omega t\right) + \left(k x - \omega t\right)}{2} \\ &= 2 a \cos \omega t \sin kx \end{align*}$
The boundary conditions for a standing wave is that the ends, at x = 0 and x = L, are fixed discontinuities (i.e. the ends of the standing wave cannot move):

D(0, t) = D(L, t) = 0 (for any t).
Solving this equation yields the allowable wavenumbers and wavelengths (λ) for standing waves:

$\begin{align*} D (L, t) &= 0 \\ A \cos \omega t \sin kL &= 0 \\ \sin kL &= 0 \\ kL &= n \pi, \quad n \in \mathbb{Z} \\ k &= \frac{n \pi}{L} \quad \Leftrightarrow \quad \lambda = \frac{2 L}{n} \end{align*}$
Coherent interference:
k₁ = k = k₂
ω₁ = ω = ω₂
The displacement equation is:

$\begin{align*} & \quad \: D (x, t) \\ &= A \cos \frac{\Delta \phi}{2} \sin \langle \phi \rangle \\ &= A \cos \frac{\left(k x_1 - \omega t + \phi_{0,1}\right) - \left(k x_2 - \omega t + \phi_{0,2}\right)}{2} \\ & \quad \: \times \sin \frac{\left(k x_1 - \omega t + \phi_{0,1}\right) + \left(k x_2 - \omega t + \phi_{0,2}\right)}{2} \\ &= A \cos \frac{k \left(x_1 - x_2\right) + \left(\phi_{0,1} - \phi_{0,2}\right)}{2} \\ & \quad \: \times \sin \left(k \frac{x_1 + x_2}{2} - \omega t + \frac{\phi_{0,1} + \phi_{0,2}}{2}\right) \\ &= A \cos \frac{k \Delta x + \Delta \phi_0}{2} \sin \left(k_2 \langle x \rangle - \omega t + \langle \phi_0 \rangle\right) \end{align*}$
The condition for destructive interference is:

$\begin{align*} D (x, t) &= 0 \\ 2 a \cos \frac{\Delta \phi}{2} \sin \langle \phi \rangle &= 0 \\ \cos \frac{\Delta \phi}{2} &= 0 \\ \frac{\Delta \phi}{2} &= n \pi + \frac{\pi}{2}, \quad n \in \mathbb{Z} \\ \Delta \phi &= \left(2 n + 1 \right)\pi \end{align*}$
Similarly, for constructive interference:

$\begin{align*} D (x, t) &= \pm 2a \sin \langle \phi \rangle \\ 2 a \cos \frac{\Delta \phi}{2} \sin \langle \phi \rangle &= \pm 2a \sin \langle \phi \rangle \\ \cos \frac{\Delta \phi}{2} &= \pm 1 \\ \frac{\Delta \phi}{2} &= n \pi, \quad n \in \mathbb{Z} \\ \Delta \phi &= 2 n\pi \end{align*}$
Beat interference:
x₁ = x = x₂
k₁ = k = k₂
ϕ_0,1 = ϕ₀ = ϕ_0,2
The displacement equation is:

$\begin{align*} & \quad \: D (x, t) \\ &= A \cos \frac{\Delta \phi}{2} \sin \langle \phi \rangle \\ &= A \cos \frac{\left(k x - \omega_1 t + \phi_0\right) - \left(k x - \omega_2 t + \phi_0\right)}{2} \\ & \quad \: \times \sin \frac{\left(k x - \omega_1 t + \phi_0\right) + \left(k x - \omega_2 t + \phi_0\right)}{2} \\ &= A \cos \left(\frac{\Delta \omega}{2} t\right) \sin \left(k x - \langle\omega\rangle t + \phi_0\right) \\ &= A \cos \omega_\text{mod} t \sin \left(k x - \langle\omega\rangle t + \phi_0\right) \end{align*}$
where ω_mod is the modulation angular frequency.

Now, at last, I can show the general equations (approximations) for multiple-slit interference.

For single-slits, the minima angles (θ_min) can be found using:
a sin θ_min = n λ where n is a positive integer and a is the slit width.

For double-slits, multiple-slits, or diffraction gratings, the primary maxima angles can be found using:
d sin θ_max = n λ where n is zero or a positive integer and d is the slit separation.

The two results actually apply to any slit interference: θ_min will be the primary minima angles and θ_max will be the primary maxima angles. The shape of the intensity plot will always have an envelope curve that is shaped like the intensity curve of a single-slit diffraction, and the minima of this envelope curve is the "primary minima".

So that was my brief summary of light interference.

Saturday, November 1, 2008

reciprocal distribution

Here's another post on math stuff, although it's also implicitly related to computing.

Suppose I generate a random real number X between 0 and 1 with a uniform probability distribution. Let Y = 1/X. What is the probability distribution of Y? What is the expected value of Y?

To solve this question (which I asked myself), I would start by finding the probability distribution of X.

For me, the main problem is that this is a continuous probability distribution rather than a discrete one. Therefore, it makes no sense to ask the probability of getting, say, exactly 0.5. In discrete probability distributions, we can describe a function called "Pr" or the probability function. Let's take a dice distribution as an example. A dice can have 6 possible values: D = {1, 2, 3, 4, 5, 6}. Each value has an equal probability of 1/6. We can define a probability function for it easily:

Pr(D = d) = 1/6 (this is read as: "the probability of random variable D getting a specific value of d is 1/6")

This is a discrete uniform distribution, because the probability is the same for every d. The equivalent way of saying this is:
$\Pr (D = d) = \begin{cases} d = 1: & \frac{1}{6} \\ d = 2: & \frac{1}{6} \\ d = 3: & \frac{1}{6} \\ d = 4: & \frac{1}{6} \\ d = 5: & \frac{1}{6} \\ d = 6: & \frac{1}{6} \end{cases}$
While this form is rather redundant, it does explicitly show that we can define a specific probability for every value of d. For a continuous distribution, it no longer makes sense to assign a probability to every possible value, since there is now an infinite number of possible values (continuous).

However, it still makes sense to define a probability for every "chunk" of values. Consider the random real number X between 0 and 1. It is absurd to ask for the probability of getting X = 0.5, but it still makes sense to ask for the probability of 0.45 ≤ X < 0.55. This is a range of possible X values from 0.45 to 0.55. The width (Δx) of this range is 0.55 - 0.45 = 0.1. We can make use of something called probability density function, f(x), as follows:
$\Pr(0.4 \le X < 0.6) = \int_{0.4}^{0.6} f(x) \operatorname{d}\! x$
This is in many ways similar to the frequency density used in histograms.

Since the probability is defined by the integral of the probability density function (i.e. area under the f(x)), the total area under the probability density graph must be 1 in order for the total probability to be 1. This is done by normalizing the function using a constant multiple.

For example, a continuous uniform distribution between x = 0 and x = 1is rather simple:
$f(x) = \begin{cases} 1, & x \in \left[0, 1\right] \\ 0, & x \notin \left[0, 1\right] \end{cases}$
It's easy to see that this graph is already normalized since the area is just a 1 by 1 square. This is just the distribution of the variable X in the question posed in the beginning of this post. Now that we know the probability density of X, how can we find the probability density of Y?

Suppose I pick 2 values a and b (and a < b) What is the probability of X between these two values? It is simply Pr(a ≤ X < b). By intuition, it should make sense that the probability of Y between 1/b and 1/a (I swapped the order because 1/b < 1/a) should be equal to the probability of X between a and b, since x = a implies y = 1/a and x = b implies y = 1/b. Therefore, we have the equation:
$\Pr\left(a \le X \le b \right) = \Pr\left(\frac{1}{b} \le Y \le \frac{1}{a} \right)$
Using the integral form:
$\int_a^b f(x)\operatorname{d}\!x = \int_{1/b}^{1/a} g(y)\operatorname{d}\!y$
where f(x) and g(y) are the probability densities of X and Y respectively. Differentiate both sides with respect to a and simplify the expression:
$\begin{align} \frac{\operatorname{d}}{\operatorname{d}\!a} \int_a^b f(x)\operatorname{d}\!x &= \frac{\operatorname{d}}{\operatorname{d}\!a} \int_{1/b}^{1/a} g(y)\operatorname{d}\!y \\ - \frac{\operatorname{d}}{\operatorname{d}\!a} \int_b^a f(x)\operatorname{d}\!x &= \frac{\operatorname{d}\!\left(1/a\right)}{\operatorname{d}\!a} \left[\frac{\operatorname{d}}{\operatorname{d}\!\left(1/a\right)} \int_{1/b}^{1/a} g(y)\operatorname{d}\!y \right] \\ - f(a) &= -\frac{1}{a^2}g\left(\frac{1}{a}\right) \\ a^2 f(a) &= g\left(\frac{1}{a}\right) \\ \frac{1}{y^2} f\left(\frac{1}{y}\right) &= g\left(y\right) \end{align}$
Explanation:
Step (2) left side: swap the limits of integral, thus swapping the sign.
Step (2) right side: use chain rule to get the differentiation with respect to 1/a instead of a.
Step (3) left side: apply fundamental theorem of calculus.
Step (3) right side: differentiate 1/a and apply the fundamental theorem of calculus to the integral part.
Step (4): multiply both sides by -a².
Step (5): substitute a = 1/y so that g is in terms of y (the y here does not have much to do with the y in (1) and (2)).

Now we have an expression of g(y), in terms of f(x) (which is the same as f(1/y), since x = 1/y). Since f(x) is a uniform distribution with f(x) = 1 for all values of x within [0, 1], we can define g(y) for all values of y within [1/1, 1/0), or [1, ∞), to be 1/y². Because f(x) is zero elsewhere, g(y) is zero elsewhere as well. So this is the formula of g(y):
$g\left(y\right) &= \begin{cases} \frac{1}{y^2}, & x \in \left[1, \infty\right) \\ 0, & x \notin \left[1, \infty\right) \end{cases}$
While not obvious from the equation, this function is definitely normalized, and can be easily proven to be so using improper integrals.

So, the aim was to find an expected value for the distribution. What is the formula for a continuous distribution? The expected value for a discrete random variable X is defined by:
$\operatorname{E}(X) = \sum_i x_i \Pr\left(X = x_i\right)$
Therefore, by analogy, the expected value for a continuous random variable (with probability density f(x))should be:
$\operatorname{E}(X) = \int_{-\infty}^{\infty} x f(x) \operatorname{d}\!x$
So the next step would be to evaluate this integral for Y:
$\begin{align} \operatorname{E}(X) &= \int_{-\infty}^{\infty} y\, g(y) \operatorname{d}\!y \\ &= \int_{-\infty}^1 y \cdot 0 \operatorname{d}\!y + \int_1^{\infty} y \cdot \frac{1}{y^2} \operatorname{d}\!y \\ &= 0 + \int_1^{\infty} \frac{\operatorname{d}\!y}{y} \\ &= \infty \end{align}$
Explanation:
Step (2): Split the integral into 2 parts, one for x ≥ 1 and one for x < 1. This is done because g(x) is defined piecewise.
Step (3): Integral of zero is zero. Simplification.
Step (4): Integral of 1/y is divergent because the natural logarithmic function is unbounded.

It looks like the expected value is infinite. What does that mean? Basically, if one were to simulate this on a real computer program (tried that myself) to calculate expected value, one could get utterly random values stretching from very small to very large.

The expected value of g(y) is infinite - this is a boring conclusion, isn't it? Well, over the steps listed here, I have found a technique that allows me to derive the probability density of one variable based on a formula that links two continuous variables, with the second variable given. So it was not a complete waste of time. This is the general solution for a variable Y linked to X by a formula x = u(y):
$g\!\left(y\right) = f\!\left[u\!\left(y\right)\right] \left|\frac{\operatorname{d}\!u}{\operatorname{d}\!y}\right|$
where f(x) is the probability density of X and g(y) is the probability density of Y. The formula does not work if u(x) is multivalued. Still, I think this formula should work under ordinary conditions. If there's any mistakes in this post, feel free to leave a comment.