optimointeja ja vanha fft takas

This commit is contained in:
2025-08-01 14:37:45 +03:00
parent 1d8d47b374
commit 266e4bd85d
4 changed files with 100 additions and 64 deletions

View File

@ -1,28 +1,78 @@
#include "fft.h" #include "fft.h"
#include <math.h> #include <math.h>
constexpr float PI = 3.14159; static float window[FFT_SIZE];
// In-place FFT on array of Complex numbers constexpr float M_PI = 3.14159;
void fft(Complex* x, int N, Complex* buffer) {
if (N <= 1) return;
Complex* even = buffer;
Complex* odd = buffer + N / 2;
for (int i = 0; i < N / 2; ++i) { // Call once before use
even[i] = x[i * 2]; void init_hamming_window() {
odd[i] = x[i * 2 + 1]; for (int i = 0; i < FFT_SIZE; i++) {
window[i] = 0.54f - 0.46f * cosf(2.0f * (float)M_PI * i / (FFT_SIZE - 1));
}
}
static unsigned int bit_reverse(unsigned int x, int log2n) {
unsigned int n = 0;
for (int i = 0; i < log2n; i++) {
n <<= 1;
n |= (x & 1);
x >>= 1;
}
return n;
}
void compute_fft(float* time_data, float* freq_out) {
static float real[FFT_SIZE];
static float imag[FFT_SIZE];
int log2n = 0;
for (int t = FFT_SIZE; t > 1; t >>= 1) ++log2n;
// Apply Hamming window
for (int i = 0; i < FFT_SIZE; i++) {
real[i] = time_data[i];// *window[i];
imag[i] = 0.0f;
} }
fft(even, N / 2, buffer + N); // deeper even // Bit reversal
fft(odd, N / 2, buffer + N + N / 2); // deeper odd for (int i = 0; i < FFT_SIZE; ++i) {
int j = bit_reverse(i, log2n);
if (j > i) {
float tmp_re = real[i], tmp_im = imag[i];
real[i] = real[j]; imag[i] = imag[j];
real[j] = tmp_re; imag[j] = tmp_im;
}
}
for (int k = 0; k < N / 2; ++k) { // Cooley-Tukey FFT
double angle = -2 * PI * k / N; for (int s = 1; s <= log2n; ++s) {
Complex twiddle(cos(angle), sin(angle)); int m = 1 << s;
Complex t = twiddle * odd[k]; for (int k = 0; k < FFT_SIZE; k += m) {
x[k] = even[k] + t; for (int j = 0; j < m / 2; ++j) {
x[k + N / 2] = even[k] - t; int t = k + j;
int u = t + m / 2;
float angle = -2.0f * (float)M_PI * j / m;
float w_real = cosf(angle);
float w_imag = sinf(angle);
float re = w_real * real[u] - w_imag * imag[u];
float im = w_real * imag[u] + w_imag * real[u];
real[u] = real[t] - re;
imag[u] = imag[t] - im;
real[t] += re;
imag[t] += im;
}
}
}
for (int i = 0; i < FFT_SIZE / 2; ++i) {
float mag = sqrtf(real[i] * real[i] + imag[i] * imag[i]) / FFT_SIZE;
float db = 20.0f * log10f(mag + 1e-6f); // Decibels
float normalized = (db + 60.0f) / 60.0f; // [0,1]
freq_out[i] = mag;
} }
} }

View File

@ -3,16 +3,6 @@
#include <math.h> #include <math.h>
#define FFT_SIZE 2048 #define FFT_SIZE 2048
// Simple complex number struct void init_hamming_window();
struct Complex { void compute_fft(float* time_data, float* freq_out);
float re, im;
Complex(float r = 0, float i = 0) : re(r), im(i) {}
Complex operator+(const Complex& o) const { return { re + o.re, im + o.im }; }
Complex operator-(const Complex& o) const { return { re - o.re, im - o.im }; }
Complex operator*(const Complex& o) const {
return { re * o.re - im * o.im, re * o.im + im * o.re };
}
};
void fft(Complex* x, int N, Complex* buffer);

View File

@ -40,8 +40,8 @@ static int pidPost;
#ifndef EDITOR_CONTROLS #ifndef EDITOR_CONTROLS
#pragma code_seg(".main") #pragma code_seg(".main")
// FFT buffers // FFT buffers
static Complex signal[FFT_SIZE]; static float fft_input[FFT_SIZE];
static Complex buffer[3 * FFT_SIZE]; static float fft_output[FFT_SIZE / 2];
static float fft_uniform[FFT_SIZE / 4]; static float fft_uniform[FFT_SIZE / 4];
void entrypoint(void) void entrypoint(void)
@ -123,6 +123,9 @@ int __cdecl main(int argc, char* argv[])
// Play sound // Play sound
direct_sound_buffer->Play(0, 0, 0); direct_sound_buffer->Play(0, 0, 0);
// Init FFT
init_hamming_window();
PFNGLUNIFORM1FVPROC glUniform1fvProc = ((PFNGLUNIFORM1FVPROC)wglGetProcAddress("glUniform1fv")); PFNGLUNIFORM1FVPROC glUniform1fvProc = ((PFNGLUNIFORM1FVPROC)wglGetProcAddress("glUniform1fv"));
const ULONGLONG targetIntervalMs = 1000 / 60; // For 60 FPS FFT updates const ULONGLONG targetIntervalMs = 1000 / 60; // For 60 FPS FFT updates
@ -178,7 +181,7 @@ int __cdecl main(int argc, char* argv[])
{ {
SUsample* samples = (SUsample*)audio_ptr; SUsample* samples = (SUsample*)audio_ptr;
for (int i = 0; i < FFT_SIZE; ++i) { for (int i = 0; i < FFT_SIZE; ++i) {
signal[i] = Complex((float)samples[i], 0.0); fft_input[i] = (float)samples[i];
} }
} }
@ -186,20 +189,18 @@ int __cdecl main(int argc, char* argv[])
} }
// Calculate FFT // Calculate FFT
fft(signal, FFT_SIZE, buffer); compute_fft(fft_input, fft_output);
// Normalize output // Normalize output
for (int i = 0; i < (FFT_SIZE / 4); i++) for (int i = 0; i < (FFT_SIZE / 4); i++)
{ {
float gain = 0.05f; float gain = 50.0f;
float alpha = 0.10f; // "Hidastaa" FFT:n piikkejä float alpha = 0.10f; // "Hidastaa" FFT:n piikkejä
float threshhold = 0.00015f; // Alin arvo mik<EFBFBD> p<EFBFBD><EFBFBD>stet<EFBFBD><EFBFBD>n shaderille (v<EFBFBD>hent<EFBFBD><EFBFBD> "noisea") float threshhold = 0.05f; // Alin arvo mikä päästetään shaderille (vähentää "noisea")
// float magnitude = sqrt(signal[i].re * signal[i].re + signal[i].im * signal[i].im); // signal strength float x_t = fft_output[i] * gain;
float magnitude = (float)signal[i].re;
float x_t = (magnitude < threshhold) ? 0.f : magnitude * gain;
// Exponential smoothing kaava // Exponential smoothing kaava
// s(t) = alpha*x(t)+(1-alpha)*s(t-1) // s(t) = alpha*x(t)+(1-alpha)*s(t-1)
fft_uniform[i] = alpha * (x_t)+(1 - alpha) * fft_uniform[i]; fft_uniform[i] = (x_t < threshhold) ? 0.f : alpha * (x_t)+(1 - alpha) * fft_uniform[i];
} }
} }

View File

@ -37,7 +37,8 @@ float noise(in vec2 xy, in float seed) {
// shorted functions // shorted functions
vec3 no(vec3 v) { return normalize(v); } vec3 no(vec3 v) { return normalize(v); }
float cl(float a, float b, float c) { return clamp(a,b,c); } float cl(float a, float b, float c) { return clamp(a,b,c); }
float le(vec3 s) { return length(s); } float lev3(vec3 s) { return length(s); }
float lev2(vec2 s) { return length(s); }
///////////////// /////////////////
@ -97,7 +98,7 @@ float hexPylon(vec3 p, vec2 h) {
p.xz = vec2(p.x * .866025 + p.z * .5, p.z); p.xz = vec2(p.x * .866025 + p.z * .5, p.z);
// The ".015" is a subtle rounding factor. Zero gives sharp edges, // The ".015" is a subtle rounding factor. Zero gives sharp edges,
// and larger numbers give a more rounded look. // and larger numbers give a more rounded look.
return le(max(abs(p) - b + .15, 0.)) - .15; return lev3(max(abs(p) - b + .15, 0.)) - .15;
} }
////////////// //////////////
@ -110,20 +111,20 @@ vec2 opU(vec2 d1, vec2 d2) {
} }
float sdSphere(vec3 p, float r){ float sdSphere(vec3 p, float r){
return length(p) -r; return lev3(p) -r;
} }
// Scene mapping with occlusion-aware SDF blending // Scene mapping with occlusion-aware SDF blending
vec2 mapScene(vec3 p) { vec2 mapScene(vec3 p) {
float dist = 20.; float dist = 20.;
int repeat = 0; int repeat = 0;
if((length(p.xz)) < 2*dist){ if((lev2(p.xz)) < 2*dist){
repeat = 2; repeat = 2;
} }
if((length(p.xz)) < 1.5*dist){ if((lev2(p.xz)) < 1.5*dist){
repeat = 3; repeat = 3;
} }
if((length(p.xz)) < dist){ if((lev2(p.xz)) < dist){
repeat = 5; repeat = 5;
} }
@ -135,15 +136,13 @@ vec2 mapScene(vec3 p) {
vec3 hexpos = vec3(p.x-dx, p.y-8.0, p.z-dy); vec3 hexpos = vec3(p.x-dx, p.y-8.0, p.z-dy);
HexData hex = hexTile(hexpos, 1.1); HexData hex = hexTile(hexpos, 1.1);
float distFromCenter = hexDistance(hex.axial); float distFromCenter = hexDistance(hex.axial);
int fftIndex = int(cl(distFromCenter + 1.0, 0.0, 511.0)); int fftIndex = int(cl(distFromCenter + 1., 0.0, 511.0));
float fftVal = fft_output[fftIndex];
float noise = mix(noise(hex.axial+1., 0.1), noise(hex.axial+1., 0.2), sin(syncs[0] * 4.)); float noise = mix(noise(hex.axial+1., 0.1), noise(hex.axial+1., 0.2), sin(syncs[0] * 4.));
float hexHeight = clamp(1.0 + fftVal * 5.0 + noise, 0., 15.); float hexHeight = cl(1.0 + (fft_output[fftIndex] * 2.0 + noise), 0., 15.);
vec3 r = vec3(hex.local.x + offset.x,hex.local.y,hex.local.z+offset.y); vec3 r = vec3(hex.local.x + offset.x,hex.local.y,hex.local.z+offset.y);
// r.yz *= rot2D(PI * 0.5); // r.yz *= rot2D(PI * 0.5);
r.xz *= rot2D(0.5); r.xz *= rot2D(0.5);
vec3 cellPos = r; float d = fHexagonCircumcircle(r, vec2(0.85, hexHeight));
float d = fHexagonCircumcircle(cellPos, vec2(0.85, hexHeight));
minDist = min(minDist, d); minDist = min(minDist, d);
} }
} }
@ -157,20 +156,16 @@ vec2 mapScene(vec3 p) {
vec3 castRay(vec3 ro, vec3 rd, inout vec3 pos) { vec3 castRay(vec3 ro, vec3 rd, inout vec3 pos) {
float mat = 0.; float mat = 0.;
float hit = 0.; float hit = 0.;
vec3 d; float t = 0.; // total distance travelled
float t = 0.,ad,tmax=200.; // total distance travelled
const float tolerance = 0.0001;
const float Z_REPEAT_DIST = 1.5; const float Z_REPEAT_DIST = 1.5;
vec2 res;
// Raymarching // Raymarching
for (int i = 0; i < 50; i++) { for (int i = 0; i < 50; i++) {
pos = ro + rd * t; pos = ro + rd * t;
res = mapScene(pos); // Get distance to objects vec2 res = mapScene(pos); // Get distance to objects
ad = abs(res.x);
mat = res.y; mat = res.y;
if (t > tmax) break; if (t > 200.) break;
if (ad < tolerance*(t*0.00125 + 1.0)) { if (abs(res.x) < 1e-4*(t*0.00125 + 1.0)) {
hit = 1.0; hit = 1.0;
break; break;
} }
@ -208,7 +203,7 @@ float softshadow(in vec3 ro, in vec3 rd, float mint, float maxt, float w) {
break; break;
float h = mapScene(ro + t * rd).x; float h = mapScene(ro + t * rd).x;
res = min(res, h / (w * t)); res = min(res, h / (w * t));
t += clamp(h, 0.1, 0.80); t += cl(h, 0.1, 0.80);
if(res < -1.0) if(res < -1.0)
break; break;
} }
@ -220,8 +215,8 @@ vec3 addPointLight(vec3 lightPos, vec3 lightColor, float intensity, vec3 worldPo
// Light vector from surface to light // Light vector from surface to light
float roughness = 1.0; float roughness = 1.0;
vec3 lightDir = lightPos - worldPos; vec3 lightDir = lightPos - worldPos;
float lightDistance = length(lightDir); float lightDistance = lev3(lightDir);
lightDir = normalize(lightDir); lightDir = no(lightDir);
// Attenuation (quadratic falloff) // Attenuation (quadratic falloff)
float attenuation = intensity / (1.0 + 0.09 * lightDistance + 0.032 * lightDistance * lightDistance); float attenuation = intensity / (1.0 + 0.09 * lightDistance + 0.032 * lightDistance * lightDistance);
@ -231,7 +226,7 @@ vec3 addPointLight(vec3 lightPos, vec3 lightColor, float intensity, vec3 worldPo
vec3 diffuse = lightColor * NdotL * attenuation; vec3 diffuse = lightColor * NdotL * attenuation;
// Specular lighting (Blinn-Phong) // Specular lighting (Blinn-Phong)
vec3 halfDir = normalize(lightDir + (-viewDir)); vec3 halfDir = no(lightDir + (-viewDir));
float NdotH = max(dot(normal, halfDir), 0.0); float NdotH = max(dot(normal, halfDir), 0.0);
float shininess = mix(128.0, 8.0, roughness); // Convert roughness to shininess float shininess = mix(128.0, 8.0, roughness); // Convert roughness to shininess
vec3 specular = lightColor * pow(NdotH, shininess) * attenuation; vec3 specular = lightColor * pow(NdotH, shininess) * attenuation;
@ -260,7 +255,7 @@ vec3 pal(float color) {
vec3 c = palette(color, 0.25, 0.63); vec3 c = palette(color, 0.25, 0.63);
vec3 c2 = palette(color, 0.5 ,0.2); vec3 c2 = palette(color, 0.5 ,0.2);
float t = clamp((syncs[0] - 129.0) / 2.0, 0.0, 1.0); // Smoothly ramps from 0 to 1 after 12s float t = cl((syncs[0] - 129.0) / 2.0, 0.0, 1.0); // Smoothly ramps from 0 to 1 after 12s
return mix(c, c2, t); // Blend between c and c2 over ~2 seconds return mix(c, c2, t); // Blend between c and c2 over ~2 seconds
//return c2; //return c2;
} }
@ -295,7 +290,7 @@ vec3 shading(vec3 p, vec3 n, vec3 dir, vec3 camPos) {
float particleStartPos = particlePos * 2.0 * PI; float particleStartPos = particlePos * 2.0 * PI;
// Direction from center to initial circle position (in XY plane) // Direction from center to initial circle position (in XY plane)
vec3 particleDir = normalize(vec3(cos(particleStartPos), 0.0, sin(particleStartPos))); // XZ direction vec3 particleDir = no(vec3(cos(particleStartPos), 0.0, sin(particleStartPos))); // XZ direction
float r = max(maxRadius - 0.0, maxRadius); float r = max(maxRadius - 0.0, maxRadius);