|
|
OPTION BASE 0
OPTION ARITHMETIC COMPLEX
PUBLIC NUMERIC KK,XX(8500),YY(8500),X,Y,ALPHA
DIM F(8500),TX(8500),TY(8500)
CALL GINIT(800,800)
SET LINE COLOR 7
LET N=5
LET L=800/2^N
CALL JUMP(5,790)
CALL RECURSIVE(N,L,90)
LET LV=INT(LOG2(KK))+1
MAT TX=XX
MAT TY=YY
DO
MAT XX=TX
MAT YY=TY
FOR I=1 TO KK-1
LET F(I+1)=COMPLEX(XX(I+1)-XX(I),YY(I+1)-YY(I))
NEXT I
CALL FFT(2^LV,F)
LOCATE VALUE ,RANGE 2 TO 2^(LV-1)-1:NCUT
LET NCUT=INT(NCUT)
FOR I=NCUT+1 TO 2^LV-NCUT-1
LET F(I)=0
NEXT I
CALL IFFT(2^LV,F)
CLEAR
PLOT LINES
FOR I=1 TO KK
LET XX(I+1)=XX(I)+RE(F(I+1))
LET YY(I+1)=YY(I)+IM(F(I+1))
IF I=1 THEN PLOT LINES:XX(I),YY(I);
PLOT LINES:XX(I+1),YY(I+1);
NEXT I
LOOP
END
EXTERNAL SUB RECURSIVE(LEV,L,RR)
OPTION ARITHMETIC COMPLEX
IF LEV<>0 THEN
CALL TURN(RR)
CALL RECURSIVE(LEV-1,L,-RR)
CALL MOVE(L)
CALL TURN(-RR)
CALL RECURSIVE(LEV-1,L,RR)
CALL MOVE(L)
CALL RECURSIVE(LEV-1,L,RR)
CALL TURN(-RR)
CALL MOVE(L)
CALL RECURSIVE(LEV-1,L,-RR)
CALL TURN(RR)
END IF
END SUB
EXTERNAL SUB MOVE(L)
OPTION ARITHMETIC COMPLEX
LET X=X+L*COS(ALPHA*PI/180)
LET Y=Y-L*SIN(ALPHA*PI/180)
PLOT LINES:X,Y;
LET KK=KK+1
LET XX(KK)=X
LET YY(KK)=Y
END SUB
EXTERNAL SUB TURN(R)
OPTION ARITHMETIC COMPLEX
LET ALPHA=MOD(ALPHA+R+360,360)
END SUB
EXTERNAL SUB JUMP(XA,YA)
OPTION ARITHMETIC COMPLEX
LET X=XA
LET Y=YA
PLOT LINES:X,Y;
LET KK=KK+1
LET XX(KK)=X
LET YY(KK)=Y
END SUB
EXTERNAL SUB GINIT(XSIZE,YSIZE)
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE XSIZE,YSIZE
SET WINDOW 0,XSIZE-1,YSIZE-1,0
SET POINT STYLE 1
SET COLOR MODE "REGULAR"
SET COLOR MIX(0) 0,0,0
SET COLOR MIX(1) 0,0,1
SET COLOR MIX(2) 1,0,0
SET COLOR MIX(3) 1,0,1
SET COLOR MIX(4) 0,1,0
SET COLOR MIX(5) 0,1,1
SET COLOR MIX(6) 1,1,0
SET COLOR MIX(7) 1,1,1
CLEAR
END SUB
EXTERNAL SUB FFT0(N,Q,X(),Y())
OPTION ARITHMETIC COMPLEX
LET J=COMPLEX(0,1)
LET N0=0
LET N1=N/4
LET N2=N/2
LET N3=N1+N2
LET THETA0=2*PI/N
IF N=1 THEN
ELSEIF N=2 THEN
LET A=X(Q+0)
LET B=X(Q+1)
LET X(Q+0)=A+B
LET X(Q+1)=A-B
ELSEIF N>2 THEN
FOR P=0 TO N1-1
LET W1P=COMPLEX(COS(P*THETA0),-SIN(P*THETA0))
LET W2P=W1P*W1P
LET W3P=W1P*W2P
LET A=X(Q+P+N0)
LET B=X(Q+P+N1)
LET C=X(Q+P+N2)
LET D=X(Q+P+N3)
LET APC=A+C
LET AMC=A-C
LET BPD=B+D
LET JBMD=J*(B-D)
LET Y(Q+P+N0)=APC+BPD
LET Y(Q+P+N1)=W1P*(AMC-JBMD)
LET Y(Q+P+N2)=W2P*(APC-BPD)
LET Y(Q+P+N3)=W3P*(AMC+JBMD)
NEXT P
CALL FFT0(N/4,Q+N0,Y,X)
CALL FFT0(N/4,Q+N1,Y,X)
CALL FFT0(N/4,Q+N2,Y,X)
CALL FFT0(N/4,Q+N3,Y,X)
FOR P=0 TO N1-1
LET X(Q+4*P+0)=Y(Q+P+N0)
LET X(Q+4*P+1)=Y(Q+P+N1)
LET X(Q+4*P+2)=Y(Q+P+N2)
LET X(Q+4*P+3)=Y(Q+P+N3)
NEXT P
END IF
END SUB
EXTERNAL SUB IFFT0(N,Q,X(),Y())
OPTION ARITHMETIC COMPLEX
LET J=COMPLEX(0,1)
LET N0=0
LET N1=N/4
LET N2=N/2
LET N3=N1+N2
LET THETA0=2*PI/N
IF N=1 THEN
ELSEIF N=2 THEN
LET A=X(Q+0)
LET B=X(Q+1)
LET X(Q+0)=A+B
LET X(Q+1)=A-B
ELSEIF N>2 THEN
FOR P=0 TO N1-1
LET W1P=COMPLEX(COS(P*THETA0),SIN(P*THETA0))
LET W2P=W1P*W1P
LET W3P=W1P*W2P
LET A=X(Q+P+N0)
LET B=X(Q+P+N1)
LET C=X(Q+P+N2)
LET D=X(Q+P+N3)
LET APC=A+C
LET AMC=A-C
LET BPD=B+D
LET JBMD=J*(B-D)
LET Y(Q+P+N0)=APC+BPD
LET Y(Q+P+N1)=W1P*(AMC+JBMD)
LET Y(Q+P+N2)=W2P*(APC-BPD)
LET Y(Q+P+N3)=W3P*(AMC-JBMD)
NEXT P
CALL IFFT0(N/4,Q+N0,Y,X)
CALL IFFT0(N/4,Q+N1,Y,X)
CALL IFFT0(N/4,Q+N2,Y,X)
CALL IFFT0(N/4,Q+N3,Y,X)
FOR P=0 TO N1-1
LET X(Q+4*P+0)=Y(Q+P+N0)
LET X(Q+4*P+1)=Y(Q+P+N1)
LET X(Q+4*P+2)=Y(Q+P+N2)
LET X(Q+4*P+3)=Y(Q+P+N3)
NEXT P
END IF
END SUB
EXTERNAL SUB FFT(N,X())
OPTION ARITHMETIC COMPLEX
OPTION BASE 0
DIM Y(N)
CALL FFT0(N,0,X,Y)
MAT X=(1/N)*X
END SUB
EXTERNAL SUB IFFT(N,X())
OPTION ARITHMETIC COMPLEX
OPTION BASE 0
DIM Y(N)
CALL IFFT0(N,0,X,Y)
END SUB
|
|