home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / 1991 / 05 / praxis / spektrum.pas < prev   
Encoding:
Pascal/Delphi Source File  |  1991-04-04  |  9.5 KB  |  330 lines
  1. (* ------------------------------------------------------- *)
  2. (*                      SPEKTRUM.PAS                       *)
  3. (*       Programm zur Aufnahme vom Signalen und zur        *)
  4. (*            Berechnung ihres Frequenzspektrums           *)
  5. (*                                                         *)
  6. (*        (c) 1991 by Andreas Bartels  &  toolbox          *)
  7. (* ------------------------------------------------------- *)
  8. PROGRAM Spektrum;
  9.  
  10. {$IFDEF CPU87} {$N+}  {$ELSE} {$N-}  {$ENDIF}
  11.  
  12. USES
  13.       crt, dos, graph;
  14.  
  15. Const
  16.       MaxLaenge =  256;               (* = 512 SHR ShiftL *)
  17.       dXOben    =  250;
  18.       dX        =  128;         (* X-Achsen-Strichabstand *)
  19.       XStrOben  = '[ ms ]';
  20.       XStrUnten = '[ Hz ]';
  21.       YStrOben  =  '[ V ]';
  22.       YStrUnten =  '[/Hz]';
  23.  
  24. {$IFDEF CPU87}
  25.   TYPE REAL = EXTENDED;
  26. {$ENDIF}
  27.  
  28. TYPE
  29.       IntegerFeld = ARRAY[0..MaxLaenge] of INTEGER;
  30.       RealFeld    = ARRAY[0..MaxLaenge] of REAL;
  31.  
  32. VAR
  33.       GraphDriver,
  34.       GraphMode, MaxYWert,
  35.       XLinks, XRechts,
  36.       Y0Oben, Y0Unten,
  37.       Y1Oben, Y1Unten,
  38.       IMR21, IMRA1,
  39.       ShiftL, dXUnten, MYH,
  40.       MaxY, i, j, Hilf     : INTEGER;
  41.       Faktor, MaxWert      : REAL;
  42.       ch                   : CHAR;
  43.       HStr                 : STRING;
  44.       BitRev               : IntegerFeld;
  45.       Sinus, Cosinus,
  46.       Wert,
  47.       RealTeil, ImagTeil   : RealFeld;
  48.  
  49.  
  50. FUNCTION Fenster( i, FNr : INTEGER) : REAL;
  51.  
  52. BEGIN
  53.   CASE
  54.     FNr OF
  55.     0 :  { Rechteck }
  56.       Fenster := 1;
  57.     1 :  { Dreieck  }
  58.       Fenster := 1 - 2*abs(i-(maxlaenge SHR 1)) / maxlaenge;
  59.     2 :  { Hanning  }
  60.       Fenster := 0.5  - 0.5  * cos(2*Pi*i/maxlaenge);
  61.     3 :  { Hamming  }
  62.       Fenster := 0.54 - 0.46 * cos(2*Pi*i/maxlaenge);
  63.   ELSE END;
  64. END; { Fenster }
  65.  
  66.  
  67.  
  68.  
  69. PROCEDURE DatenLesen;
  70.  
  71. var   i, k : INTEGER;
  72.       c    : REAL;
  73.  
  74. BEGIN
  75.  
  76. (* Interrupts sperren *)
  77.   IMR21 := Port[$21];
  78.   IMRA1 := Port[$A1];
  79.   Port[$21]:= $FF;
  80.   Port[$A1]:= $FF;
  81.  
  82.                       (* Simulation aus Sinusschwingungen *)
  83.                                            (* 50 Hz-Sinus *)
  84. (*  FOR i := 0 TO Pred(MaxLaenge) DO
  85.     Wert[i] := 0.4*sin(dXOben*Pi*i/(2.5*Maxlaenge))
  86.                * Fenster(i,0); *)
  87.                                     (* Verrauschter Sinus *)
  88. (*  c := (63 + random) * 2 * Pi / MaxLaenge;
  89.   FOR i := 0 TO Pred(MaxLaenge) DO
  90.       Wert[i] := ((2*sin(c*i) + Random - Random) / 6)
  91.                  * Fenster(i,2); *)
  92.  
  93.                         (* 3 Sinusschwingungen überlagert *)
  94. (*  FOR i := 0 TO Pred(MaxLaenge) DO
  95.     Wert[i] := ( 60*Sin(i* 8*2*pi/MaxLaenge)
  96.                 +15*Sin(i*24*2*pi/MaxLaenge)
  97.                 + 5*Sin(i*40*2*pi/MaxLaenge))
  98.                 * Fenster(i,0)/120; *)
  99.  
  100.                                   (* A/D-Wandler abfragen *)
  101.   FOR i := 0 TO Pred(MaxLaenge) DO BEGIN
  102.     Wert[i] := PORT[$208] * Fenster(i,0);
  103.  
  104.         (* Verzögerung zur Einstellung der Samplefrequenz *)
  105.     FOR k := 0 TO 9 DO;
  106.   END;
  107.  
  108.                        (* Alte Interrupts wieder zulassen *)
  109.   Port[$21]:= IMR21;
  110.   Port[$A1]:= IMRA1;
  111.  
  112. END; { DatenLesen }
  113.  
  114.  
  115. PROCEDURE LookUpTable;
  116. (* Erstellt eine Tabelle für Sinus,Cosinus und Bit-Umkehr *)
  117.  
  118. VAR
  119.      Laenge, L12,
  120.      AdrNorm, AdrBRev  : INTEGER;
  121.      WinkelEinheit, Wi : REAL;
  122.  
  123. BEGIN
  124.  
  125. (* Sinus und Cosinus *)
  126.   WinkelEinheit := 2*Pi/MaxLaenge;
  127.   L12 := MaxLaenge SHR 1;
  128.   FOR i := 0 TO L12 DO BEGIN
  129.     Wi := WinkelEinheit*i;
  130.  
  131.         (* Symmetrieausnutzung bringt weiteren Zeitgewinn *)
  132.     Sinus[i]             := -sin(Wi);
  133.     Sinus[L12-i]         :=  Sinus[i];
  134.     Sinus[L12+i]         := -Sinus[i];
  135.     Sinus[MaxLaenge-i]   := -Sinus[i];
  136.     Cosinus[i]           :=  cos(Wi);
  137.     Cosinus[L12-i]       := -Cosinus[i];
  138.     Cosinus[L12+i]       := -Cosinus[i];
  139.     Cosinus[MaxLaenge-i] :=  Cosinus[i];
  140.   END; { FOR }
  141.  
  142.                                             (* Bit-Umkehr *)
  143.   AdrBRev   := 0;
  144.   BitRev[0] := 0;
  145.   FOR AdrNorm := 1 TO Pred(MaxLaenge) DO BEGIN
  146.     Laenge := MaxLaenge SHR 1;
  147.     WHILE AdrBRev+Laenge > Pred(MaxLaenge) DO
  148.        Laenge := Laenge SHR 1;
  149.     AdrBRev := AdrBRev Mod Laenge + Laenge;
  150.     IF AdrBRev > Pred(AdrNorm) Then BEGIN
  151.        BitRev[AdrNorm] := AdrBRev;
  152.        BitRev[AdrBRev] := AdrNorm;
  153.     END; { IF }
  154.   END; { FOR }
  155.  
  156. END; { LookUpTable }
  157.  
  158.  
  159.  
  160. PROCEDURE FFT( RealT : RealFeld );
  161.                         (* Sogenannte Hintransformation ! *)
  162.  
  163. VAR
  164.      TempReal, TempImag,
  165.      WichtungReal, WichtungImag : REAL;
  166.      TabNr, l, m, iSchritt      : INTEGER;
  167.  
  168. BEGIN
  169.  
  170.   FOR i := 0 TO Pred(MaxLaenge) DO BEGIN
  171.     RealTeil[i] := RealT[BitRev[i]];
  172.     ImagTeil[i] := 0;
  173.   END; { FOR }
  174.  
  175.                                      (* FFT - Algorithmus *)
  176.   l := 1;
  177.   WHILE l < MaxLaenge DO BEGIN
  178.     iSchritt := l SHL 1;
  179.     FOR m := 1 to l DO BEGIN
  180.       TabNr := ROUND(MaxLaenge SHR 1 Div l) * Pred(m);
  181.       WichtungReal := cosinus[TabNr];
  182.       WichtungImag := sinus[TabNr];
  183.       i := Pred(m);
  184.       REPEAT
  185.         j  := i+l;
  186.         TempReal  :=   WichtungReal * RealTeil[j]
  187.                      - WichtungImag * ImagTeil[j];
  188.         TempImag  :=   WichtungReal * ImagTeil[j]
  189.                      + WichtungImag * RealTeil[j];
  190.         RealTeil[j] := RealTeil[i] - TempReal;
  191.         ImagTeil[j] := ImagTeil[i] - TempImag;
  192.         RealTeil[i] := RealTeil[i] + TempReal;
  193.         ImagTeil[i] := ImagTeil[i] + TempImag;
  194.         i := i + iSchritt;
  195.       UNTIL i >= MaxLaenge
  196.     END; { FOR }
  197.     l := iSchritt;
  198.   END; { WHILE }
  199.  
  200. END; { FFT }
  201.  
  202.  
  203.  
  204. BEGIN { PROGRAM Spektrum }
  205.  
  206.   XLinks   :=  50;
  207.   XRechts  := 562;
  208.   Y0Oben   :=  17;
  209.   Y1Oben   := 160;
  210.   MaxYWert := Y1Oben - Y0Oben;
  211.   MYH      := MaxYWert DIV 2;
  212.   Y0Unten  := 182;
  213.   Y1Unten  := Y0Unten + MaxYWert;
  214.  
  215.   ShiftL   := 0;
  216.   Hilf     := MaxLaenge;
  217.   WHILE Hilf < 512 DO BEGIN
  218.     Hilf := Hilf SHL 1;
  219.     INC(ShiftL);
  220.   END;
  221.   dXUnten := (64 SHR ShiftL) * 250 DIV dXOben;
  222.  
  223.   FOR i := 0 TO MaxLaenge DO Wert[i] := 0;
  224.  
  225.   Writeln('Erstelle Look-up-table ... ');
  226.   LookUpTable;
  227.  
  228.   graphdriver := detect;
  229. (* falls nicht definiert, über "SET BGIPATH=PfadName" setzen *)
  230.   initgraph(graphdriver, graphmode, '');
  231.   RectAngle(  0, 0, 600, 347 );
  232.   OutTextXY( 55, 6, 'S P E K T R U M      Ver. 1.03'
  233.                   + '  (c) 1991 A. Bartels & toolbox');
  234.  
  235.                                      (* Beschriftung oben *)
  236.   RectAngle( XLinks,  Y0Oben, XRechts,  Y1Oben );
  237.   FOR i := 0 TO 4 DO BEGIN
  238.     Line( XLinks+i*dX, Y1Oben, XLinks+i*dX, Y1Oben+5 );
  239.     Str( (i*dXOben):3, HStr );
  240.     OutTextXY( XLinks+i*dX-15, Y1Oben+10, HStr ) ;
  241.   END;
  242.   OutTextXY( XLinks+ROUND(3.5*dX)-15, Y1Oben+10, XStrOben );
  243.   FOR i := -2 TO 2 DO BEGIN
  244.     Line( XLinks-5, Y1Oben-MYH-1-ROUND(i*MaxYWert/5),
  245.             XLinks, Y1Oben-MYH-1-ROUND(i*MaxYWert/5) );
  246.     Str( (i*0.5):4:1, HStr );
  247.     OutTextXY( XLinks-40, Y1Oben-MYH-ROUND((i*MaxYWert/5)+1),
  248.                HStr ) ;
  249.   END;
  250.   OutTextXY( XLinks-45, Y0Oben+ROUND(MaxYWert/6), YStrOben );
  251.  
  252.                                     (* Beschriftung unten *)
  253.  
  254.   RectAngle( XLinks, Y0Unten, XRechts, Y1Unten );
  255.   FOR i := 0 TO 4 DO BEGIN
  256.     Line( XLinks+i*dX, Y1Unten, XLinks+i*dX, Y1Unten+5 );
  257.     Str( (i*dXUnten):3, HStr );
  258.     OutTextXY( XLinks+i*dX-15, Y1Unten+10, HStr ) ;
  259.   END;
  260.  
  261.                                     (* 50 Hz - Markierung *)
  262.  
  263.   Line( ROUND(XLinks+(50 SHL Succ(ShiftL) * (dXOben/250))),
  264.         Y1Unten+2,
  265.         ROUND(XLinks+(50 SHL Succ(ShiftL) * (dXOben/250))),
  266.         Y1Unten+7 );
  267.   OutTextXY( XLinks+ROUND(3.5*dX)-15, Y1Unten+10,
  268.              XStrUnten );
  269.   FOR i := 0 TO 5 DO BEGIN
  270.     Line( XLinks-5, Y1Unten-ROUND(i*MaxYWert/5),
  271.             XLinks, Y1Unten-ROUND(i*MaxYWert/5) );
  272.     Str( i:2, HStr );
  273.     OutTextXY( XLinks-40, Y1Unten-ROUND((i*MaxYWert/5)+1),
  274.                HStr ) ;
  275.   END;
  276.   OutTextXY( XLinks-45, Y0Unten+ROUND(MaxYWert/10),
  277.              YStrUnten );
  278.  
  279. REPEAT
  280.  
  281.                       (* Aufnahme und Ausgabe des Signals *)
  282.     DatenLesen;
  283.     MaxWert := 1;
  284.     FOR i := 1 TO MaxLaenge-1 DO
  285.       IF Wert[i] > MaxWert THEN MaxWert := Wert[i];
  286.     SetViewPort( XLinks+1, Y0Oben+1,
  287.                  XRechts-1, Y1Oben-1, false );
  288.     ClearViewPort;
  289.     SetViewPort( XLinks, Y0Oben, XRechts, Y1Oben, false );
  290.     Faktor := MaxYWert / MaxWert;
  291.     Line( 0, MYH, 512, MYH );
  292.     MoveTo( 0, MaxYWert-MYH-ROUND(Wert[0]*Faktor) );
  293.     FOR i := 1 TO MaxLaenge DO
  294.       LineTO( i SHL ShiftL,
  295.               MaxYWert-MYH-ROUND(Wert[i]*Faktor) );
  296.  
  297.                          (* FFT und Ausgabe des Spektrums *)
  298.     FFT( Wert );
  299.     RealTeil[0] := 0;        (* Gleichanteil unterdrücken *)
  300.     ImagTeil[0] := 0;
  301.  
  302.    (* Aus dem komplexen- das Leistungs-Spektrum berechnen *)
  303.  
  304.     FOR i := 1 TO MaxLaenge SHR 1 DO
  305.       Wert[i] := SQRT( SQR(RealTeil[i]) + SQR(ImagTeil[i]) );
  306.     MaxWert := 1;
  307.     FOR i := 1 TO MaxLaenge SHR 1 DO
  308.       IF Wert[i] > MaxWert THEN MaxWert := Wert[i];
  309.     SetViewPort( XLinks+1, Y0Unten+1,
  310.                  XRechts-1, Y1Unten-1, false );
  311.     ClearViewPort;
  312.     SetViewPort( XLinks, Y0Unten, XRechts, Y1Unten, false );
  313.     Faktor := MaxYWert / MaxWert;
  314.     MoveTo( 0, MaxYWert - ROUND(Wert[0]*Faktor) );
  315.     FOR i := 1 TO MaxLaenge SHR 1 DO
  316.       LineTO( i SHL (ShiftL+1),
  317.               MaxYWert-ROUND(Wert[i]*Faktor) );
  318.  
  319.   UNTIL KeyPressed;
  320.  
  321.   ch := ReadKey;
  322.   REPEAT UNTIL KeyPressed;
  323.  
  324.   CloseGraph;
  325.  
  326. END.   { Spektrum }
  327.  
  328. (***********************************************************)
  329. (*                  Ende von SPEKTRUM.PAS                  *)
  330.