#version 3.8; // Fast Fourier Transform (FFT), adapted to SDL from Paul Bourke's c++ code // Bill Walker "Bald Eagle" global_settings {assumed_gamma 1.0} #declare Red = rgb <1, 0, 0>; #declare Green = rgb <0, 1, 0>; #declare Blue = rgb <0, 0, 1>; #declare Yellow = rgb <1,1,0>; #declare Orange = rgb <1, 0.5, 0>; #declare Cyan = rgb <0, 1, 1>; #declare Magenta = rgb <1, 0, 1>; #declare Clear = rgbf 1; #declare White = rgb 1; #declare Black = rgb 0; // ################################################################################ #declare Zoom = 0.98; #declare AdjustCamera = 500; //max(Resolution.x, Resolution.y)*Y; camera { orthographic location right x*image_width/Zoom up y*image_height/Zoom look_at } // ################################################################################ plane {z, 0 pigment { gradient y translate -y*0.5 scale 500 color_map { [0.0 rgb <0.5, 0, 0>] [1 rgb <0, 0.5, 0>] } } } light_source { color rgb <1, 1, 1>} #declare SineData = array [64]; #declare Points = dimension_size (SineData, 1); #declare Cycles = tau/Points; #declare AA = 1*Cycles; #declare BB = 2*Cycles; #declare CC = 3*Cycles; #declare DD = 4*Cycles; #declare EE = 5*Cycles; #for (X, 0, Points-1 ) #declare SineData[X] = sin(AA*X) + sin(BB*X) + sin(CC*X) + sin(DD*X) + sin(EE*X); #end // end for X // determine if data array is a power of 2 // If yes, continue, if not, then expand array size to next higher power of 2 and fill with zeroes #declare PowerOf2 = int(log(Points) / log(2)); #debug concat ("Data input is ", str(Points, 8, 0), " points, which is 2 ^", str(log(Points)/log(2), 3, 0), "\n") #if (log(Points)/log(2) = PowerOf2) #debug concat(str(Points, 8, 0), " is an even power of 2. \n") #declare Count = Points; #declare Data = array [2][Count]; // copy data into array, with zeroes for imaginary components #for (X1, 0, Points-1) #declare Data [0][X1] = SineData [X1]; #declare Data [1][X1] = 0; #end // end for X1 #else #debug concat(str(Points, 8, 0), " is too small: expanding array to next largest power of 2. \n") #for (Test, 1, PowerOf2) #declare Count = pow(2, Test); #if (Count > Points) #debug concat("Creating an array with ", str(Count, 8, 0), " points... \n\n") #declare Data = array [2][Count]; // copy data into array, with zeroes for imaginary components #for (X1, 0, Points-1) #declare Data [0][X1] = SineData [X1]; #declare Data [1][X1] = 0; #end // end for X1 // fill out array with zeroes #for (X2, Points, PowerOf2-1) #declare Data [0][X2] = 0; #declare Data [1][X2] = 0; #end // end for X2 #end // end if #end // end for Test #end #declare cnt = Count; // size of [expanded] data array #declare sig = int(log(cnt) / log(2) + 0.9999); #declare real1 = pow (2, sig); #declare real = (real1 - 1); #declare real2 = int(real1 / 2); #declare real4 = int(real1 / 4); #declare real3 = real4 + real2; #debug concat ("real1 = ", str(real1, 2, 0), " real = ", str(real, 2, 0), " real2 = ", str(real2, 2, 0), " real4 = ", str(real4, 2, 0), " real3 = ", str(real3, 2, 0), "\n") #declare rel = array [real1]; #debug concat ("Array rel is ", str(dimension_size(rel, 1), 2, 0), " elements. 0 - ", str(dimension_size(rel, 1)-1, 2, 0), "\n") #declare img = array [real1]; #debug concat ("Array img is ", str(dimension_size(img, 1), 2, 0), " elements. 0 - ", str(dimension_size(img, 1)-1, 2, 0), "\n") #declare cmp = array [real3+1]; #debug concat ("Array cmp is ", str(dimension_size(cmp, 1), 2, 0), " elements. 0 - ", str(dimension_size(cmp, 1)-1, 2, 0), "\n\n") #for (i, 0, (cnt-1) ) #declare rel[i] = Data [0][i]; #declare img[i] = Data [1][i]; #end #declare sig2 = int(sig / 2); #declare sig1 = (sig - sig2); #declare cnt1 = pow(2, sig1); #declare cnt2 = pow(2, sig2); #debug concat ("sig2 = ", str(sig2, 2, 0), " sig1 = ", str(sig1, 2, 0), " cnt1 = ", str(cnt1, 2, 0), " cnt2 = ", str(cnt2, 2, 0), "\n") #declare V = array [cnt1]; #debug concat ("Array V is ", str(dimension_size(V, 1), 2, 0), " elements. 0 -", str(dimension_size(V, 1)-1, 2, 0), "\n\n") #declare V[0] = 0; #declare dv = 1; #declare ptr = cnt1; #for (j, 1, sig1) #debug concat ("j = ", str(j, 2, 0), "\n") #debug "==========\n" #declare hlfPtr = int(ptr/2); #declare pt = (cnt1-hlfPtr); #debug concat ("for i = ", str(hlfPtr, 2, 0), " to ", str(pt, 2, 0), " step ", str(ptr, 2, 0), " \n") #for (i, hlfPtr, pt, ptr) #debug concat (" i = ", str(i, 2, 0), " hlfPtr = ", str(hlfPtr, 2, 0), " i-hlfPtr = ", str(i-hlfPtr, 2, 0), " \n") #debug concat (" #declare V[", str(i, 2, 0), "] = V[", str(i-hlfPtr, 2, 0), "] + ", str(dv, 2, 0), " \n") #declare V[i] = V[i-hlfPtr] + dv; #end // end for i #declare dv = dv + dv; #declare ptr = hlfPtr; #debug "\n\n" #end // end for j #declare k = tau/real1; #debug concat ("for X = 0 to ", str(real4, 2, 0), " \n") #for (X, 0, real4) #debug concat ("X = ", str(X, 2, 0), "\n") #debug concat (" #declare cmp[", str(X, 2, 0), "] = cos(", str(k, 2, 0), "*", str(X, 2, 0), ") \n") #declare cmp[X] = cos(k*X); #debug concat (" #declare cmp[", str(real2, 2, 0), "-", str(X, 2, 0), "] = -", str(cmp[X], 2, 0), " \n") #declare cmp[real2 - X] = 0-cmp[X]; #debug concat (" #declare cmp[", str(real2, 2, 0), "+", str(X, 2, 0), "] = -", str(cmp[X], 2, 0), " \n") #declare cmp[real2 + X] = 0-cmp[X]; #end #debug "FFT: bit reversal \n \n" #for (i, 0, (cnt1-1) ) #declare ip = (i*cnt2); #for (j, 0, (cnt2-1) ) #declare h = (ip+j); #declare g = (V[j]*cnt2)+V[i]; #if (g > h) #declare temp = rel[g]; #declare rel[g] = rel[h]; #declare rel[h] = temp; #declare temp = img[g]; #declare img[g] = img[h]; #declare img[h] = temp; #end // end if #end // end for j #end // end for i #declare T = 1; #for (stage, 1, sig) #debug concat ("stage", str(stage, 2, 0), "\n\n") #declare d = int(real2 / T); #debug concat (" T = ", str(T, 2, 0), " d = ", str(d, 2, 0), "\n") #for (ii, 0, (T-1) ) #declare l = d*ii; #declare ls = l+real4; #debug concat (" l = ", str(l, 2, 0), " ls = ", str(ls, 2, 0), "\n") #for (i, 0, (d-1) ) #declare a = (2*i*T)+ii; #declare b = a+T; #debug concat (" a = ", str(a, 2, 0), " b = ", str(b, 2, 0), "\n") #declare f1 = rel[a]; #declare f2 = img[a]; #debug concat (" f1 = ", str(f1, 2, 0), " f2 = ", str(f2, 2, 0), "\n") #declare cnt1 = cmp[l] * rel[b]; #declare cnt2 = cmp[ls] * img[b]; #declare cnt3 = cmp[ls] * rel[b]; #declare cnt4 = cmp[l] * img[b]; #debug concat (" cnt1 = ", str(cnt1, 2, 0), " cnt2 = ", str(cnt2, 2, 0), " cnt3 = ", str(cnt3, 2, 0), " cnt4 = ", str(cnt4, 2, 0), "\n") #declare rel[a] = f1 + cnt1 - cnt2; #declare img[a] = f2 + cnt3 + cnt4; #declare rel[b] = f1 - cnt1 + cnt2; #declare img[b] = f2 - cnt3 - cnt4; //#debug concat ("sig2 = ", str(sig2, 2, 0), " sig1 = ", str(sig1, 2, 0), " cnt1 = ", str(cnt1, 2, 0), " cnt2 = ", str(cnt2, 2, 0), "\n") #end // end for i #end // end for ii #declare T = T + T; #end // end for stage #debug concat ("Coefficients:\n") #debug concat ("AA = ", str(AA*Points, 2, 3), "\n") #debug concat ("BB = ", str(BB*Points, 2, 3), "\n") #debug concat ("CC = ", str(CC*Points, 2, 3), "\n") #debug concat ("DD = ", str(DD*Points, 2, 3), "\n") #debug concat ("EE = ", str(EE*Points, 2, 3), "\n") #debug "\n\n Fourier Transform Results: \n" #debug " Num real imaginary \n" #for (i, 0, real) #if (abs(rel[i]) < pow(10,-5) ) #declare rel[i] = 0; #end // end if #if (abs(img[i]) < pow(10,-5) ) #declare img[i] = 0; #end // end if #debug concat (" ", str(i, 2, 0), " ", str(rel[i], 2, 8), " ", str(img[i], 2, 8), "\n" ) #end // end for i //======================================================================================================== #declare ScaleX = (image_width/Points); #declare ScaleY = ScaleX * 5; #declare FTScale = image_width*Cycles; #declare R = (image_width/Points)/6; #declare R2 = R/2; cylinder {<0, 0, 0>, , R2 pigment {rgb x*0}} cylinder {<0, 0, 0>, <0, image_height/2, 0>, R2 pigment {rgb y*0}} #for (X, 0, Points, 0.002) sphere { R2 texture {pigment {Red} } } sphere { R2 texture {pigment {Orange} } } sphere { R2 texture {pigment {Yellow} } } sphere { R2 texture {pigment {Green} } } sphere { R2 texture {pigment {Blue} } } sphere { R texture {pigment {Black} } } #end // end for X cylinder {, , R2 pigment {Yellow*0.5}} cylinder {, , R2 pigment {Yellow*0.5}} cylinder {, , R2 pigment {Yellow*0.5}} cylinder {, , R2 pigment {Yellow*0.5}} cylinder {, , R2 pigment {Yellow*0.5}} #for (i, 0, real) cylinder {<(image_width/i), 0, 0>, , R2 pigment {Magenta}} #end // end for i