Math 236  •  9 May 2001

Diagnosis and cure of hum

Sinuoidal hum function

Clear[ω,  λ,  δ,  A,  h,  s,  f]

The typical sinusoidal function of frequency ω has, for suitable parameters A and δ, the form:

h[t_]  :=  
      A Sin[2 π (ω t  +  δ)]

For a 60-cycle hum (which will represent the "noise" added to an input signal), we use parameters:

\(\(\(ω\  = \ 60\)\(;\)\(\ \ \ \ \ \ \ \ \ \)\( (*\ 
      frequency\ \ *) \)\(
\)\(λ\  = \ 
      1/ω\)\(;\)\(\ \ \ \ \)\( (*\ 
      period\ \ \ \ \ *) \)\(
\)\(δ\  = \  .05\)\(;\)\(\ \ \ \ \ \ \)\( (*\ phase\ \ \ \ \ \ *) \)\(
\)\(A\  = \ 
      1.5\)\(;\)\(\ \ \ \ \ \ \ \)\( (*\ amplitude\ \ *) \)\)\)

One cycle:

Plot[h[t],  {t,  0,  λ}];

Several cycles:

Plot[h[t],  {t,  0,  5  λ}];

Listen to the hum:

Play[h[t],  {t,  0,  100  λ}];

Significance of phase δ:

h[t]  /.  t  ->  -δ/ω

Plot[h[t],  {t,  -δ/ω,  λ}];

Signal with hum

In this first example, we use the following input signal (which we shall soon treat as if we did not know it):

s[t_]  :=  
      0.01  Sin[120 π t]  +  0.02  Cos[120 π t]  +  
        4 (t  -  0.5)^2  - t

What we will then "observe" or "measure" is the following output signal, in which the signal has been corrupted by the addition of the hum as "noise":

f[t_]  :=  s[t]  +  h[t]

Compare plots of the input (non-noisy signal) and output (noisy—signal plus noise) functions:

Plot[s[t],  {t,  0,  1}];

Plot[f[t],  {t,  0,  1}];

Compare sounds of the input and output  functions:

Play[s[t],  {t,  0,  100  λ}];

Play[f[t],  {t,  0,  100  λ}];

The problem

We are given only the output signal f and the hum frequency ω = 60.  The overall problem is to determine what the input signal s was.

The solution splits into two parts:

Part 1: Diagnosis

Find the function h(t).

Part 2: Cure

Reconstruct s(t)  =  f(t)  -  h(t).

Representation of hum as linear combination of basic functions

First we change the form of the hum ("noise") function h in order to express it in terms of more basic functions.  We shall do this in general, without any particular values of the parameters.

Clear[ω,  λ,  δ,  A]

h[t]  /.  a_ (b_  +  c_)  ->  a b  +  a c

%  /.  
        Sin[a_  +  b_]  ->  Sin[a] Cos[b]  +  Cos[a] Sin[b]  //  
      Expand

ϕ[t_]  :=  Sin[2 π ω t];     ψ[t_]  :=  
      Cos[2 π ω t]

α  =  A Cos[2 π δ];  β  =  
      A Sin[2 π δ];
h[t]  ==  α ϕ[t]  +  β ψ[t]  //  
      Simplify

Thus h is a linear combination of the two functions ϕ and ψ.

Restatement of the problem of determining the hum

For lack of any better approach, we assume that the (unknown) hum function h is that linear combination of ϕ and ψ that is closest to the measured signal f.  What does "closest" mean here?  To answer that, we have to work in an appropriate inner product space, namely, C[0,  T] for some suitable time T.  (In other words, we shall be considering the input signal s as being generated for times 0  ≤  t  ≤  T—for a total time duration of T  -  0  =  T—and the output signal f being measured for those same times.)

The inner product space C[0,  T]

We already know how to define addition and multiplication of scalars for functions in such a way that the set C[0,  T] of all continuous functions with domain [0,  T] becomes a vector space.

Defining the innner product and norm

Defining inner product function

One-dimensional Mathematica input form:

ip[u_,  v_] := Integrate[
        u[t] v[t],  {t,  0,  T}]

Mathematica standard form for same thing:

ip[u_,  
        v_]  :=  (∫_0^T) u[t] v[
              t] \[DifferentialD]t

⚠  Dangerous way to define inner product using ordinary dot product notation:

Unprotect[Dot];
Dot[u_,  v_]  :=  ip[u,  v]

Dot[u,  v]

u  .  v

Defining norm function

norm[v_]  :=  Sqrt[v ~ip~ v]

The two basic period 2 π ω sinusoidal functions

ϕ ~ip~ ψ

We just used infix notation for ip in order to come as close as possible to the usual "dot" notation for the standard Euclidean inner product in (R^n).  The same thing may be entered in the usual prefix form as ip[ϕ, ψ].

norm[ϕ]^2

norm[ψ]^2

Assumption about duration T

From now on we assume that the duration T is an integral multiple of the period, that is, 

	 T  =  k ((1/ω)) with k some positive integer.

The formulas above for the inner product of ϕ with ψ and for the norms of these two functions now reduce to

	 ϕ · ψ  =  0,	 \(TraditionalForm\`\(\(\ \)\(‖ϕ‖\^2\  = \ \ \(‖ψ‖\^2\  = \ T\/2\)\)\)\).  

Thus the functions ϕ and ψ in C[0,  T] are orthogonal but do not necessarily have norm 1.

Using the assumption ω T  =  k with Mathematica Version 4

Mathematica Version 4 allows an optional argument to Simplify that specifies assumptions to be used in the simplification.  The assumption here is that ω T is an integer k.  (This does not work with earlier versions of Mathematica, which do not "know" as a built-in fact that the sine and cosine of an integral multiple of 2  π are 0 and 1, respectively!)

Simplify[ϕ ~ip~ ψ  /.  ω T  ->  k,  
      k  ∈  Integers]

Simplify[norm[ϕ]^2  /.  ω T  ->  k,  
      k  ∈  Integers]

Simplify[norm[ψ]^2  /.  ω T  ->  k,  
      k  ∈  Integers]

Normalizing to get orthonormal basis

We normalize ϕ and ψ to get orthonormal functions:

u[t_]  :=  (Sqrt[(2/(T ))]) Sin[120 π t]
v[t]  :=  (Sqrt[(2/(T ))]) Cos[120 π t]

Final restatement of the problem

Given f ϵ C[0,  T], find the hum as if it were the projection p of f onto Span(u,  v).

According to the Two-dimensional Projection Theorem, this projection p is given by the formula

	p  =  (f ·
                u)  u  +  (f·v)  v
     
(where the inner product, now being denoted with the usual dot, is the one defined above).

Examples

General setup for these examples

Clear[ω,  λ,  δ,  A,  T]

h[t_]  :=  
      A Sin[2 π (ω t  +  δ)]

f[t_]  :=  s[t]  +  h[t]

Below we include Re to take the real part of the inner products in order to try to prevent the spurious appearance of non-real complex numbers with very small imaginary parts.

p[t_]  :=  Re[f ~ip~ u] u[t]  +   Re[f ~ip~ v] v[t]

Example 1

This is our original example.

\(ω\  = \ 60; \ \ \ \ \ \ \ \ \  (*\ 
      frequency\ \ *) 
λ\  = \ 1/ω; \ \ \ \  (*\ 
      period\ \ \ \ \ *) 
δ\  = \  .05; \ \ \ \ \ \  (*\ 
      phase\ \ \ \ \ \ *) 

    A\  = \ 1.5; \ \ \ \ \ \  (*\ amplitude\ \ *) \)
T  =  1;

s[t_]  :=  
      0.01  Sin[120 π t]  +  0.02  Cos[120 π t]  +  
        4 (t  -  0.5)^2  - t

Our estimate of the hum is:

p[t]

Is this projection p really the original hum h?  And then does the "cure" after the "diagnosis" actually reconstruct the original signal s?  Let's compare graphs.  (We look just on the interval [0,  1/12] instead of the entire interval [0,  1] in order to cut down the number of cycles in the interval.)

Plot[Evaluate[p[t]],  {t,  0,  1/12}];

Notice that we wrapped the Mathematica function Evaluate around p[t] before plotting the graph.  This is so that the expression for p[t] will be calculated once and for all in terms of t, rather than calculated for each individual discrete value of t that Mathematica chooses at which to sample the function for plotting.  (Try it without the Evaluate to see what would happen.)

Plot[h[t],  {t,  0,  1/12}];

At first glance p and h look the same.  If the are the same, then f  -  p ought to be the same as the original ("unknown") input s.  Again, let's look to see:

Plot[Evaluate[f[t]  -  p[t]],  {t,  0,  1/12}];

Plot[s[t],  {t,  0,  1/12}];

Uh-oh!  The supposed cure didn't quite work: we have a reconstructed signal f  -  p that is only very roughly like the actual input signal s.  How far off are we?  Look at the distance between s and f  -  p in C[0,  T]:

norm[Function[t,  s[t]  -  (f[t]  -  p[t])]]

[We had to use a pure function there—which we could have denoted instead by

   (s[#] - (f[#] - p[#]))&

—owing to the way norm and ip are defined in terms of functions rather than in terms of expressions in t.)

Don't worry especially about the sudden appearance of complex numbers in the preceding result; it's a consequence of the particular way that Mathematica does integration.  If you wish, get rid of the (nearly) zero imaginary part by taking the real part of the result:

Re[%]

Is the calculated projection p really the unknown hum h?  Let's have a closer look.  First, plot them on the same axes:

\(Plot[Evaluate[{h[t], \ p[t]}], \ {t, \ 0, \ 1/12}, 

      PlotStyle\  -> \ {RGBColor[1, 0, 0], \ 
          RGBColor[0, 0, 1]}]; \)

If you look closely, you can see that the two graphs do not completely coincide.  To see this better, look at the difference:

Plot[Evaluate[h[t]  -  p[t]],  {t,  0,  1/12}];

Or, look at the distance between h and p in C[0,  T]:

norm[Function[t,  h[t]  -  p[t]]]

So—the projection p is NOT the same as the hum h—it differs from it by relatively small amounts (by what looks like another sinusoidal function, but one with an amplitude of only around 0.025 as contrasted with the amplitude 1.5 of the actual hum).

But the projection p is the best we can do—the closest sinusoidal function of frequency ω  =  60 that we can find to the actual hum h.  And hence the calculated f  -  p is the closest we can come to reconstructing the actual signal s.

Why could we not do better?  You'll have to think about this as you examine more examples.  Another example follows.

Example 2

This is example uses the same hum h as in the first example but a different input signal s.

\(\(\(ω\  = \ 60\)\(;\)\(\ \ \ \ \ \ \ \ \ \)\( (*\ 
      frequency\ \ *) \)\(
\)\(λ\  = \ 
      1/ω\)\(;\)\(\ \ \ \ \)\( (*\ 
      period\ \ \ \ \ *) \)\(
\)\(δ\  = \  .05\)\(;\)\(\ \ \ \ \ \ \)\( (*\ phase\ \ \ \ \ \ *) \)\(
\)\(A\  = \ 1.5\)\(;\)\(\ \ \ \ \ \ \)\( (*\ 
      amplitude\ \ *) \)\(
\)\(T\  = \ 
      1\)\(;\)\(\ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ 
      duration\ \ \ *) \)\)\)

s[t_]  :=  
      0.01  Sin[130 π t]  +  0.02  Cos[130 π t]  +  
        4 (t  -  0.5)^2  - t

Our estimate of the hum is:

p[t]

How good is our diagnosis—and hence the cure—in this example?

Plot[Evaluate[f[t]  -  p[t]],  {t,  0,  1/12}];

Plot[s[t],  {t,  0,  1/12}];

At first glance, things look pretty good!  Let's have a further look:

\(Plot[Evaluate[{s[t], \ f[t]\  - \ p[t]}], \ {t, \ 0, \ 1/12}, 

      PlotStyle\  -> \ {RGBColor[1, 0, 0], \ 
          RGBColor[0, 0, 1]}]; \)

Hardly perfect.  How badly or well did we do this time?

Plot[
      Evaluate[s[t]  -  (f[t]  -  p[t])],  {t,  0,  
        1/12}];

norm[Function[t,  s[t]  -  (f[t]  -  p[t])]]

Not perfect—but pretty good after all!  Certainly much better than in Example 1.