LES #1, #2 en #3 SOGM-theorie, Pieter Suurmond, 5-17 sept 2003.
    Blok 4, 3de jaar AD (Thomas en Sam + Tjeerd).
    ===============================================================

Je kunt een digitaal filter doorrekenen met als _enige_ kennis de 
rekenregels van machten:

 * Bij vermenigvuldiging van twee machten met gelijke grondtallen 
   mag je exponenten optellen.

Oftewel, bij complexe getallen:

 * Bij vermenigvuldiging van twee complexe getallen mag je hoeken 
   (oftewel argumenten) optellen en lengtes (oftewel moduli) 
   vermenigvuldigen.

Een concreet voorbeeld met een niet-recursief eerste orde 
laagdoorlaatfilter en 2 sinoide testsignalen:

 * Lineair, tijd-invariant systeem wordt beschreven door differentie-
   vergelijking  y[n] = (x[n] + x[n-1]) / 2.
 * Invoer testsignaal is het complexe wisselsignaal op halve Nyquist-
   frequentie (oftewel kwart samplerate)
                                                     (n PI j / 2)
                                             x[n] = e
   (waarbij j kwadraat = -1)
 * Gevraagd: In hoeverre wordt dit invoersignaal door dit systeem
   verzwakt of versterkt?
 * En daarna: Dezelfde vraag voor testsignaal op Nyquistfrequentie
 
                                                     (n PI j)
                                             x[n] = e

-----------------------------------------------------------------------------
De oplossing voor kwart samplerate:

  Substitutie x[n] in de differentievergelijking door het invoersignaal geeft

                 n PI j / 2    (n-1) PI j / 2
                e           + e
        y[n] = -------------------------------
                            2


                           j(n-1)    jn-j     -j  jn
  Volgens rekenregels is  e       = e      = e   e    zodat we 't kunnen her-
  schrijven als

                 n PI j / 2      -PI j / 2  n PI j / 2
                e             + e          e
        y[n] = ----------------------------------------
                              2

            -PI j / 2
  Zie dat  e          constant is, niet van n afhangt. Stel je het getal voor
  in het complexe vlak en zie, het is gelijk aan -j.
  
                 n PI j / 2        n PI j / 2              n PI j / 2 
                e             - j e                 (1-j) e 
        y[n] = -------------------------------  =  -------------------
                              2                               2


                          n PI j / 2
             = 0.5 (1-j) e 

  Het originele invoersignaal had magnitude 1, het uitvoersignaal heeft 
  magnitude (modulus)
                       | 0.5 (1-j) |  =  0.5 wortel(2)

  en is dus 3 deciBel in sterkte afgenomen.
  Als je deze laatste stap niet kunt maken, stel je dan het getal 0.5 (1-j)
  voor in het complexe vlak en gebruik de stelling van Pythagoras om de 
  lengte tot de oorsprong te berekenen.


-----------------------------------------------------------------------------
De oplossing voor halve samplerate:

  Substitutie x[n] in de differentievergelijking door het invoersignaal geeft

                 n PI j    (n-1) PI j        n PI j    n PI j
                e       + e                 e       - e
        y[n] = -----------------------  =  -------------------  = 0
                        2                           2

  Met andere woorden, een oneindige demping in deciBellen, Nyquist komt niet 
  meer uit het filter.

-----------------------------------------------------------------------------
Huiswerk voor les 2:

  * Bereken de demping voor beide frequenties (Nyquist en halve Nyquist) voor 
    het tweede orde laagdoorlaatfilter gekarakteriseerd door differentie-
    vergelijking
                    y[n] = (x[n] + 2 x[n-1] + x[n-2]) / 4.

  De antwoorden moeten zijn: een half en nul.
  Gebruik papier, potlood en gum, maak alle tussenstappen en _teken_ vooral 
  veel in het complexe vlak.

-----------------------------------------------------------------------------
Conclusie: zelfs zonder iets over 'polen', 'zeros' en 'Fourier' te weten kun 
je filters doorrekenen, volledig in het tijd-domein.
Wel is het misschien even wennen dat we als testsignaal geen pure reele sinus
nemen maar een 'complex wisselsignaal'. (Er zwabbert een imaginair deel 
mee, met precies dezelfde frequentie en amplitude, maar 90 graden verschoven 
in fase.)
-----------------------------------------------------------------------------

Oplossing huiswerk les 2:

Gegeven systeem  y[n] = (x[n] + 2 x[n-1] + x[n-2]) / 4
en complex wisselsignaal op kwart samplerate

                         (n PI j / 2)
                 x[n] = e             .

    Substitutie van x[n] in de differentievergelijking geeft

                          n PI j / 2      (n-1) PI j / 2    (n-2) PI j / 2
                         e           + 2 e               + e
                 y[n] = ---------------------------------------------------  .
                                                 4

    Uitvermenigvuldiging van de laatste twee exponenten geeft
    
                  n PI j/2      (n PI j/2 - PI j/2)    (n PI j/2 - PI j)
                 e         + 2 e                    + e
         y[n] = --------------------------------------------------------- .
                                         4

    Volgens rekenregel (vermenigvuldigen van machten met gelijke grondtallen 
    is gelijk aan optellen exponenten) kunnen we dat herschrijven als
    
                  n PI j/2      n PI j/2  -PI j/2    n PI j/2  -PI j
                 e         + 2 e         e        + e         e
         y[n] = --------------------------------------------------------- .
                                         4


                                n PI j/2
    De gemeenschappelijke term e         in de teller kan buiten haken worden 
    gebracht, wat 't geheel vereenvoudigt tot

                           n PI j/2         -PI j/2     -PI j
                   y[n] = e         (1 + 2 e        +  e     ) / 4 .

    De laatste twee e-machten zijn onafhankelijk van n en zien er in 
    cartesische notatie een stuk eenvoudiger uit. Het uitgangssignaal 
    
                           n PI j/2
                   y[n] = e         (1 - 2j -1) / 4


                           n PI j/2
                        = e         (-2j) / 4
                        
                                 n PI j/2
                        = -0.5j e

    heeft dezelfde frequentie als het ingangssignaal x[n].
    De amplitude echter, en ook de fase, zijn aangetast. De amplitude was 1 
    in het ingangssignaal. De amplitude van het uitgangssignaal is
    
                        | -0.5j | = 0.5

    wat overeenkomt met een versterking van 20 log  (1/2) = -6.0 dB.
                                                  10


-----------------------------------------------------------------------------
Huiswerk voor les #3:

  * Schrijf in een of andere programmeertaal (of pseudocode) een 
    LANGZAME Fourier transform routine.
    Eitje, 20 regels: geneste for-lus: eerste for-lus itereert harmonischen 
    tot aan Nyquist; tweede for itereert hapje audio. Reeel in, complex uit.
    Bereken magnitudes.
  
  * Vind (public) FFT broncode.

-----------------------------------------------------------------------------
Oplossing huiswerk voor les #3

  // SLOW FOURIER TRANSFORM in C-like pseudocode (without using comlex 
  // datatypes and complex.h). Exercise: rewrite to proper C++.

  const short N = 1024;     // Number of input samples, must be even.
  const short M = (N/2)+1;  // Number of output values (513).
  short       x[N];         // Array of (1024) 16-bit real input samples.
  long        y_re[M];      // Array containing (513) real output values.
  long        y_im[M];      // Array containing (513) imaginary output values.
  short       h, n;

  // Get input samples:
  for (n=0; n<N; n++)
      x[n] = read_from_somewhere(n);
      
  // Calculate Fourier transform:
  for (h=0; h<M; h++)             // Iterate harmonics (h), from DC to Nyquist.
      {
      y_re[h] = y_im[h] = 0;
      for (n=0; n<N; n++)         // Iterate discrete time (n).
          {
          y_re[h] += x[n] * cos(h * n * 2 * PI / N);
          y_im[h] += x[n] * sin(h * n * 2 * PI / N);
          }
      }

  // Deliver linear output magnitudes (no phases):
  for (h=0; h<M; h++)
      print_to_somewhere(sqrt((y_re[h] * y_re[h]) + (y_im[h] * y_im[h])));

-----------------------------------------------------------------------------
SOGM-theory #4:

INVERSE FOURIER (has no direct relation to LPC though, Sam):
The slow Fourier routine shown above, can NOT be used 
for inverse Fourier because that requires COMLEX input, not 
just REAL input. For inverse Fourier you need it the other 
way around: complex in --> real out.
When one uses complex-i/o for both one does not have to 
bother about this, it 's only that one performs twice as 
many multiplications and additions as strictly needed.

BONUS EXERCISE:
Analyze first order highpass filter that is characterised by 
linear difference equation
                            y[n] = (x[n] - x[n-1]) / 2

for 3 different test-frequencies with amplitude of 1:

                                (n PI j 0)
1) Zero Herz:         x[n] = 1 e

                                (n PI j / 2)
2) Quart samplerate:  x[n] = 1 e

                                (n PI j)
3) Half samplerate:   x[n] = 1 e

The outcomes should mirror the first order lowpass that we've 
previously studied.



-------------------------------------------------------------------------