Re: Please check my article on DSP basics

On Mar 18, 2:32 am, "Humanise" <k...@xxxxxxxxxxxx> wrote:
Below is my attempt at describing the needed maths for programming the DSP
part of the Humanise.org project. (Seehttp://humanise.org).  Humanise.org
is an open source project to create a high quality software oriented metal
detector suitable for humanitarian demining.  I'm trying to do this article
with simple school type maths (no integration, no radians, no complex maths
etc).

In particular I will be trying to describe how a Phase Lock Loop (PLL) and
how a Two Channel Lock-In Amplifier work.  The actual PLL and Lock-In
Amplifier is described in another article.  This is just the maths behind
them.  I expect it could also be used as an introduction to DFT and then
FFT.  I looked for something as simply as this that still gave exact maths
sufficient for it to be programmed from but although I found many attempts
by others, I never liked any.  By doing my own version, it means that I can
directly host it instead of just linking to it and then finding later that
what I had linked to is no longer there.

I don't consider myself good at this type of maths so I'm looking for
others to check out the maths.  When I convert it to HTML I will convert a
lot of the symbols to their proper symbols and get rid of all the * used as
multiplication signs.  I also hope to illustrate it with sine wave graphs
etc.  I'm not good at making diagrams so I'm hoping that others wanting to
help in the Humanise.org will make some diagrams for it.

Can you please read it over and tell me any mistakes that I have made.

---------------------------------------------------------------------------------------------------
Note on AC/DC terminology:

With Electronics we talk about currents being DC (Direct Current) or AC
(Alternating Current) or a mixture of both.  When converted to numbers and
graphed on a typical X and Y axis graph, the DC is represented as a shift
of the graph on the Y axis while the AC is a continuing alternation on the
Y axis without any overall shift.

Throughout the discussion and mathematics on this website we will continue
to use the term DC and AC for similar actions on streams of numbers
representing waveforms even though there is no actual current still
involved.  Similarly we assume that the typical X and Y axis graphs exist
for each waveform with the Y axis showing the wave's Amplitude and the X
axis being used for Time.

Gentle Maths Version:

Now lets try describing the maths behind the Phase Lock Loop and the Two
Channel Lock-In Amplifier using nothing more than simple high school
maths.

It has been stated that the heart of all Fourier analysis is, amazingly,
all based on the following single high school trigonometry formula for the
product of two sines:

sin(X1) * sin(X2) = 1/2 * cos(X1-X2) - 1/2 * cos(X1+X2)

For sine waveforms, the above X1 and X2 are based on time and frequency.

For simplicity, let us assume that Frequency is given to us in units of
degrees per second where a complete cycle of the sine wave is 360 degrees..
This is sometimes called Angular Velocity.  We will use F to represent this
Frequency.  Let us assume that the sine table or sine lookup we use (sin
for short) is based on degrees.  Let us also assume that time represented
by t is given in seconds.  For these sine waves, the X1 and X2 above are as
follows:

   (F * t)  + Initial_Phase_Angle_In Degrees

As F is an Angle per time (degrees per second), at each time point we
multiply the Angle per time by the current time (in seconds) we get an
angle (in degrees) which is what we need for looking up the the sin
lookup.

Unfortunately, frequencies are not normally specified as degrees per
second.  Rather, they are specified in cycles per second.   This creates a
mishmash between the units used for angle with the the sin lookup and the
units used for angle per time in the frequency specification.  Let us
represent these cycles per second frequency by f.  Then

    F = 360 * f

With the sin lookup still being based on angles that have 360 units per
waveform unit, X1 and X2 now needs to be as follows:

    (360 * f * t) + Initial_Phase_Angle_In_Degrees

We have had to place in the factor of 360 to take account of the mishmash
with 360 degrees in each cycle in the sin lookup but frequency being
specified with whole cycles per second.

Note that with each of the above definitions of X1 and X2, once we have a
fixed input frequency and phase, the only thing that is varying is time.  
For frequencies, the initial equation relates the Amplitude of these waves
to time.

Now lets start using the initial equation with waveforms.

sin(X1) * sin(X2) = 1/2 * cos(X1-X2) - 1/2 * cos(X1+X2)

When this equation is used on waveforms, it tells us what happens when we
have one sin(X1) waveform input plus another another sin(X2) waveform input
and we multiply them together as they are being input.  We are just
multiplying the same point in time on one waveform to the same point in
time of the other waveform.  The resultant waveform is (1/2 * cos(X1-X2)) -
(1/2 * cos(X1+X2)).

Lets see what happens when X1 and X2 are different frequencies.  Lets
assume that these frequencies are f1 and f2 and that both of these are of
an amplitude of 1.  If we look at the first resultant term

 1/2 * cos(X1-X2)

With waveform X1 being a frequency f1 and X2 being a frequency f2, this
will be equal to

 1/2 * cos((360 * f1 * t) - (360 *f2 * t))

This is simplified to

 1/2 * cos(360 * (f1-f2) * t)

This is a waveform with a frequency that is the difference between the two
input frequencies. As we know these frequencies are different (we made that
assumption at the beginning) then this must produce a wave form that has a
non zero frequency that is the difference between the two frequencies.  As
this is a standard sine waveform, it is symmetric about the X axis.  This
is, its an AC only waveform with an amplitude of one half (due to the 1/2
previous to the cos(X1-X2)).

Similarly with the 1/2 * cos(X1+X2) term, the result

 1/2 * cos((360 * (f1+f2) * t) + 2 * Initial_Phase_Angle)

  is a waveform that is the addition of the two frequencies.  This
produces a second waveform of half the amplitude of the input waveforms.
In this case it has a frequency that is the sum of the two frequencies.  As
this is also a standard sine waveform, it is also symmetric about the X
axis.

The addition of two pure AC complete waveforms always produces a resultant
waveform that is pure AC only.  The word 'complete' is here used to
indicate that we are looking at complete cycles, not just part of a cycle
of the waveform.  As can be seen from the above, so long as the frequencies
are different, there is no DC offset in the result.

Lets assume that X1 and X2 above are exactly equal.  Start again with the
initial equation:

sin(X1) * sin(X1) = 1/2 * cos(X1-X1) - 1/2 * cos(X1+X1)

As Cos(X1-X1) = Cos(0) = 1

Sin(X1) * Sin(X1) = 1/2 - (1/2 * cos(2 * X1))

In this case, the same input sine waves multiplied by each other produce a
DC offset (equal to half the amplitude of the input sine waves) plus a wave
at double the frequency i.e. (2 * X1) and half the amplitude of the input
sine waves.

One strange aspect of sine waves is that whether a wave is added or
subtracted doesn't change its position with respect to the Y axis.  That
is, it doesn't give or change the wave's DC offset.  As the wave is minus,
it is equivalent to the wave being 180 degrees out of phase.  That is, the
following equation is equivalent to the last one:

Sin(X1) * Sin(X1) = 1/2 + (1/2 * cos((2 * X1) + 180 degrees))

Lets assume that X1 and X2 above are equal except that X2 is phase shifted
forward by 90 degrees.

Sin(X1) * Sin(X1 + 90 degrees) = 1/2 * cos(0 - 90 degrees) - 1/2 * cos((2
* X1) + 90 degrees)

As Cos(-90 degrees) = 0  

Sin(X1) * Sin(X1 + 90 degrees) = - 1/2 * cos((2 * X1) + 90 degrees)

                       =   1/2 * cos((2 * X1) - 90 degrees)

In the case of waveforms this means that when input sine waves of the same
frequency but with a phase offset of 90 degrees are multiplied by each
other, a wave at double the frequency i.e. (2 * X1) and half the amplitude
of the input sine waves is produced.  As the wave is minus, it is
equivalent to the wave being 180 degrees out of phase but the extra 90
degrees in ((2 * X1) + 90 degrees) then shifts it back 90 degrees.  Note
that this wave is symmetrical about the X axis.  There is no DC offset.

If we did the same as above but with X2 phase shifted back by 90 degrees
(the same as phase shifted forward by 270 degrees) then we would get a
similar result in that the resultant wave has no DC offset.

Sin(X1) * Sin(X1 - 90 degrees) = - 1/2 * cos((2 * X1) - 90 degrees)

Lets assume that X1 and X2 above are equal except that X2 is phase shifted
by 180 degrees.  This is the equivalent to being inverted.   Starting at
our initial equation we see that:

sin(X1) * (-1 * sin(X1)) = -1 * ( 1/2 * cos(X1-X1) - 1/2 * cos(X1+X1))
                     =  1/2 * cos(X1+X1) - 1/2 * cos(X1-X1)

 Again, as cos(0) = 1

                     =  (1/2 * cos(2*X1)) - 1/2

In the case of waveforms this means that when input sine waves that are
equal except one is inverted from the other, are multiplied by each other,
they produce a minus DC offset (equal to half the amplitude of the input
sine waves) plus a wave at double the frequency i.e. (2 * X1) and half the
amplitude of the input sine waves.

From these equations we can see that we only obtain a DC offset when we
multiply two waveforms if they are exactly at the same frequency.  Even
then, we only obtain the DC offset if the waves are not at 90 or 270
degrees phase difference from each other.  The DC offset reaches a maximum
positive value when the waves are exactly in phase with each other and
reach a maximum negative value when the waves are 180 degrees out of
phase.

Also note that we have simplified the basic equation by leaving off the
amplitudes.  We have assumed that the amplitudes of the two input waves are
both 1.  Assume the amplitude of the sin(X1) input wave is A and the
amplitude of the sin(X2) input wave is B then the original equation would
be

(A * sin(X1)) * (B * sin(X2)) = (A * B ) * (1/2 * cos(X1-X2) - 1/2 *
cos(X1+X2))

If the amplitudes of the input sin(X1) and sin(X2) waves is not 1 then all
the previous formulas can be corrected by simply multiply the results by a
factor of (A * B).

The previous equations should have explained the workings of the Phase
Lock Loop.  The following expands this to show the workings of the Two
Channel Lock-In Amplifier.

Lets see what happens if we create an X2 waveform that is the same as an
input X1 except that it is phase shifted an angle Theta from X1.

sin(X1) * sin(X1 + Theta degrees) = (1/2 * cos(X1-X1 - Theta degrees)) -
(1/2 * cos(2*X1 + Theta degrees))
                 =  (1/2 * cos(- Theta degrees)) - (1/2 * cos(2*X1 + Theta
degrees))

The second term (1/2 * cos(2*X1 + Theta degrees) by itself is a pure sine
wave that is symmetrical about the X axis.  It has no DC component.  If we
are only interested in the resultant DC component then we can throw away
this term.  The first term gives us the value of the resultant waveform if
we averaged it, filtering out the AC component.  The DC component would be

1/2 * cos(- Theta degrees) = 1/2 * cos(Theta degrees)

cos(Theta degrees) is a constant so we have a constant DC offset of 1/2 *
cos(Theta degrees).  If we perform this with an input wave amplitude of A
and an internally created wave of amplitude B, and measure the resultant DC
offset as Y1 then

   Y1 = 1/2 * A * B * cos(Theta degrees)

If we now create an X3 waveform that is the same as X2 except that it is
phase sifted by (Theta plus 90) degrees from X1 and use that to multiply
against another copy of the input X1 waveform.

sin(X1) * sin(X1 + (Theta + 90) degrees) = (1/2 * cos(X1-X1 - (Theta + 90)
degrees)) - (1/2 * cos(2 *X1 + (Theta + 90) degrees))

Again we are only interested in the DC component so we throw away the
second term.  The first term gives us

     1/2 * cos((Theta + 90) degrees) = 1/2 sin(Theta degrees)

If we perform this with an input wave amplitude of A and an internally
created wave of amplitude of B, (assume that X3 is the same amplitude as
X2) and measure the resultant DC offset as Y2 then

  Y2 = 1/2 * A * B * sin(Theta degrees).

If we now square Y1 and Y2, add these squares together and then take the
square root

(Square root of (Y1**2 + Y2**2))
 = square root of  ((1/2 * A * B * cos(Theta degrees))**2 + (1/2 * A * B *
sin(Theta degrees))**2)
 = square root of  (1/4 *A*A*B*B* (cos(Theta degrees)**2)  +  1/4
*A*A*B*B* sin(Theta degrees)**2)
 = square root of  (1/4 *A*A*B*B* (cos(Theta degrees)**2)  +  1/4
*A*A*B*B* sin(Theta degrees)**2)
 = square root of  (1/4 *A*A*B*B * ((sin(Theta degrees)**2) + (cos(Theta
degrees)**2)))

     Note: sin(X)**2 + cos(X)**2 = 1
     This is a trigonometrical identity that is easily derived from the
main trigonometrical identity we have been using.  Using this reduces the
above to

 = square root of  (1/4*A*A*B*B*)
 = 1/2 * A * B

Consequently, if we know B then the Amplitude A can be calculated from

  A = 2 * (Square root (Y1**2 + Y2**2)) / B

Using the above, if we have an input waveform X1 that contained a known
frequency but we didn't know the phase or amplitude of this frequency, we
can find the amplitude by creating two reference waveforms, both of
Amplitude B and frequency the same as X1.  The second internally created
reference waveform would be created with a 90 degree phase shift from the
first reference waveform.  We also create a second copy of the input
waveform.

We now multiply the input waveform by the first reference waveform and
average the result so that we measure the resultant Y1.  We also multiply
the copy of the input waveform by the second reference waveform and average
the result so that we measure the resultant Y2.  The amplitude A of the
known frequency within the input waveform can now be calculated as

  A = 2 * (Square root (Y1**2 + Y2**2)) /  B

Lets denote by Y3 what would have been measured had the internally created
reference waveform been exactly in phase with the input frequency.  This
gets rid of the 1/2 that's needed when comparing to the input waveform.
The measured Y1 and Y2 are the sine and cosine components of the Y3
amplitude.  They have the same relationship as the clock hand (of Y3
length) has to the X and Y axis or the same relationship as the hypotenuse
(of Y3 length) of a right angle triangle.

Consequently, we can find the true amplitude or "magnitude" as the
hypotenuse of this right triangle by taking the square root of the sum of
the squares of the sine and cosine components. This is called the "vector
sum" of the components.  Similarly, we can find the phase as the angle
whose tangent is equal to the sine component divided by the cosine
component.

The sine and cosine components are often referred to as the "in-phase" and
"quadrature" components. Some engineering and mathematicians prefer to use
"real" and "imaginary" for the sine and cosine components.   From this they
create some very intimidating expressions using an imaginary number (square
root of -1) in a style of maths called "complex notation".

Also note that the first internally created reference frequency can be at
any phase difference to the actual signal.  Given any particular pure
sinewave, without a DC component, at a particular frequency f1, it can be
fully described by specifying its amplitude and phase.  Alternatively, the
sine wave can be described by specifying two amplitudes, these being the
equivalent in-phase and quadrature components of that frequency.  As we can
change where (which phase angle) our first reference frequency is, there is
effectively an infinite number of pairs of amplitudes that can be used as
an equivalent to the phase plus amplitude pair.  As long as we keep in mind
where the the reference is in each case, all these pairs of amplitudes are
equivalents.  

With VLF metal detectors, where we are both transmitting the frequency on
one coil and receiving the resultant frequency on another coil, the sine
wave being output is normally used as our reference (sine component)
frequency for the primary received frequency.  Consequently, in this case
it clearly makes sense to call the sine wave component the in-phase
component.  Changes in the received amplitude of this component are
normally caused by resistive property effects in the soil and anomalies
near the coil.   Note that resistive is here used to indicated properties
that change the received voltage in phase with the signal like resistance
does.  It doesn't necessarily mean resistance itself. Changes in the
received amplitude of the cosine component are normally caused by reactive
(such as capacitance and inductance) property effects of the soil and
anomalies near the coil.

When programming using any of the above equations, you need to recognise
the order in which the sine or cosine waves will be processed.  The sine
graphs with time increasing on the X axis are a representation of waves
that are travelling from right to left.  Generally, we expect things on
paper to be represented from left to right so this works out to be opposite
to what many people expect.

To confirm that it is this direction note what happens if we place the
sine or cosine graph onto a long tape with time in seconds marked as normal
along the X axis.  We cover over the tape except that we leave a small view
port that allows us to see the tape at the zero time mark (normally the Y
axis).  We now fire a starters gun to mark our zero time point.  Following
the gun firing we pull the tape so that the time marks of 1 second, 2
seconds etc are viewable at our view port at 1 second, 2 seconds etc after
our zero time.  To do this, we have to pull the tape to the left so that
the graphed sine wave is moving right to left.

The conclusion here is that if you have a sine wave and a cosine wave
being input to the processor, the cosine wave will reach the processor at
1/4 wave cycle ahead of the sine wave or 3/4 wave cycle after the sine
wave.  The equation cos(X) = sin(X + 90 degrees) indicates that cosine is
before sine by 90 degrees.

We should also point out the usefulness of the earlier discussed fact that
multiplying two frequencies produces a sine wave that is the difference of
the two frequencies plus a sine wave that is the sum of the two
frequencies.

As an example, Let us assume that we are interested in using an audio card
input to look at an input that has a frequency of 300Khz.  This is a higher
frequency than current audio cards can handle correctly.  To make it
suitable we could create an oscillator that is for example a fixed 295khz..
We then multiply this oscillator frequency to the 300Khz input.  The result
is a waveform at (300 + 295)Khz plus a wave form at (300 - 295)Khz.  We can
then use a low pass filter to filter out the 595Khz waveform leaving us
with just a 5Khz waveform which is very easy for the audio card to handle..
As the original 300Khz either changes in frequency or changes in amplitude,
the 5Khz input to the audio card will similarly change.  This principle is
used in many applications such as the super-heterodyne receiver

It also explains a common problem in many wave systems called
Intermodulation Distortion. This distortion occurs when two sine waves of
frequencies X1 and X2 in the input result in the creation of other unwanted
frequency components.  These frequencies include (X1+X2), (X1-X2),
(2*X1-x2), (2*X2-X1), and generally (mX1 ± nX2) for integer m and n.
Generally, the size of the unwanted frequency components falls rapidly as m
and n increase.

A similar explanation for some of the above is given at the following
link:http://www.daqarta.com/eex01.htm

That's not maths...that's kids stiff from school. For DSP you will
need the Fourier transform,Fourier Series,Convolution integral and all
the discrete equivalents.
z-transform of course,power series expansions,Laplace transforms etc
etc
Not for amateurs.

Hardy
.

Loading