DataMuseum.dk

top - metrics - download

    Length: 3072 (0xc00)
    Types: TextFile
    Names: »fftipowatxt«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦this⟧ »fftipowatxt« 

TextFile

fftipowa=set 4
fftipowa=algol index.no
procedure for computation of power spectra of integer data
call: fftipowa(m,R);
    R real array (0:2**(m-1))
      at call 2**m = N integer input data x(j)
      R.j = x(j)   j integer field 0,2,4,...,2N-2
      at return N/2+1 real output data
      a(0)=R(0),   a(2q) = R(q),
      a(2q+1) = R(N/2+1-q)        q = 1,2,...,N/4
    m call value integer     2**m = number of input data.
external procedure fftipowa(m,R);
value m;
integer m;
real array R;
begin
     integer t,t1,t2,i,j,l,n,nn,q,p,pp;
     real r,s,sv,s2v,c,c2v;
     integer array field iaf;
     iaf:= 0; q:= 1 shift(m-1); n:= q-1; nn:= l:= p:= q+n;
     for l:= nn-l while l<>n do
     begin
        p:= nn-p;
        if l<n then
        begin l:= l+1; i:= -1;
           for i:= i+1 while p>=q do
           p:= p shift 1 extract m;
           p:= p add q shift (-i)
        end;
        if l<p then
        begin t:= R.iaf(l); R.iaf(l):= R.iaf(p);
           R.iaf(p):= t
        end;
     end l;
     sv:= 0.0; p:= 0; pp:= 1;
     for l:= 1 step 1 until m do
     begin s2v:= sv; s:= 0.0; c:= 1.0;
        sv:= case l of
           (1.0,
            0.70 710 678 1152,    
            0.38 268 343 2368,   
            0.19 509 032 2008,    
            0.98 017 140 3296' -1,
            0.49 067 674 3248' -1,
            0.24 541 228 5224' -1,
            0.12 271 538 2856' -1,
            0.61 358 846 4896' -2,
            0.30 679 567 6288' -2,
            0.15 339 801 8624' -2,
            0.76 699 031 8720' -3,
            0.38 349 518 7568' -3,
            0.19 174 759 7304' -3,
            0.95 873 799 0976' -4,
            0.47 936 899 6016' -4);
        c2v:= -2*sv*sv+c;
        q:= p-1; p:= pp; pp:= p+p;
        for i:= 0 step pp until nn do
        begin t:= R.iaf(i+p); R.iaf(i+p):= R.iaf(i)-t;
           R.iaf(i):= R.iaf(i)+t;
        end;
        for n:= 1 step 1 until q do
        begin r:= -s*s2v+c*c2v; s:= c*s2v+s*c2v; c:= r;
           for i:= n step pp until nn do
           begin j:= i-n-n+p; t:= R.iaf(i+p);
              t1:= R.iaf(j+p); t2:= t*c-t1*s; 
              t1:= t*s+t1*c; R.iaf(j+p):= R.iaf(j)+t1;
              R.iaf(i+p):= t1-R.iaf(j); R.iaf(j):= R.iaf(i)-t2;
              R.iaf(i):= R.iaf(i)+t2;
           end;
        end;
     end;
     nn:= nn+1; r:= nn; r:= r*nn; q:= p-1; t:= 0;
     for i:= 1 step 1 until q do
     begin t:= p-t+(i-1)extract 1;
        c:= R.iaf(i); s:= R.iaf(nn-i);
        R(t):= (c*c+s*s)/r;
     end;
     c:= R.iaf(0); R(0):= c*c/r;
     c:= R.iaf(p); R(t-1):= c*c/r;
end;
end
▶EOF◀