1 # Radix-2 DIT FFT R implementation
 2 # John Bryan, 2018
 3 rm(list=ls(all=TRUE))
 4 
 5 
 6 bitreversal<-function(x,N)
 7 {
 8     M1=2
 9     r=c()
10     r[1]=0
11     r[2]=1
12     # initialize r() in while loop
13     while  (M1 < N)
14     {
15          k = 0;
16          while (k < M1)
17          {
18                T1 = 2*r[k+1];
19                r[k+1] = T1;
20                r[k+M1+1] = T1+1;
21                k = k+1;
22          }
23          M1 = 2*M1;
24     }
25     cat("r  is " ,r)
26     cat("\n")
27     # end of initializing r()
28     # r(n1) is the bitreversal of n1.
29     # For n1=1,2,...,N-2 do: if r(n1) > n1, swap x(n1),x(r(n1)).
30     n1=1;
31     while (n1 < (N-1))
32     {
33         if ((r[n1+1])>n1)
34         {
35            temp=x[n1+1];
36            x[n1+1]=x[r[n1+1]+1];
37            x[r[n1+1]+1]=temp;
38         }
39         n1=n1+1;
40     }
41     return(x)
42 }
43 
44 
45 fft1<-function(x,T,N)
46 {
47     tss=N/2;        # initial twiddle step size
48     p=log2(N);
49     Bp=N/2;
50     Np=2;
51     for (P in 0:(p-1))        # pass loop
52     {
53       Npp=Np/2;               # butterflies per sub-block
54       BaseT=0;
55       for (b in 0:(Bp-1))     # sub=block loop
56       {
57          BaseB=BaseT+Npp;
58          for (n in 1:Npp)     # butterfly loop
59          {
60             top=x[BaseT+n];
61             bot=x[BaseB+n]*T[((n-1)*tss)+1];
62             x[BaseT+n]=top+bot;
63             x[BaseB+n]=top-bot;
64          }
65          BaseT=BaseT+Np;
66       }
67       Bp=Bp/2;        # subblocks halved
68       Np=Np*2;        # butterflies/sub-block doubles
69       tss=tss/2;      # twiddle step size halved
70    }
71    return(x)
72 }
73 
74 
75 options(digits=4)
76 N=64
77 t<-seq(from=0, to=(N-1)/N, by=1/N)
78 t0<-seq(from=0, to=((N/2)-1), by=1)
79 n<-seq(from=0, to=(N-1), by=1)
80 f=8
81 x=sin(2*pi*f*t)
82 x0 <- bitreversal(x,N)
83 T=c()
84 T[1] = 1+0i
85 for (m in 2:(N/2))
86     {T[m]=exp(2*pi*1i*(m-1)/N)}   # size N/2 twiddle table
87 X=fft1(x0,T,N)
88 X_norm=abs(X)/N
89 length(X_norm) <- N/2
90 jpeg('r1.jpeg')
91 par(mfrow=c(2,1),bg="grey")
92 plot(n,x,main="N=64 sampled 8 hz sine wave input x[n]",
93      xlab="n", ylab="x[n]", cex=0.8, cex.lab=0.9, cex.main=0.9, pch=1)
94 plot(t0,X_norm, main="|X[k]|",
95      xlab="k",ylab="|X[k]|", cex=0.8, cex.lab=0.9, cex.main=0.9,  pch = 2)
96 dev.off()