Turbo Toolbox/Turbo_Toolbox.iso/1988/02/fourier.mod

Modula Implementation | 1987-12-10 | 5.4 KB | 235 lines

IMPLEMENTATION MODULE Fourier ;
(*
   (C) Michael Kowalik & PASCAL INTERNATIONAL

   x   - vector of points projected on abscissa
   y   - vector of points projected on ordinate
   a,b - vectors of Fourier's coefficients
   a[i]*cos(2*pi*x[i]/T)
   b[i]*sin(2*pi*x[i]/T)
*)

FROM Conversions IMPORT ConvertAddrDec;
FROM MathLib0    IMPORT sin, cos ;
FROM SYSTEM      IMPORT ADDRESS, CODE, REGISTER, SETREG;

CONST
  MAXSIZE=512;       (* max number of vector elements *)
  S=MAXSIZE-1;
  pi=3.141592654;

PROCEDURE Interpolation( VAR x, y, a, b : ARRAY OF REAL);

VAR
  i, k, l, n, q : CARDINAL ;
BEGIN
  n := HIGH(y) + 1 ;
  q := n DIV 2 ;
  FOR i := 1 TO n DO
    x[i-1] := FLOAT(i-1)/FLOAT(n)*2.0*pi
  END ;
  a[0] := 0.0 ;
  FOR k := 0 TO n-1 DO
    a[0] := y[k] + a[0]
  END ;
  a[0] := 1.0/FLOAT(n)*a[0] ;
  FOR l := 1 TO q-1 DO
    a[l] := 0.0 ;
    b[l] := 0.0 ;
    FOR k := 0 TO n-1 DO
      a[l] := y[k]*cos(FLOAT(l)*x[k]) + a[l] ;
      b[l] := y[k]*sin(FLOAT(l)*x[k]) + b[l]
    END ;
    a[l] := 2.0/FLOAT(n)*a[l] ;
    b[l] := 2.0/FLOAT(n)*b[l]
  END ;
  a[q] := 0.0 ;
  FOR k := 0 TO n-1 DO
    a[q] := y[k]*cos(FLOAT(q)*x[k]) + a[q]
  END ;
  a[q] := 1.0/FLOAT(n)*a[q]
END Interpolation ;


PROCEDURE RungeFaltung(VAR x, y, a, b : ARRAY OF REAL);
CONST S=MAXSIZE;
      q=S DIV 2;
      p=q DIV 2;
VAR
     s                   : ARRAY [0..q] OF REAL;
     d                   : ARRAY [1..q-1] OF REAL;
     ss                  : ARRAY [0..p] OF REAL;
     sd                  : ARRAY [1..p] OF REAL;
     ds                  : ARRAY [0..p-1] OF REAL;
     dd                  : ARRAY [1..p-1] OF REAL;
     N, P, Q, n, i, k, l : CARDINAL;

BEGIN
  N:=HIGH(y)+1;
  Q:=N DIV 2;
  P:=Q DIV 2;
  FOR   i := 1 TO N DO
    x[i-1] := FLOAT(i-1)/FLOAT(N)*2.0*pi
  END;
  (* first convolution *)
  s[0]    := y[0];
  s[Q]    := y[Q];
  FOR i   := 1 TO Q-1 DO
    s[i]  := y[i] + y[N-i];
    d[i]  := y[i] - y[N-i]
  END;
  (* second convolution *)
  ss[P]   := s[P];
  sd[P]   := d[P];
  FOR i   := 0 TO P-1 DO
    ss[i] := s[i] + s[Q-i];
    ds[i] := s[i] - s[Q-i]
  END;
  FOR i   := 1 TO P-1 DO
    sd[i] := d[i] + d[Q-i];
    dd[i] := d[i] - d[Q-i]
  END;
  (* first and last cosinus coefficients *)
  a[0]    := 0.0;
  FOR   k := 0 TO P DO
    a[0]  := a[0] + ss[k]
  END;
  a[0]    := a[0]/FLOAT(N); a[Q]    := 0.0;
  FOR k   := 0 TO P DO
    IF ODD(k) THEN
      a[Q] := a[Q] - ss[k]
    ELSE
      a[Q] := a[Q] + ss[k]
    END
  END;
  a[Q]    := a[Q]/FLOAT(N);
  (* coefficients with odd  index *)
  FOR l   := 1 TO Q-1 BY 2 DO
    a[l]  := 0.0;
    FOR k := 0 TO P-1 DO
      a[l]:= a[l] + ds[k]*cos(FLOAT(l)*x[k])
    END;
    b[l]  := 0.0;
    FOR k := 1 TO P DO
      b[l]:= b[l] + sd[k]*sin(FLOAT(l)*x[k])
    END;
    a[l]  := 2.0*a[l]/FLOAT(N);
    b[l]  := 2.0*b[l]/FLOAT(N)
  END;
  (* coefficients with even index *)
  FOR l   := 2 TO Q-2 BY 2 DO
    a[l]  := 0.0;
    FOR k := 0 TO P DO
      a[l]:= a[l] + ss[k]*cos(FLOAT(l)*x[k])
    END;
    b[l]  := 0.0;
    FOR k := 1 TO P-1 DO
      b[l]:= b[l] + dd[k]*sin(FLOAT(l)*x[k])
    END;
    a[l]  := 2.0*a[l]/FLOAT(N);
    b[l]  := 2.0*b[l]/FLOAT(N)
  END
END RungeFaltung;


PROCEDURE DFT(VAR x, y, a, b : ARRAY OF REAL);
VAR
  alpha, c, s, suma, sumb, h, he : REAL;
  N, k, n                        : CARDINAL;
  Ar, Ai                         : ARRAY [0..S] OF REAL;
BEGIN
  N := HIGH(y);
  h := FLOAT(N+1);
  FOR k   := 0 TO N DO
    he    := FLOAT(k);
    suma  := 0.0; sumb  := 0.0;
    FOR n := 0 TO N DO
      alpha := 2.0*pi*he*FLOAT(n)/h;
      suma  := suma + y[n]*cos(alpha);
      sumb  := sumb + y[n]*sin(alpha);
    END;
    Ar[k] := 2.0*suma/h;
    Ai[k] := 2.0*sumb/h
  END;
  a[0] := Ar[0];
  b[0] := Ai[0];
  FOR k  := 1 TO N DIV 2 DO
    a[k] := Ar[k]+Ar[N+1-k];
    IF k<((N+1) DIV 2) THEN
      b[k] := Ai[k]-Ai[N+1-k]
    END
  END
END DFT;


PROCEDURE FFT(VAR Re, Im, a, b : ARRAY OF REAL);
VAR
  alpha, tr, ti, wr, wi, h        : REAL;
  L, N, i, j, m, mr, mn, l, istep : CARDINAL;
  Ar, Ai                          : ARRAY [0..S] OF REAL;


PROCEDURE BitRevShuffle(i, length : CARDINAL) : CARDINAL;
VAR
  m : CARDINAL;
BEGIN
  CODE(04287H);                    (* clr.l d7          *)
  SETREG(6, i);
  FOR m:=1 TO length DO
    CODE(0E24EH);                  (* lsr.w  #1,d6      *)
    CODE(0E357H);                  (* roxl.w #1,d7      *)
  END;
  RETURN CARDINAL(REGISTER(7))
END BitRevShuffle;


BEGIN
  N := HIGH(Re);
  h := FLOAT(N+1);
  L := 0;
  m := N;

  WHILE m>0 DO (* estimate number of bits *)
    m := m DIV 2;
    INC(L)
  END;

  FOR m := 0 TO N DO
    Ar[BitRevShuffle(m, L)] := Re[m];
    Ai[BitRevShuffle(m, L)] := Im[m]
  END;

  l := 1;
  WHILE l<N DO
    istep := 2*l;
    FOR m := 1 TO l DO
      alpha := pi*FLOAT(m-1)/FLOAT(l);
      wr  := cos(alpha);
      wi  := sin(alpha);
      i   := m - 1;
      REPEAT
        j := i+l;
        tr:= wr*Ar[j]-wi*Ai[j];
        ti:= wr*Ai[j]+wi*Ar[j];
        Ar[j] := Ar[i]-tr;
        Ai[j] := Ai[i]-ti;
        Ar[i] := Ar[i]+tr;
        Ai[i] := Ai[i]+ti;
        i:= i + istep
      UNTIL i>=N
    END;
    l := istep
  END;
  a[0] := Ar[0]/h;
  b[0] := Ai[0]/h;

  FOR m := 1 TO N DIV 2 DO
    a[m] := (Ar[m]+Ar[N+1-m])/h;
    IF m<((N+1) DIV 2) THEN
      b[m] := (Ai[m]-Ai[N+1-m])/h
    END
  END
END FFT

END Fourier.