Spectrum by Similarity

Craig Stuart Sapp <craig@ccrma.stanford.edu>
Peabody Conservatory
14 November 2002

This Mathematica Notebook demonstrates the basic principle used to calculate a spectrum from a signal.  A spectrum is a plot of sinusoid frequencies on the x-axis versus the amplitude of that sinusoid displayed on the y-axis.  This notebook demonstrates how the amplitude of a sinusoid is calculated from any arbitrary signal.  The basic principle is this: (1) multiply the signal by a sinusoid at the desired frequency along the x-axis, and (2) sum the discrete elements in this multiplied signal together to determine the amplitude of the sinusoid present in the original signal.  This process is checking to see how similar the spectrum is to any given test sinusoid.

Part I: The Real World

First we define some variables to create a test signal.  The Amp variable stores the list of amplitudes for each sinusoid component in the signal.  Freq contains a list of the frequencies in the signal.  Periods indicates how many fundamental periods (of 100 Hz in this case) are to be displayed in the plots and used in the calculations.  Srate is the sampling rate of the signal.

The next cell defines a signal with the variable t representing time in seconds.

Here is a plot of signal[t] where the horizontal axis is in milliseconds:

The previous picture is the continuous picture of the signal, but if a computer is processing this sound, it has to do it in discrete steps.  The following command samples the signal at the sampling period (which is 1 millisecond in this case).

The following  plot shows the discrete samples just calculated in the above cell:

The following plot demonstrates that the sampled values really are derived from the continuous signal.  Note that the samples line up with the original continuous function.

Now that samples from the test signal have been generated, let's make a table of the individual sinusoids which we know are present in the signal.  The following cell makes a list of the samples for each sinusoid present in the signal.

Here is a plot of the individual frequencies.  Note that all of the amplitudes are one.  Our job is to figure out what were the original amplitudes of each sinusoid in the original signal.

In order to determine the amplitudes of the sinusoids in the original signal, you must do two steps (1) mulitply the signal by each test sinusoid, and (2) add up all of the elements in the resulting signal and then divied by the normalization factor.  These sets of number will be the original amplitudes of the sinusoids that we started out adding to create the original signal.
Below is the first step which involves multiplying the original signal by each of the test sinusoids:

Here is a plot of the data generated above.  These plots do not have any physical meaning, but are a step on the way to calculating the amplitudes of the sinewaves in the original signal:

Next, each of these multiplied signals need to be added, sample by sample:

Finally, normalize the values:

Now compare this result to the amplitudes of the sinusoids present in the original signal:

Part II: The Complex World

You saw in the previous section that the amplitudes of the original sinusoids can be extracted by multiplying the spectrum by the individual sinusoids that we want to find the amplitude of in the signal.  In practice, we will need a more comprehensive sinusoid, called the complex sinusoid.  This sinusoid can handle any phase shift of the real sinusoids found in the signal and will be demonstracted in the next section.
A complex sinusoid is a combination of a sinewave and a cosinewave which are at right angles to each other.  Here is a plot of a complex sinusoid:

The horizontal axis represent the real numbers, while the vertical axis represents the Imaginary numbers.  The origin of the complex sinusoid is on the left, and you can see that the green projection of the sinusoid creates a cosine.  The blue projection of the complex sinusoid (onto the imaginary axis) creates a sinewave.  Thus the mathematical equation for the complex sinusoid is:    csin(x) = cos(x) + i sin(x), where i is the square root of negative one, and csin is the complex sinusoid.
The complex sinusoid itself can be represented as an imaginary power of the number e (approximately equal to 2.71828): e^(i x).  Here is a demonstration that there are two equal ways of calculating the value of a complex sinusoid (Euler's Identity):

Since the difference between ans1 and ans2 is zero, the methods are equal.  Try it for any other value of x to verify that the calculations are still equal.

Now let us use the complex sinusoids to measure the amplitudes of the original sinusoids used to make the signal.  Here is the original signal:

Analygous to part I, create a set of test sinusoids which will be used to examine the signal:

Finally, normalize:

Now compare this to the amplitudes of the sinusoids present in the original signal:

These amplitudes also exactly match the amplitudes found by using only real sinusoids.

Part III: Phase and why the complex world is better than the real world

The complex sinusoid is much more accurate at finding the original sinusoid amplitudes.  This section will demonstrate this fact.  A problem occurs when you try to find the amplitude of a sinusoid component in a signal which has an unknown phase.  In part I, the phase was exactly 0 (or 90 degrees if you are thinking of sinewaves).  But what happens if the phase of each sinusoid is random?  Real sinusoids cannot extract the full amplitude of the original sinusoid in the signal, and this amplitude will more than likely be underreported.
Lets test this theory by creating a signal which contains the same sineusoids as in the previous signal which also have the same amplitudes, but have arbitrary phases:

The next cell defines a signal with the variable t representing time in seconds.

Here is a plot of signal[t] where the horizontal axis is in milliseconds:

Notice that the above signal with different phases for each sinusoid  looks different from the original signal:

The previous picture is the continuous picture of the signal, but if a computer is processing this sound, it has to do it in discrete steps.  The following command samples the signal at the sampling period (which is 1 millisecond in this case).

The following  plot shows the discrete samples just calculated in the above cell:

The following plot demonstrates that the sampled values really are derived from the continuous signal.  Note that the samples line up with the original continuous function.

Finally, normalize (divide) by the number of samples per fundamental frequency (10 samples from 100 Hz):

Now compare this to the amplitudes of the sinusoids present in the original signal:

Notice that the amplitudes this time are not equal to the amplitudes of the sinusoids.  Depending on the phase of the sinusoids, the measured amplitudes will range from the full amplitude to the negative full amplitude.

Measuring the correct original amplitude

The only way to extract the original amplitudes from the signal is to use complex sinusoids which resist the effects of a change in phase.  This section demonstrates that you can measure the correct original amplitudes of the sinusoids from the original signal.

The amplitude is taken as the magnitude of the complex number, which is also the absolute value of the complex number:

Finally, normalize (divide) by the number of samples per fundamental frequency (10 samples from 100 Hz):

Now compare this to the amplitudes of the sinusoids present in the original signal:

These amplitudes also exactly match the amplitudes of the arbitrary-phased signal, while using the real test sinusoids did not give the correct amplitudes.

That is not all -- you can also determine the phase of the orignal sinusoids in the signal:

Here is the original phase, which matches exactly:

Part IV: Symmetry

Even and Odd functions

The problem with using only a real sinsoid is that you have to choose a sinusoid which will not correctly identify the amplitude of another sinusoid with the same frequency.  This section demonstrates this fact.

If you multiply the sine and the cosine functions together, you will get a sinewave with twice the frequency of the original:

This implies that there is no Cosine quality in a Sine and no Sine in a Cosine.  However, if you multiply a sine by a sine, all of the function will be positive:

Thus, a sinewave cannot "see" a cosinewave and a cosinewave cannot "see" a sinewave.  This is related to the symmetry of each function.  A cosine was has even symmetry because it looks the same on each side of the line x=0 if you place a mirror on that line.  A sine function has odd symmetry because it can be rotated 180 from the point (x,y) = (0,0) and still look the same.   All functions can be broken down into two pieces: one which has even symmetry, and one which has odd symmetry.

An arbitrary phased sinusoid which is neither a sinewave or a cosine wave can be broken down into two pieces, one which is symmetric and one which is antisymmetric.  The symmetric piece corresponds to a cosine and the antisymmetric piece corresponds to a sine.  The amplitudes of the cosine and sine add up to the amplitude of the original arbitrary-phased sinusoid.

For example, try a phase of Pi/4 which is half-way between a cosine with 0 phase and a sine with 0 phase (or a cosein with Pi/2 phase):

This sinusoid can be broken down into two sinusoids: one sinewave and one cosine which can be added together to create the original phased sinusoid:

Check to see that the functions are equal by plotting their difference:

Here is a function which will extract the sine and cosine of any sinusoid:

With one sine and one cosine, you can extract the amplitude of a sinusoid with any phase in a signal.  Remember that the complex sinusoid is a cosine on the real axis and a sine on the imaginary axis, so you now know why the complex sinusoid can extract the amplitude of a sinusoid which has any phase.

V Spectrum plots

The spectrum is usually given as a plot of amplitude v frequency of sinusoids, but there is also the phase plot which is the other part of the spectrum.

Here are functions to display the sinusoid frequencies, amplitudes and phases all in one spectral plot:

It is interesting to view the continuous version of the spectrum in 3D: