DSP 및 실시간 유틸리티 #

#include <std/dsp/dsp_all.h>  // master header: pulls in every module below

// complex_dsp.h
struct Complex32 {
    float re; float im;
    /* + - * / neg abs arg conj, as Complex */
    static Complex32 new(float re, float im);
    static Complex32 polar(float magnitude, float theta);
    static Complex32 from(Complex c);     // round each part to float
    Complex          to_complex() const;  // back to double (exact)
};

// fft.h
int fft_f32([]Complex32 data);
int ifft_f32([]Complex32 data);
int fft_serial_f32([]Complex32 data);
int ifft_serial_f32([]Complex32 data);
int fft_prepare_f32(unsigned long n);
int fft_real_f32([]const float in, []Complex32 out);
int fft_batch_f32([]Complex32 data, unsigned long n, int inverse);   // data.len / n transforms
int fft_real_batch_f32([]const float in, []Complex32 out, unsigned long n);
int fft_real_half_batch_f32([]const float in, []Complex32 out, unsigned long n);   // and fft_real_half_f32

// dsp.h
float dsp_dot_f32([]const float a, []const float b);   // f32x8 FMA accumulators

// filter.h, biquad.h, comb.h: each filter has process(x), reset() and
// process_block(in, out) on float samples, as its double counterpart
static struct FirFilter32 FirFilter32::new([]const float coeffs, []float history);
static struct IirFilter32 IirFilter32::new([]const float b, []const float a, []float xHistory, []float yHistory);
static struct Biquad32 Biquad32::new(float b0, float b1, float b2, float a1, float a2);
static struct Biquad32 Biquad32::from(struct Biquad bq);   // a double design, rounded once
static struct CombFF32 CombFF32::new([]float history, float gain);   // delay = history.len
static struct CombFB32 CombFB32::new([]float history, float gain);
static struct KarplusStrong32 KarplusStrong32::new([]float buffer, float damping);   // next(), pluck()

// convolution.h, window.h, stft.h, resample.h
int conv_direct_f32([]const float x, []const float h, []float out);
int conv_fft_f32([]const float x, []const float h, []float out);
void window_hann_f32([]float w);       // also _rectangular_, _hamming_, _blackman_f32
void window_apply_f32([]float x, []const float w);
int stft_forward_f32([]const float signal, unsigned long frameSize, unsigned long hopSize, int windowType, []Complex32 out);
int stft_inverse_f32([]const Complex32 frames, unsigned long frameSize, unsigned long hopSize, int windowType, []float outSignal);
void resample_sinc_f32([]const float in, []float out, unsigned long kernelHalfWidth);   // also resample_nearest_f32, resample_linear_f32

// dct.h, dft.h, imaging.h
int dct2_f32([]const float in, []float out);
int idct3_f32([]const float in, []float out);
int dft_f32([]const Complex32 in, []Complex32 out);   // also dft_real_f32, idft_f32
int conv2d_f32([]const float img, unsigned long h, unsigned long w, []const float kernel, unsigned long kh, unsigned long kw, []float out);
int fft2d_f32([]Complex32 img, unsigned long h, unsigned long w);    // and ifft2d_f32

#include <std/dsp/fixed.h>

Fixed  fixed_from_int(int x);       // x << 24
Fixed  fixed_from_float(double f);  // (Fixed)(f * 16777216.0)
int    fixed_to_int(Fixed x);       // x >> 24 (산술 시프트, -inf 방향으로 반올림)
double fixed_to_float(Fixed x);     // x / 16777216.0

Fixed fixed_add(Fixed a, Fixed b);
Fixed fixed_sub(Fixed a, Fixed b);
Fixed fixed_mul(Fixed a, Fixed b);   // (a * b) >> 24, 64비트 중간값
Fixed fixed_div(Fixed a, Fixed b);   // (a << 24) / b, 64비트 중간값
Fixed fixed_abs(Fixed x);
Fixed fixed_neg(Fixed x);
Fixed fixed_sqrt(Fixed x);           // 뉴턴-랩슨, 4회 반복; x는 0 이상이어야 함

#include <std/dsp/fixed.h>
#include <std/io.h>

int main() {
    Fixed a = fixed_from_float(3.5);
    Fixed b = fixed_from_float(2.0);
    Fixed product = fixed_mul(a, b);       // 7.0
    Fixed root    = fixed_sqrt(product);   // ~2.6458

    print("3.5 * 2.0 = ");
    println_float(fixed_to_float(product));
    print("sqrt of that ~= ");
    println_float(fixed_to_float(root));
    return 0;
}

#include <std/dsp/dsp.h>

Fixed dsp_dot([]const Fixed a, []const Fixed b);           // a[i]*b[i]의 합
void dsp_scale([]Fixed a, Fixed scale);            // a[i] *= scale
unsigned long dsp_add([]const Fixed a, []const Fixed b, []Fixed out); // out[i] = a[i]+b[i]
void dsp_clip([]Fixed a, Fixed lo, Fixed hi);      // [lo, hi]로 클램프
Fixed dsp_peak([]const Fixed a);                          // max |a[i]|
Fixed dsp_rms([]const Fixed a);                           // sqrt(mean(a[i]^2))

// SIMD FMA로 가속되는 두 *double* 배열의 내적 (simd.h 참고) --
// simd_fma_f64x4를 통해 반복당 4개 원소씩, n % 4에 대한 스칼라 나머지 처리.
// Fixed 대응 버전은 없다: Q8.24 곱셈에는 원소마다 64비트 확장 시프트가
// 필요한데, std::simd의 정수 FMA는 이를 제공하지 않는다. convolution.h의
// conv_direct가 커널을 뒤집어 연속적인 접근 패턴으로 만든 뒤 내부적으로
// 사용한다; 직접 사용하는 것도 가능하다.
double dsp_dot_f64([]const double a, []const double b);

// q < nOut에 대해 out[q] = dot(a, x[q .. q+n)): 여러 출력에 대한 FIR/상관.
// 한 번에 출력 32개(float는 64개), 탭마다 값을 네 개의 누산 벡터에
// 브로드캐스트: 수평 합이 없고, 곱셈-덧셈마다 로드가 하나입니다.
int dsp_sliding_dot_f64([]const double a, []const double x, []double out);
int dsp_sliding_dot_f32([]const float a, []const float x, []float out);

#include <std/dsp/complex_dsp.h>

struct Complex {
    double re; double im;
    Complex operator+(Complex other) const;
    Complex operator-(Complex other) const;
    Complex operator*(Complex other) const;
    Complex operator/(Complex other) const;
    Complex neg() const;    // 단항 부호 반전 (연산자 오버로드는 인자 개수로 구분되지 않음)
    double  abs() const;    // sqrt(re^2 + im^2)
    double  arg() const;    // atan2(im, re)
    Complex conj() const;   // (re, -im)

    static Complex new(double re, double im);
    static Complex polar(double magnitude, double theta);  // (mag*cos(theta), mag*sin(theta))
};

struct FComplex {   // Q8.24 대응 버전 -- Fixed re/im, operator/를 제외하면 형태 동일
    Fixed re; Fixed im;
    FComplex operator+(FComplex other) const;
    FComplex operator-(FComplex other) const;
    FComplex operator*(FComplex other) const;
    FComplex neg() const;
    Fixed     abs() const;
    FComplex  conj() const;

    static FComplex new(Fixed re, Fixed im);
    static FComplex from(Complex c);   // 설계 시점의 double 결과를 Q8.24로 양자화
    Complex         to_complex() const;
};

#include <std/dsp/dft.h>

int dft([]const Complex in, []Complex out);       // 순변환, 복소수 -> 복소수
int dft_real([]const double in, []Complex out);          // 순변환, 실수 -> 복소수
int idft([]const Complex in, []Complex out);      // 역변환 (1/N 포함)
int fdft([]const FComplex in, []FComplex out);    // Q8.24 순변환
int fidft([]const FComplex in, []FComplex out);   // Q8.24 역변환

#include <std/dsp/fft.h>

int dsp_is_pow2(unsigned long n);                                    // n >= 1이고 2의 거듭제곱이면
int fft([]Complex data);                      // 제자리 순변환, 1=성공/0=2의거듭제곱아님
int ifft([]Complex data);                     // 제자리 역변환 (1/N 포함)
int fft_serial([]Complex data);              // 호출한 스레드에서만 fft
int ifft_serial([]Complex data);
int fft_prepare(unsigned long n);                             // n 이하 길이의 트위들 테이블 생성
int fft_real([]const double in, []Complex out); // 실수 -> n점 전체 스펙트럼
// 연속으로 저장된 'batch'개 변환 (행 r은 data[r*n]), 스레드 사용
int fft_batch([]Complex data, unsigned long n, int inverse);   // data.len / n transforms
int fft_real_batch([]const double in, []Complex out, unsigned long n);   // batch = in.len / n
// 같은 작업을 제자리에서: 행 r의 앞쪽 n/2칸에 실수 n개가 들어 있음
int fft_real_batch_inplace([]Complex data, unsigned long n);
// 빈 0 .. n/2만(값 n/2 + 1개; 행 r은 out[r*(n/2+1)]부터)
int fft_real_half([]const double in, []Complex out);
int fft_real_half_batch([]const double in, []Complex out, unsigned long n);
int ffft([]FComplex data);                    // Q8.24 순변환
int fifft([]FComplex data);                   // Q8.24 역변환

#include <std/dsp/fft.h>
#include <std/io.h>

int main() {
    struct Complex data[8];
    for (int i = 0; i < 8; i++) data[i] = Complex::new((double)i, 0.0);

    fft(data[0..8]);     // in-place forward transform
    ifft(data[0..8]);    // in-place inverse — round-trips back to the original signal

    print("data[1].re after round-trip (expect ~1.0): ");
    println_float(data[1].re);
    return 0;
}

#include <std/dsp/convolution.h>

// O(len_x * len_h) -- 짧은 커널에 가장 적합. 출력 하나가 SIMD 내적 하나이며,
// 긴 입력은 스레드 풀에 나뉩니다.
int conv_direct([]const double x, []const double h, []double out);
int fconv_direct([]const Fixed x, []const Fixed h, []Fixed out);

// fft.h를 통한 O(n log n) -- 직접 컨볼루션의 O(len_x*len_h) 비용이 부담이
// 되는 긴 입력/커널에 가장 적합. double 정밀도 전용. 커널이 신호의 1/16
// 이하면 overlap-add(커널의 약 8배 블록을 fft_batch로), 그 밖에는 x와 h를
// 복소 FFT 하나에 함께 담아 변환합니다.
int conv_fft([]const double x, []const double h, []double out);

#include <std/dsp/window.h>

void window_rectangular([]double w);  // w[i] = 1 (테이퍼링 없음 -- 누설 기준선)
void window_hann([]double w);         // 0.5*(1 - cos(2*pi*i/(n-1)))
void window_hamming([]double w);       // 0.54 - 0.46*cos(2*pi*i/(n-1))
void window_blackman([]double w);      // 0.42 - 0.5*cos(2*pi*i/(n-1)) + 0.08*cos(4*pi*i/(n-1))
void window_apply([]double x, []const double w);  // x[i] *= w[i], 제자리 연산

#include <std/dsp/filter.h>

struct FirFilter {
    const &double coeffs;   // b[0..numTaps-1], caller-owned
    unsigned long numTaps;
    &double       history;  // caller-provided ring buffer, numTaps elements
    unsigned long pos;

    double process(double x);
    void   reset();          // zeros the history
    void   process_block([]const double in, []double out);   // out may alias in

    static struct FirFilter new([]const double coeffs, []double history);   // numTaps = min(coeffs.len, history.len)
};

struct IirFilter {           // general-order, direct form I
    const &double b; unsigned long numB;   // feedforward
    const &double a; unsigned long numA;   // feedback, a[0] assumed 1
    &double xHistory; &double yHistory;    // numB / numA elements, zeroed
    unsigned long xPos; unsigned long yPos;

    double process(double x);
    void   reset();
    void   process_block([]const double in, []double out);

    static struct IirFilter new([]const double b, []const double a,
                                []double xHistory, []double yHistory);
};

#include <std/dsp/filter.h>
#include <std/io.h>

int main() {
    double coeffs[3] = { 0.333333, 0.333333, 0.333333 };
    double history[3] = { 0.0, 0.0, 0.0 };
    struct FirFilter mavg = FirFilter::new(coeffs[0..3UL], history[0..3UL]);

    double step[6] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
    double out[6];
    mavg.process_block(step[0..6UL], out[0..6UL]);
    // out ~= [0.333, 0.667, 1, 1, 1, 1]
    println_float(out[2]);
    return 0;
}

#include <std/dsp/biquad.h>

struct Biquad {
    double b0, b1, b2, a1, a2, z1, z2;
    double process(double x);
    void   reset();
    // 블록 단위 process(), 상태는 레지스터에 유지; in == out 허용, 실시간 안전.
    // 분산 룩어헤드(scattered look-ahead, Parhi & Messerschmitt)를 쓴 직접형 I:
    // y[i] = q[i] + (a1^2 - 2 a2) y[i-2] - a2^2 y[i-4] (q는 입력만으로 계산)이므로
    // 짝수와 홀수 출력이 서로 독립인 두 사슬이 되고, 2레인 벡터 하나로 함께
    // 계산합니다. 추가되는 극점은 필터 자신의 극점을 뒤집은 것(반지름이 같음)이라
    // 안정성이 유지됩니다. process()와는 반올림 차이만 납니다.
    void   process_block([]const double in, []double out);

    // freqHz에서의 Z영역 주파수 응답: 크기(선형)와 위상(라디안).
    double response_mag(double fs, double freqHz) const;
    double response_phase(double fs, double freqHz) const;

    static struct Biquad new(double b0, double b1, double b2, double a1, double a2);

    // Robert Bristow-Johnson의 "Audio EQ Cookbook" 공식. fs = 샘플레이트,
    // f0 = 컷오프/중심 주파수, q = Q 계수 (높을수록 더 좁고 공진이 강함).
    static struct Biquad lowpass(double fs, double f0, double q);
    static struct Biquad highpass(double fs, double f0, double q);
    static struct Biquad bandpass(double fs, double f0, double q);   // 0dB 피크 게인 일정
    static struct Biquad notch(double fs, double f0, double q);
    static struct Biquad allpass(double fs, double f0, double q);
    static struct Biquad peaking(double fs, double f0, double q, double dbGain);
    static struct Biquad lowshelf(double fs, double f0, double s, double dbGain);   // s = 셸프 기울기, 기본값 1.0
    static struct Biquad highshelf(double fs, double f0, double s, double dbGain);

    // 쌍선형 변환을 통한 범용 아날로그-디지털 2차 섹션 변환:
    // H(s) = (B2 s^2 + B1 s + B0) / (A2 s^2 + A1 s + A0) -> 대응하는 디지털 biquad로.
    static struct Biquad bilinear(double B0, double B1, double B2,
                                  double A0, double A1, double A2, double fs);
};

struct FBiquad {   // Q8.24 대응 버전
    Fixed b0, b1, b2, a1, a2, z1, z2;
    Fixed process(Fixed x);
    void  reset();
    void  process_block([]const Fixed in, []Fixed out);

    static struct FBiquad from(struct Biquad bq);   // double로 설계된 biquad를 한 번 양자화
};

#include <std/dsp/biquad.h>
#include <std/io.h>

int main() {
    struct Biquad lp = Biquad::lowpass(48000.0, 1000.0, 0.707);

    // Unity DC gain — the standard invariant a lowpass must satisfy.
    print("|H(0)| (expect 1.0): ");
    println_float(lp.response_mag(48000.0, 0.0));

    double y = 0.0;
    for (int i = 0; i < 200; i++) y = lp.process(1.0);
    print("settled output for constant input 1.0 (expect ~1.0): ");
    println_float(y);

    struct Biquad32 lp32 = Biquad32::from(lp);   // 같은 설계를 float로 한 번 반올림
    struct FBiquad  lpq  = FBiquad::from(lp);    // ... 또는 Q8.24로 한 번 양자화
    return 0;
}

#include <std/dsp/dct.h>

int dct2([]const double in, []double out);    // X[k] = sum x[n]*cos(pi/N*(n+0.5)*k)
int idct3([]const double in, []double out);   // dct2의 정확한 역변환
int fdct2([]const Fixed in, []Fixed out);      // Q8.24 순변환
int fidct3([]const Fixed in, []Fixed out);     // Q8.24 역변환

#include <std/dsp/stft.h>

#define STFT_WIN_RECTANGULAR 0
#define STFT_WIN_HANN 1
#define STFT_WIN_HAMMING 2
#define STFT_WIN_BLACKMAN 3

unsigned long stft_num_frames(unsigned long signalLen, unsigned long frameSize, unsigned long hopSize);

// out: stft_num_frames(...) * frameSize개의 Complex 값을 담을 공간.
int stft_forward([]const double signal, unsigned long frameSize, unsigned long hopSize, int windowType, []Complex out);

// outSignal: (numFrames-1)*hopSize + frameSize개 샘플을 담을 공간, 호출자가
// 0으로 초기화 (stft_inverse가 그 안에 누적함).
int stft_inverse([]const Complex frames, unsigned long frameSize, unsigned long hopSize, int windowType, []double outSignal);

double mr_stft_loss([]const double x, []const double y, []const unsigned long frameSizes);

#include <std/dsp/resample.h>

void resample_nearest([]const double in, []double out);
void resample_linear([]const double in, []double out);

// 윈도우가 적용된 sinc(대역 제한): 가장 품질이 높은 옵션. 다운샘플링할
// 때는 커널의 컷오프가 새로운(더 낮은) 나이퀴스트 레이트에 맞춰 자동으로
// 스케일되므로, 앤티에일리어싱 필터 역할도 겸합니다. kernelHalfWidth의
// 기본값으로는 8-16이 적당합니다. 탭마다의 가중치는 각의 차 공식으로
// 계산하고, 탭 각도 테이블은 호출마다 한 번 만듭니다. 그래서 탭끼리
// 독립이 되어 벡터 하나에 네 개씩 실행되며, 출력마다의 각도 네 개는 출력에서
// 출력으로 회전으로 넘기고 64개마다 정확히 다시 계산합니다. 출력은 스레드
// 풀에서 실행됩니다.
void resample_sinc([]const double in, []double out, unsigned long kernelHalfWidth);

// Q8.24 nearest/linear. sinc 대응 버전은 없다 -- 윈도우 적용된 sinc
// 커널은 런타임의 소수 위치에 대해 탭마다 sin/cos를 계산해야 하는데,
// 이는 일반적인 설계 시점-float/런타임-fixed 분리 방식에 맞지 않는다.
void fresample_nearest([]const Fixed in, []Fixed out);
void fresample_linear([]const Fixed in, []Fixed out);

#include <std/dsp/ztransform.h>

// B, A: order+1개의 아날로그(s영역) 계수, s의 오름차순 거듭제곱.
// bOut, aOut: 음의 z 거듭제곱 표기법(aOut[0] = 1)에서의 order+1개
// 디지털 계수 -- filter.h의 IirFilter::new에 그대로 넣으면 된다.
int bilinear_nth_order([]const double B, []const double A, double fs, []double bOut, []double aOut);   // order = B.len - 1

double ztransform_response_mag([]const double b, []const double a, double fs, double freqHz);
double ztransform_response_phase([]const double b, []const double a, double fs, double freqHz);

#include <std/dsp/comb.h>

struct CombFF {
    &double history; unsigned long delay; unsigned long pos; double gain;
    double process(double x);
    void   reset();
    void   process_block([]const double in, []double out);   // out may alias in

    // 'history' is the caller's ring buffer (its length is the delay), zeroed.
    static struct CombFF new([]double history, double gain);   // delay = history.len
};
struct CombFB {
    &double history; unsigned long delay; unsigned long pos; double gain;
    double process(double x);
    void   reset();
    void   process_block([]const double in, []double out);

    static struct CombFB new([]double history, double gain);   // delay = history.len
};

struct KarplusStrong {
    &double buffer;         // ring buffer; also the delay line and the pitch period
    unsigned long length; unsigned long pos;
    double damping;         // 0..1, closer to 1 = slower decay (typical 0.99-0.999)

    double next();          // the next output sample: no input, self-sustaining
    // Refill the delay line with a new excitation (n samples copied in; size
    // it to 'length' for a fresh pluck, the rest keeps its old contents).
    void   pluck([]const double excitation);

    // 'buffer' is the caller's delay line (its length is the period); seed it (or
    // call pluck) before next().
    static struct KarplusStrong new([]double buffer, double damping);   // length = buffer.len
};

#include <std/dsp/comb.h>
#include <std/io.h>

int main() {
    double excitation[8] = { 0.5, -0.3, 0.8, -0.6, 0.2, -0.9, 0.4, -0.1 };  // noise burst
    double buf[8];
    struct KarplusStrong ks = KarplusStrong::new(buf[0..8UL], 0.995);
    ks.pluck(excitation[0..8UL]);

    for (int i = 0; i < 16; i++) {
        double sample = ks.next();
        println_float(sample);   // decaying, self-sustained "pluck"
    }

    double block[32];
    generate_block(&ks, block[0..32UL]);   // 다음 32샘플: 모든 SampleGenerator에 동작
    return 0;
}

#include <std/dsp/sample_filter.h>   // filter.h, biquad.h, comb.h도 포함합니다

#include <std/dsp/filter.h>
#include <std/io.h>

// 1극 DC 차단기: y[n] = x[n] - x[n-1] + r*y[n-1].
struct DcBlocker {
    double r; double x1; double y1;
    double process(double x);
    void   process_block([]const double in, []double out);
    void   reset();
};
double DcBlocker::process(double x) {
    double y = x - self.x1 + self.r * self.y1;
    self.x1 = x;
    self.y1 = y;
    return y;
}
void DcBlocker::process_block([]const double in, []double out) {
    filter_process_each(self, in, out);
}
void DcBlocker::reset() { self.x1 = 0.0; self.y1 = 0.0; }

// 모든 SampleFilter: 단위 계단 입력 n샘플 뒤의 출력.
generic<F: SampleFilter>
double settled(&F f) {
    double s[256];
    filter_step_response(f, s[0..256UL]);
    return s[255];
}

int main() {
    struct DcBlocker dc;
    dc.r = 0.995;
    dc.reset();

    double taps[3] = { 0.25, 0.5, 0.25 };
    double hist[3] = { 0.0, 0.0, 0.0 };
    struct FirFilter fir = FirFilter::new(taps[0..3UL], hist[0..3UL]);

    println_float(settled(&dc));    // ~0.28이며 계속 감소: DC 차단기는 DC를 제거합니다
    println_float(settled(&fir));   // 1.0: 탭의 합
    return 0;
}

generic<S, A: SampleFilter, B: SampleFilter>
struct FilterChain {
    A first;
    B second;

    S    process(S x);     // second.process(first.process(x))
    void reset();          // 두 단 모두
    // 각 단이 자신의 process_block으로 블록 전체를 처리하고(각자의 빠른 경로 유지),
    // 'out'을 중간 버퍼로 씁니다: 별도 스크래치 없음, out은 in과 겹쳐도 됨.
    void process_block([]const S in, []S out);

    static struct FilterChain<S, A, B> of(A first, B second);
};

#include <std/dsp/biquad.h>
#include <std/dsp/filter.h>
#include <std/io.h>

int main() {
    // 80 Hz - 8 kHz 대역: biquad 두 개로 만든 4차 응답.
    struct Biquad hp = Biquad::highpass(48000.0, 80.0, 0.7071);
    struct Biquad lp = Biquad::lowpass(48000.0, 8000.0, 0.7071);
    struct FilterChain<double, struct Biquad, struct Biquad> band = FilterChain::of(hp, lp);

    // 체인은 SampleFilter이므로 다른 체인의 한 단이 될 수 있습니다.
    double taps[3] = { 0.25, 0.5, 0.25 };
    double hist[3] = { 0.0, 0.0, 0.0 };
    struct FirFilter smooth = FirFilter::new(taps[0..3UL], hist[0..3UL]);
    struct FilterChain<double, struct FilterChain<double, struct Biquad, struct Biquad>, struct FirFilter> eq =
        FilterChain::of(band, smooth);

    double buf[256];
    for (int i = 0; i < 256; i++) buf[i] = (i % 32 < 16) ? 1.0 : -1.0;   // 사각파
    eq.process_block(buf[0..256UL], buf[0..256UL]);                              // 제자리 처리
    println_float(buf[255]);
    return 0;
}

#include <std/dsp/minphase.h>

// x: n개 샘플 (n <= fftLen); fftLen은 2의 거듭제곱이어야 하며, 정확한
// 켑스트럴 리프터링을 위해 n보다 충분히 커야 함 (경험적으로 4-8배 정도가
// 적당함). outMinPhase: fftLen개 샘플을 담을 공간 -- 에너지가 시작 부분에
// 집중되며, 호출자는 보통 처음 n개 정도의 샘플만 남겨서 사용함.
int minimum_phase([]const double x, unsigned long fftLen, []double outMinPhase);

#include <std/dsp/cqt.h>

// signal은 최소한 가장 낮은(가장 긴 커널을 갖는) 빈에 필요한 만큼의
// 샘플을 가져야 한다 -- 대략 Q*fs/fmin, Q = 1/(2^(1/binsPerOctave)-1) --
// 그렇지 않으면 해당 빈의 커널이 signalLen에 맞춰 조용히 잘린다
// (오류가 아니라 해당 빈의 해상도 저하). out: numBins개의 Complex 값을
// 담을 공간.
void cqt_forward([]const double signal, double fs, double fmin, unsigned long binsPerOctave, []Complex out);

#include <std/dsp/cwt.h>

// scales: numScales개 항목 (각각 0보다 커야 함). w0: Morlet 중심 주파수
// 파라미터 (5-6이 표준). out: numScales*signalLen개의 Complex 값을 담을
// 공간, 스케일 기준 행 우선(row-major) 순서.
int cwt_morlet([]const double signal, []const double scales, double w0, []Complex out);

#include <std/dsp/imaging.h>

// img: h*w개, 행 우선; kernel: kh*kw개, 행 우선 (kh/kw는 홀수여야 함);
// out: h*w개의 double을 담을 공간. 범위를 벗어난 커널 탭은 0으로 취급.
int conv2d([]const double img, unsigned long h, unsigned long w, []const double kernel, unsigned long kh, unsigned long kw, []double out);

int fft2d([]Complex img, unsigned long h, unsigned long w);   // 제자리 연산, h와 w는 2의 거듭제곱이어야 함
int ifft2d([]Complex img, unsigned long h, unsigned long w);

int gaussian_kernel([]double kernel, unsigned long size, double sigma);  // 정규화됨 (합이 1)
int sobel_x_kernel([]double kernel);   // 고정 3x3, 9개의 double을 담을 공간
int sobel_y_kernel([]double kernel);

#include <std/dsp/audio_buffer.h>

struct AudioBuffer {
    Fixed         data[AUDIO_BUF_FRAMES * AUDIO_MAX_CHANNELS];   // interleaved frames
    unsigned long cap_frames;
    unsigned long channels;
    volatile unsigned long head;   // write cursor (frames)
    volatile unsigned long tail;   // read cursor (frames)

    unsigned long readable_frames() const;
    unsigned long writable_frames() const;
    unsigned long write_frames([]const Fixed src);   // returns frames written
    unsigned long read_frames([]Fixed dst);    // returns frames read
    unsigned long peek_frames([]Fixed dst) const;
    unsigned long mix_frames([]const Fixed src);
    void          clear();
    int           is_empty() const;
    int           is_full() const;

    // 'channels'개 채널의 비어 있고 0으로 채운 버퍼 (AUDIO_BUF_FRAMES 프레임).
    static struct AudioBuffer new(unsigned long channels);
};

#include <std/dsp/audio_buffer.h>
#include <std/io.h>

static struct AudioBuffer g_ab;

int main() {
    g_ab = std::AudioBuffer::new(2UL);        // stereo

    Fixed frame[2];
    unsafe {
        frame[0] = std::fixed_from_float(0.5f);
        frame[1] = std::fixed_from_float(-0.25f);
    }
    g_ab.write_frames(frame[0..1]);           // producer
    std::print("readable: ");
    std::println_int((long long)g_ab.readable_frames());   // 1

    Fixed out[2];
    g_ab.read_frames(out[0..1]);              // consumer
    unsafe {
        std::print("L = "); std::println_float((double)std::fixed_to_float(out[0]));   //  0.5
        std::print("R = "); std::println_float((double)std::fixed_to_float(out[1]));   // -0.25
    }
    return 0;
}

#include <std/dsp/timer_wheel.h>

#define WHEEL_SLOTS   256
#define WHEEL_TIMERS  64

struct TimerWheel {
    struct WheelTimer timers[WHEEL_TIMERS];
    unsigned long     current_tick;

    void          init();            // in place: all slots free, tick 0. Call before anything else.

    // One-shot at absolute tick `tick`. Returns the timer id (1-based) or 0 if full.
    generic<T> int add(unsigned long tick, fn void(?&T) cb, ?&T ctx);
    // Every `period` ticks, first at now() + period. Returns the id or 0 if full.
    generic<T> int add_periodic(unsigned long period, fn void(?&T) cb, ?&T ctx);

    void          cancel(int id);    // ignores invalid ids
    int           tick();            // advance one tick; returns how many timers fired
    unsigned long now() const;
};

#include <std/dsp/timer_wheel.h>
#include <std/io.h>

struct Led { int toggles; };
static int g_done = 0;

static void blink(?&Led led) {
    match (led) { case some(l): l.toggles = l.toggles + 1; case null: { } }
}
static void stop(?&int unused) { g_done = 1; }

static struct TimerWheel g_tw;   // ~2 KB of slots: keep it off small stacks

int main() {
    g_tw.init();
    struct Led led;
    led.toggles = 0;
    unsafe { g_tw.add_periodic(10UL, blink, (?&Led)&led); }  // every 10 ticks
    g_tw.add(55UL, stop, null);                               // once, at tick 55
    while (!g_done) { g_tw.tick(); }
    std::print("blinks: ");
    std::println_int(led.toggles);   // 5 (ticks 10, 20, 30, 40, 50)
    return 0;
}