10 DEFDBL A-C,S,U-Z
20 DIM A(1025,2),U(5,2)
30 DIM W(5,2),T(5,2)
40 REM FFT DEMONSTRATION PROGRAM FOR N POINT TRANSFORM
50 PRINT "THIS PROGRAM EVALUATES THE MAGNITUDE AND PHASE"
60 PRINT "RESPONSE OF A NON-RECURSIVE Z TRANSFORM FILTER"
70 PRINT "USING THE FAST FOURIER TRANSFORM"
80 PRINT
90 PRINT "THE NUMBER OF POINTS IN THE FFT WILL BE 2^M"
100 PRINT
110 INPUT "TYPE M, OR  FOR M = 10; N = 2^M ";M
120 INPUT "TYPE SAMPLING FREQUENCY IN HERTZ ";SFREQ
130 INPUT "TYPE THE NUMBER OF TAPS FOR THE FIR FILTER ";NTAPS
140 REM THE COEFFICIENTS OF THE FIR FILTER ARE IN A(1,1), A(2,1),....
150 REM THE REST OF THE A(I,1) COEFFICIENTS ARE PADDED WITH ZEROES
160 REM AFTER THE FFT, A(I,1) CONTAINS THE REAL PART OF THE 
170 REM RESPONSE AT FREQUENCY W = (I-1)*3.14159/512
180 REM AND A(I,2) CONTAINS THE IMAGINARY PART OF THE RESPONSE
190 REM
200 FOR I = 1 TO NTAPS
210 PRINT "TYPE COEFFICIENT OF TAP ";I;TAB(2);
220 INPUT A(I,1)
230 NEXT I
240 REM
250 IF M = 0 THEN M=10
260 PRINT "MAG";TAB(16);"MAGDB";TAB(32);"HZ";TAB(48);"PHASE"
1000 N=2^M
1010 NV2 = N/2
1020 NM1 = N - 1
1030 J = 1
1040 FOR I = 1 TO NM1
1050 IF (I >= J) THEN GOTO 1090                  
1060 T(1,1) = A(J,1):T(1,2) = A(J,2)
1070 A(J,1)=A(I,1):A(J,2)=A(I,2)
1080 A(I,1)=T(1,1):A(I,2)=T(1,2)
1090 K = NV2
1100 IF ( K >= J) THEN GOTO  1140                   
1110 J = J - K
1120 K = K/2
1130 GOTO 1100                                  
1140 J = J + K
1150 NEXT I
1160 PI = 3.141592654#
1170 FOR L = 1 TO M
1180 LE=2^L
1190 LE1 = LE/2
1200 U(1,1)=1:U(1,2)=0
1210 W(1,1)=COS(PI/LE1):W(1,2)=SIN(PI/LE1)
1220 FOR J = 1 TO LE1
1230 FOR I = J TO N STEP LE
1240 IP = I + LE1
1250 SAVER = A(IP,1)*U(1,1)-A(IP,2)*U(1,2)
1260 SAVEI = A(IP,2)*U(1,1)+A(IP,1)*U(1,2)
1270 T(1,1)=SAVER:T(1,2)=SAVEI
1280 A(IP,1)=A(I,1)-T(1,1):A(IP,2)=A(I,2)-T(1,2)
1290 A(I,1)=A(I,1)+T(1,1):A(I,2)=A(I,2)+T(1,2)
1300 NEXT I
1310 SAVER = U(1,1)*W(1,1)-U(1,2)*W(1,2)
1320 SAVEI = U(1,2)*W(1,1)+U(1,1)*W(1,2)
1330 U(1,1)=SAVER:U(1,2)=SAVEI
1340 NEXT J
1350 NEXT L
5000 FOR I = 1 TO N/2+1
5010 RE = A(I,1):IM=A(I,2)
5020 MAG = SQR(ABS( (A(I,1)^2+A(I,2)^2)))
5030 IF MAG > 0 THEN MAGDB = 20*LOG(MAG)/LOG(10) ELSE MAGDB = -999
5040 IF RE > 0 THEN PHASE = ATN(IM/RE)*180/3.14159
5050 IF RE = 0 THEN PHASE = 90*SGN(IM)
5060 IF RE < 0 THEN PHASE = ATN(IM/RE)*180/3.14159+180
5070 IF PHASE > 180 THEN PHASE = PHASE -360:GOTO 5070
5080 IF PHASE < -180 THEN PHASE = PHASE + 360:GOTO 5080
5090 PRINT CSNG(MAG);TAB(16);CSNG(MAGDB);TAB(32);CSNG((I-1)*(SFREQ/N));TAB(48);CSNG(PHASE)
5100 IF I MOD 16 = 0 THEN GOSUB 6000
5110 NEXT I
5120 A$=INKEY$
5130 IF A$="" THEN GOTO 5120
5140 END
6000 A$=INKEY$
6010 IF A$="" THEN GOTO 6000
6020 RETURN