add functional lpf, remake hilbert

lib/filters.c

@@ -16,6 +16,82 @@ float apply_pre_emphasis(ResistorCapacitor *rc, float sample) {
     return audio;
 }
 
+void init_lpf(FrequencyFilter* filter, float cutoffFreq, float sampleRate) {
+    float nyquist = sampleRate / 2.0f;
+    float fc = cutoffFreq / nyquist;
+    
+    // Blackman window for sharp transition
+    for (int n = 0; n < FILTER_TAPS; n++) {
+        float m = n - (FILTER_TAPS - 1.0f) / 2.0f;
+        
+        // Sinc function
+        float sinc = (m == 0) ? 1.0f : sinf(PI * m * fc) / (PI * m);
+        
+        // Blackman window
+        float window = 0.42f - 0.5f * cosf(2.0f * PI * n / (FILTER_TAPS - 1)) 
+                     + 0.08f * cosf(4.0f * PI * n / (FILTER_TAPS - 1));
+        
+        filter->coeffs[n] = sinc * window;
+    }
+    
+    // Normalize
+    float sum = 0;
+    for (int i = 0; i < FILTER_TAPS; i++) sum += filter->coeffs[i];
+    for (int i = 0; i < FILTER_TAPS; i++) filter->coeffs[i] /= sum;
+    
+    // Clear delay line
+    memset(filter->delay, 0, sizeof(filter->delay));
+    filter->index = 0;
+}
+
+void init_hpf(FrequencyFilter* filter, float cutoffFreq, float sampleRate) {
+    float nyquist = sampleRate / 2.0f;
+    float fc = cutoffFreq / nyquist;
+    
+    // Blackman window for sharp transition
+    for (int n = 0; n < FILTER_TAPS; n++) {
+        float m = n - (FILTER_TAPS - 1.0f) / 2.0f;
+        
+        // Sinc function
+        float sinc = (m == 0) ? -1.0f : sinf(PI * m * fc) / (PI * m);
+        
+        // Blackman window
+        float window = 0.42f - 0.5f * cosf(2.0f * PI * n / (FILTER_TAPS - 1)) 
+                     + 0.08f * cosf(4.0f * PI * n / (FILTER_TAPS - 1));
+        
+        filter->coeffs[n] = sinc * window;
+    }
+
+    filter->coeffs[FILTER_TAPS/2] += 1.0f;
+    
+    // Normalize
+    float sum = 0;
+    for (int i = 0; i < FILTER_TAPS; i++) sum += filter->coeffs[i];
+    for (int i = 0; i < FILTER_TAPS; i++) filter->coeffs[i] /= sum;
+    
+    // Clear delay line
+    memset(filter->delay, 0, sizeof(filter->delay));
+    filter->index = 0;
+}
+
+float apply_freqeuncy_filter(FrequencyFilter* filter, float input) {
+    // Shift delay line
+    filter->delay[filter->index] = input;
+    
+    // Compute output
+    float output = 0;
+    int j = filter->index;
+    for (int i = 0; i < FILTER_TAPS; i++) {
+        output += filter->coeffs[i] * filter->delay[j];
+        j = (j + 1) % FILTER_TAPS;
+    }
+    
+    // Update index
+    filter->index = (filter->index + FILTER_TAPS - 1) % FILTER_TAPS;
+    
+    return output;
+}
+
 void init_delay_line(DelayLine *delay_line, int max_delay) {
     delay_line->buffer = (float *)calloc(max_delay, sizeof(float));
     delay_line->size = max_delay;

lib/filters.h

@@ -1,5 +1,5 @@
 #pragma once
-
+#define FILTER_TAPS 256
 #include <math.h>
 #include <string.h>
 #include <stdlib.h>
@@ -14,6 +14,15 @@ void init_rc(ResistorCapacitor *pe, float alpha);
 void init_rc_tau(ResistorCapacitor *pe, float tau, float sample_rate);
 float apply_pre_emphasis(ResistorCapacitor *pe, float sample);
 
+typedef struct {
+    float coeffs[FILTER_TAPS];
+    float delay[FILTER_TAPS];
+    int index;
+} FrequencyFilter;
+
+void init_lpf(FrequencyFilter* filter, float cutoffFreq, float sampleRate);
+float apply_freqeuncy_filter(FrequencyFilter* filter, float input);
+
 typedef struct {
     float *buffer;
     int write_idx;  // Write position

lib/hilbert.c

@@ -1,24 +1,43 @@
 #include "hilbert.h"
 
-void init_hilbert(HilbertTransformer *hilbert) {
-    hilbert->delay = calloc(D_SIZE, sizeof(float));
-    hilbert->dptr = 0;
+void compute_hilbert_coeffs(float* coeffs, int taps) {
+    int mid = taps / 2;
+    for (int i = 0; i < taps; i++) {
+        if ((i - mid) % 2 == 0) {
+            coeffs[i] = 0.0f;
+        } else {
+            coeffs[i] = 2.0f / (PI * (i - mid));
+        }
+    }
 }
 
-void apply_hilbert(HilbertTransformer *hilbert, float sample, float *output_0deg, float *output_90deg) {
-    float hilb;
+// Initialize the Hilbert transformer
+void init_hilbert(HilbertTransformer* filter) {
+    compute_hilbert_coeffs(filter->coeffs, HILBERT_TAPS);
+    memset(filter->delay, 0, sizeof(filter->delay));
+    filter->index = 0;
+}
 
-    hilbert->delay[hilbert->dptr] = sample;
-    hilb = 0.0f;
-    for(int i = 0; i < NZEROS/2; i++) {
-        hilb += (xcoeffs[i] * hilbert->delay[(hilbert->dptr - i*2) & (D_SIZE - 1)]);
+// Apply the Hilbert transformer
+void apply_hilbert(HilbertTransformer* filter, float input, float* inphase, float* quadrature) {
+    // Insert the new sample into the circular buffer
+    filter->delay[filter->index] = input;
+
+    // Compute the in-phase and quadrature components
+    float i_sum = 0.0f; // In-phase (0-degree output)
+    float q_sum = 0.0f; // Quadrature (90-degree output)
+
+    int coeff_index = 0;
+    for (int i = filter->index; coeff_index < HILBERT_TAPS; coeff_index++) {
+        i_sum += filter->delay[i] * (coeff_index == HILBERT_TAPS / 2 ? 1.0f : 0.0f);
+        q_sum += filter->delay[i] * filter->coeffs[coeff_index];
+
+        i = (i > 0) ? i - 1 : HILBERT_TAPS - 1;
     }
 
-    *output_0deg = hilbert->delay[(hilbert->dptr - 99) & (D_SIZE - 1)];
-    *output_90deg = hilb;
-    hilbert->dptr = (hilbert->dptr + 1) & (D_SIZE - 1);
-}
+    // Update the index for the next sample
+    filter->index = (filter->index + 1) % HILBERT_TAPS;
 
-void exit_hilbert(HilbertTransformer *hilbert) {
-    free(hilbert->delay);
+    *inphase = i_sum;
+    *quadrature = q_sum;
 }

lib/hilbert.h

@@ -1,44 +1,17 @@
 #pragma once
 
 #include <stdlib.h>
+#include <string.h>
+#include "constants.h"
+#include <math.h>
 
-#define D_SIZE 256
-#define NZEROS 200
-
-/* The non-zero taps of the Hilbert transformer */
-static float xcoeffs[] = {
-     +0.0008103736f, +0.0008457886f, +0.0009017196f, +0.0009793364f,
-     +0.0010798341f, +0.0012044365f, +0.0013544008f, +0.0015310235f,
-     +0.0017356466f, +0.0019696659f, +0.0022345404f, +0.0025318040f,
-     +0.0028630784f, +0.0032300896f, +0.0036346867f, +0.0040788644f,
-     +0.0045647903f, +0.0050948365f, +0.0056716186f, +0.0062980419f,
-     +0.0069773575f, +0.0077132300f, +0.0085098208f, +0.0093718901f,
-     +0.0103049226f, +0.0113152847f, +0.0124104218f, +0.0135991079f,
-     +0.0148917649f, +0.0163008758f, +0.0178415242f, +0.0195321089f,
-     +0.0213953037f, +0.0234593652f, +0.0257599469f, +0.0283426636f,
-     +0.0312667947f, +0.0346107648f, +0.0384804823f, +0.0430224431f,
-     +0.0484451086f, +0.0550553725f, +0.0633242001f, +0.0740128560f,
-     +0.0884368322f, +0.1090816773f, +0.1412745301f, +0.1988673273f,
-     +0.3326528346f, +0.9997730178f, -0.9997730178f, -0.3326528346f,
-     -0.1988673273f, -0.1412745301f, -0.1090816773f, -0.0884368322f,
-     -0.0740128560f, -0.0633242001f, -0.0550553725f, -0.0484451086f,
-     -0.0430224431f, -0.0384804823f, -0.0346107648f, -0.0312667947f,
-     -0.0283426636f, -0.0257599469f, -0.0234593652f, -0.0213953037f,
-     -0.0195321089f, -0.0178415242f, -0.0163008758f, -0.0148917649f,
-     -0.0135991079f, -0.0124104218f, -0.0113152847f, -0.0103049226f,
-     -0.0093718901f, -0.0085098208f, -0.0077132300f, -0.0069773575f,
-     -0.0062980419f, -0.0056716186f, -0.0050948365f, -0.0045647903f,
-     -0.0040788644f, -0.0036346867f, -0.0032300896f, -0.0028630784f,
-     -0.0025318040f, -0.0022345404f, -0.0019696659f, -0.0017356466f,
-     -0.0015310235f, -0.0013544008f, -0.0012044365f, -0.0010798341f,
-     -0.0009793364f, -0.0009017196f, -0.0008457886f, -0.0008103736f,
-};
+#define HILBERT_TAPS 95
 
 typedef struct {
-    float* delay;
-    int dptr;
+    float coeffs[HILBERT_TAPS];
+    float delay[HILBERT_TAPS];
+    int index;
 } HilbertTransformer;
 
-void init_hilbert(HilbertTransformer *hilbert);
-void apply_hilbert(HilbertTransformer *hilbert, float sample, float *output_0deg, float *output_90deg);
-void exit_hilbert(HilbertTransformer *hilbert);
+void init_hilbert(HilbertTransformer* filter);
+void apply_hilbert(HilbertTransformer* filter, float input, float* inphase, float* quadrature);