minidsp.h File Reference

A mini library of DSP (Digital Signal Processing) routines. More...

#include <string.h>
#include <assert.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <complex.h>
#include <fftw3.h>

Include dependency graph for minidsp.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	M_PI   3.14159265358979323846

Enumerations
enum	MD_GCC_WEIGHTING_TYPE { SIMP , PHAT }
	Weighting types for Generalized Cross-Correlation. More...

Functions
double	MD_dot (const double *a, const double *b, unsigned N)
	Compute the dot product of two vectors.
double	MD_entropy (const double *a, unsigned N, bool clip)
	Compute the normalized entropy of a distribution.
double	MD_energy (const double *a, unsigned N)
	Compute signal energy: sum of squared samples.
double	MD_power (const double *a, unsigned N)
	Compute signal power: energy divided by the number of samples.
double	MD_power_db (const double *a, unsigned N)
	Compute signal power in decibels: 10 * log10(power).
double	MD_scale (double in, double oldmin, double oldmax, double newmin, double newmax)
	Map a single value from one range to another.
void	MD_scale_vec (double *in, double *out, unsigned N, double oldmin, double oldmax, double newmin, double newmax)
	Map every element of a vector from one range to another.
void	MD_fit_within_range (double *in, double *out, unsigned N, double newmin, double newmax)
	Fit values within [newmin, newmax].
void	MD_adjust_dblevel (const double *in, double *out, unsigned N, double dblevel)
	Automatic Gain Control: scale a signal so its power matches the requested dB level, then clip to [-1, 1].
void	MD_magnitude_spectrum (const double *signal, unsigned N, double *mag_out)
	Compute the magnitude spectrum of a real-valued signal.
void	MD_power_spectral_density (const double *signal, unsigned N, double *psd_out)
	Compute the power spectral density (PSD) of a real-valued signal.
void	MD_phase_spectrum (const double *signal, unsigned N, double *phase_out)
	Compute the one-sided phase spectrum of a real signal.
unsigned	MD_stft_num_frames (unsigned signal_len, unsigned N, unsigned hop)
	Compute the number of STFT frames for the given signal length and parameters.
void	MD_stft (const double *signal, unsigned signal_len, unsigned N, unsigned hop, double *mag_out)
	Compute the Short-Time Fourier Transform (STFT) of a real-valued signal.
void	MD_Gen_Hann_Win (double *out, unsigned n)
	Generate a Hanning window of length n.
void	MD_sine_wave (double *output, unsigned N, double amplitude, double freq, double sample_rate)
	Generate a sine wave.
void	MD_white_noise (double *output, unsigned N, double amplitude, unsigned seed)
	Generate Gaussian white noise.
void	MD_impulse (double *output, unsigned N, double amplitude, unsigned position)
	Generate a discrete impulse (Kronecker delta).
void	MD_chirp_linear (double *output, unsigned N, double amplitude, double f_start, double f_end, double sample_rate)
	Generate a linear chirp (swept sine with linearly increasing frequency).
void	MD_chirp_log (double *output, unsigned N, double amplitude, double f_start, double f_end, double sample_rate)
	Generate a logarithmic chirp (swept sine with exponentially increasing frequency).
void	MD_square_wave (double *output, unsigned N, double amplitude, double freq, double sample_rate)
	Generate a square wave.
void	MD_sawtooth_wave (double *output, unsigned N, double amplitude, double freq, double sample_rate)
	Generate a sawtooth wave.
void	MD_shutdown (void)
	Free all internally cached FFT plans and buffers.
void	MD_get_multiple_delays (const double **sigs, unsigned M, unsigned N, unsigned margin, int weightfunc, int *outdelays)
	Estimate delays between a reference signal and M-1 other signals.
int	MD_get_delay (const double *siga, const double *sigb, unsigned N, double *ent, unsigned margin, int weightfunc)
	Estimate the delay between two signals.
void	MD_gcc (const double *siga, const double *sigb, unsigned N, double *lagvals, int weightfunc)
	Compute the full generalized cross-correlation between two signals.

Detailed Description

A mini library of DSP (Digital Signal Processing) routines.

This header declares functions for:

• Basic signal measurements (energy, power, entropy)
• Signal scaling and gain adjustment
• Window generation (Hanning window)
• Signal generators (sine, white noise, impulse, chirp, square, sawtooth)
• FFT-based magnitude spectrum, power spectral density, and STFT
• Generalized Cross-Correlation (GCC-PHAT) for delay estimation

These are the kinds of building blocks you'd use in an audio processing pipeline – for example, estimating which direction a sound came from using a pair of microphones.

Definition in file minidsp.h.

Macro Definition Documentation

M_PI

#define M_PI   3.14159265358979323846

Definition at line 27 of file minidsp.h.

Enumeration Type Documentation

MD_GCC_WEIGHTING_TYPE

enum MD_GCC_WEIGHTING_TYPE

Weighting types for Generalized Cross-Correlation.

Enumerator
SIMP	Simple 1/N weighting (basic cross-correlation)
PHAT	Phase Transform weighting (sharper peaks, more robust to noise)

Definition at line 563 of file minidsp.h.

Function Documentation

MD_adjust_dblevel()

void MD_adjust_dblevel	(	const double *	in,
		double *	out,
		unsigned	N,
		double	dblevel
	)

Automatic Gain Control: scale a signal so its power matches the requested dB level, then clip to [-1, 1].

Automatic Gain Control: scale a signal so its power matches the requested dB level, then clip to [-1, 1].

The input signal is assumed to be in the range [-1.0, 1.0]. The function computes a gain factor so that the output signal has the desired power level (in dB), then clips to [-1.0, 1.0] if any sample exceeds that range.

Based on Jaydeep Dhole's AGC implementation for MATLAB.

Definition at line 438 of file minidsp.c.

MD_chirp_linear()

void MD_chirp_linear	(	double *	output,
		unsigned	N,
		double	amplitude,
		double	f_start,
		double	f_end,
		double	sample_rate
	)

Generate a linear chirp (swept sine with linearly increasing frequency).

The instantaneous frequency sweeps linearly from f_start to f_end over N samples:

f(t) = f_start + (f_end - f_start) * t / T

where T = (N-1) / sample_rate is the sweep duration. The output is:

output[i] = amplitude * sin(2*pi * (f_start*t + 0.5*chirp_rate*t^2))

A linear chirp is the standard test signal for spectrograms – its instantaneous frequency traces a straight diagonal line in the time-frequency plane.

Parameters

output	Output buffer of length N (caller-allocated).
N	Number of samples to generate. Must be > 0.
amplitude	Peak amplitude of the chirp (e.g. 1.0).
f_start	Starting frequency in Hz.
f_end	Ending frequency in Hz.
sample_rate	Sampling rate in Hz. Must be > 0.

double sig[16000];
// 1-second linear chirp from 200 Hz to 4000 Hz at 16 kHz sample rate
MD_chirp_linear(sig, 16000, 1.0, 200.0, 4000.0, 16000.0);

Definition at line 617 of file minidsp.c.

MD_chirp_log()

void MD_chirp_log	(	double *	output,
		unsigned	N,
		double	amplitude,
		double	f_start,
		double	f_end,
		double	sample_rate
	)

Generate a logarithmic chirp (swept sine with exponentially increasing frequency).

The instantaneous frequency sweeps exponentially from f_start to f_end over N samples:

f(t) = f_start * (f_end / f_start)^(t / T)

where T = (N-1) / sample_rate is the sweep duration. The output is:

output[i] = amplitude * sin(2*pi * f_start * T * (r^(t/T) - 1) / ln(r))

where r = f_end / f_start.

A logarithmic chirp spends equal time per octave, making it ideal for measuring systems whose behaviour is best described on a log-frequency axis (e.g. audio equaliser response).

Parameters

output	Output buffer of length N (caller-allocated).
N	Number of samples to generate. Must be > 0.
amplitude	Peak amplitude of the chirp (e.g. 1.0).
f_start	Starting frequency in Hz. Must be > 0.
f_end	Ending frequency in Hz. Must be > 0 and != f_start.
sample_rate	Sampling rate in Hz. Must be > 0.

double sig[44100];
// 1-second log chirp from 20 Hz to 20 kHz at 44.1 kHz sample rate
MD_chirp_log(sig, 44100, 1.0, 20.0, 20000.0, 44100.0);

Definition at line 639 of file minidsp.c.

MD_dot()

double MD_dot	(	const double *	a,
		const double *	b,
		unsigned	N
	)

Compute the dot product of two vectors.

The dot product is the sum of element-wise products: a[0]*b[0] + a[1]*b[1] + ...

Compute the dot product of two vectors.

dot = a[0]*b[0] + a[1]*b[1] + ... + a[N-1]*b[N-1]

This is one of the most fundamental operations in DSP. The dot product measures how "similar" two signals are – if they are identical the result is large; if they are unrelated it is near zero.

Definition at line 298 of file minidsp.c.

MD_energy()

double MD_energy	(	const double *	a,
		unsigned	N
	)

Compute signal energy: sum of squared samples.

Compute signal energy: sum of squared samples.

     N-1

E = SUM a[n]^2 n=0

Energy tells you "how loud" a signal is overall. A silent signal has zero energy; a loud one has high energy. Notice that squaring makes all values positive, so negative samples contribute too.

Definition at line 318 of file minidsp.c.

MD_entropy()

double MD_entropy	(	const double *	a,
		unsigned	N,
		bool	clip
	)

Compute the normalized entropy of a distribution.

Returns a value between 0.0 (all energy concentrated in one bin) and 1.0 (energy spread equally across all bins).

Parameters

clip	If true, ignore negative values. If false, square all values first.

Compute the normalized entropy of a distribution.

Entropy measures how "spread out" a distribution is:

• 0.0 means all the energy is in a single bin (very peaky).
• 1.0 means the energy is spread equally (completely flat).

The Shannon entropy formula is: H = -SUM( p_i * log2(p_i) ) We normalise by dividing by log2(N) so the result is always in [0, 1].

Parameters

a	Array of values representing the distribution.
N	Length of the array.
clip	If true, ignore negative values (treat them as 0). If false, use a[i]^2 so all values become non-negative.

Returns

Normalised entropy in [0.0, 1.0].

Definition at line 483 of file minidsp.c.

MD_fit_within_range()

void MD_fit_within_range	(	double *	in,
		double *	out,
		unsigned	N,
		double	newmin,
		double	newmax
	)

Fit values within [newmin, newmax].

If all values already fit, they are copied unchanged. Otherwise the entire vector is rescaled.

Fit values within [newmin, newmax].

If the input range already lies inside [newmin, newmax], the values are copied as-is (no stretching). Otherwise, the full input range is mapped onto the new range.

Definition at line 402 of file minidsp.c.

MD_gcc()

void MD_gcc	(	const double *	siga,
		const double *	sigb,
		unsigned	N,
		double *	lagvals,
		int	weightfunc
	)

Compute the full generalized cross-correlation between two signals.

Parameters

siga	First signal.
sigb	Second signal.
N	Length of both signals.
lagvals	Output array of N doubles (pre-allocated). The zero-lag value is at index ceil(N/2).
weightfunc	SIMP or PHAT.

Compute the full generalized cross-correlation between two signals.

This is the core GCC-PHAT routine. It fills lagvals[] with the cross-correlation values, shifted so that the zero-lag value is at index ceil(N/2).

Parameters

siga	First signal.
sigb	Second signal.
N	Length of both signals.
lagvals	Output array of N doubles (pre-allocated by caller).
weightfunc	SIMP or PHAT.

Definition at line 1028 of file minidsp.c.

MD_Gen_Hann_Win()

void MD_Gen_Hann_Win	(	double *	out,
		unsigned	n
	)

Generate a Hanning window of length n.

A Hanning window tapers the edges of a signal to zero, which reduces spectral leakage when you later take an FFT.

Generate a Hanning window of length n.

The formula is: w[i] = 0.5 * (1 - cos(2*pi*i / (n-1)))

Windowing is essential before taking an FFT because real-world signals don't start and end at exactly zero. Without a window, the abrupt edges create false high-frequency content ("spectral leakage"). The Hanning window smoothly tapers the signal to zero at both ends, reducing this artefact.

Definition at line 543 of file minidsp.c.

MD_get_delay()

int MD_get_delay	(	const double *	siga,
		const double *	sigb,
		unsigned	N,
		double *	ent,
		unsigned	margin,
		int	weightfunc
	)

Estimate the delay between two signals.

Parameters

siga	First signal.
sigb	Second signal.
N	Length of both signals.
ent	If non-null, receives the normalised entropy of the correlation peak region (closer to 1.0 = less trustworthy).
margin	Search +/- this many samples around zero-lag.
weightfunc	SIMP or PHAT.

Returns

Delay in samples (positive = sigb lags siga).

Estimate the delay between two signals.

The function computes the GCC between the two signals, then searches within a +/- margin window around zero-lag for the peak. The offset of that peak from zero is the estimated delay.

Parameters

siga	First signal (reference).
sigb	Second signal.
N	Length of both signals.
ent	If non-null, receives the normalised entropy of the lag values in the search window. High entropy (~1.0) means a flat, unreliable correlation; low entropy (~0.0) means a sharp, trustworthy peak.
margin	How far (in samples) to search around zero-lag.
weightfunc	SIMP or PHAT.

Returns

Delay in samples.

Definition at line 978 of file minidsp.c.

MD_get_multiple_delays()

void MD_get_multiple_delays	(	const double **	sigs,
		unsigned	M,
		unsigned	N,
		unsigned	margin,
		int	weightfunc,
		int *	outdelays
	)

Estimate delays between a reference signal and M-1 other signals.

Parameters

sigs	Array of M pointers to signals (sigs[0] is the reference).
M	Number of signals.
N	Length of each signal (all must be the same length).
margin	Search +/- this many samples around zero-lag.
weightfunc	SIMP or PHAT (see MD_GCC_WEIGHTING_TYPE).
outdelays	Output array of M-1 delay values (must be pre-allocated).

Estimate delays between a reference signal and M-1 other signals.

This is a convenience wrapper that calls MD_get_delay() for each non-reference signal. A Hanning window is applied to all signals before computing the cross-correlation.

Parameters

sigs	Array of M pointers to signals. sigs[0] is the reference.
M	Number of signals (must be >= 2).
N	Length of every signal.
margin	Search +/- margin samples around zero-lag.
weightfunc	SIMP or PHAT.
outdelays	Output: M-1 delay values (pre-allocated by caller).

Definition at line 913 of file minidsp.c.

MD_impulse()

void MD_impulse	(	double *	output,
		unsigned	N,
		double	amplitude,
		unsigned	position
	)

Generate a discrete impulse (Kronecker delta).

Fills the output buffer with zeros except at position, where the value is set to amplitude. A unit impulse (amplitude 1.0 at position 0) is the identity element of convolution and has a perfectly flat magnitude spectrum.

Common uses:

• Measure a system's impulse response by feeding it through a filter.
• Verify that MD_magnitude_spectrum() returns a flat spectrum.
• Create delayed spikes for testing convolution and delay estimation.

Parameters

output	Output buffer of length N (caller-allocated).
N	Number of samples to generate. Must be > 0.
amplitude	Value of the impulse spike (e.g. 1.0 for unit impulse).
position	Sample index of the spike (0-based). Must be < N.

double sig[1024];
MD_impulse(sig, 1024, 1.0, 0);  // unit impulse at sample 0

Definition at line 608 of file minidsp.c.

MD_magnitude_spectrum()

void MD_magnitude_spectrum	(	const double *	signal,
		unsigned	N,
		double *	mag_out
	)

Compute the magnitude spectrum of a real-valued signal.

Given a signal of length N, this function computes the FFT and returns the magnitude |X(k)| for each frequency bin. Because the input is real-valued, the FFT output is conjugate-symmetric, so only the first N/2 + 1 bins are unique (from DC through Nyquist).

To convert bin index k to a frequency in Hz: freq_k = k * sample_rate / N

The output is not normalised by N – divide each value by N to get the "standard" DFT magnitude, or by N/2 (except DC and Nyquist) for single-sided amplitude.

Parameters

signal	Input signal of length N.
N	Number of samples in the signal.
mag_out	Output array, must be pre-allocated to at least N/2 + 1 doubles. On return, mag_out[k] = |X(k)|.

Note

The caller must allocate mag_out. The required size is (N / 2 + 1) * sizeof(double).

Example:

double signal[1024];
// ... fill signal with audio samples ...
 
unsigned num_bins = 1024 / 2 + 1;  // = 513
double *mag = malloc(num_bins * sizeof(double));
MD_magnitude_spectrum(signal, 1024, mag);
 
// mag[0]   = DC component magnitude
// mag[k]   = magnitude at frequency k * sample_rate / 1024
// mag[512] = Nyquist frequency magnitude
 
free(mag);
MD_shutdown();  // free cached FFT plans when done

This performs a real-to-complex FFT using FFTW, then computes the absolute value (magnitude) of each complex frequency bin.

For a real signal of length N, the FFT is conjugate-symmetric, so only the first N/2 + 1 bins are unique:

• Bin 0: DC component (zero frequency)
• Bin k: frequency = k * sample_rate / N
• Bin N/2: Nyquist frequency (sample_rate / 2)

The magnitude is computed as: |X(k)| = sqrt( Re(X(k))^2 + Im(X(k))^2 )

The FFT plan is cached and reused across calls of the same length, following the same pattern as the GCC functions.

Parameters

signal	Input signal of length N.
N	Number of samples (must be >= 2).
mag_out	Output: magnitudes for bins 0..N/2. Must be pre-allocated to N/2 + 1 doubles.

Definition at line 727 of file minidsp.c.

MD_phase_spectrum()

void MD_phase_spectrum	(	const double *	signal,
		unsigned	N,
		double *	phase_out
	)

Compute the one-sided phase spectrum of a real signal.

Returns the instantaneous phase angle \(\phi(k) = \arg X(k)\) for each DFT bin using an unnormalised real-to-complex FFT (FFTW r2c). The phase is expressed in radians in the range \([-\pi,\,\pi]\).

For a real signal of length \(N\), only the first \(N/2+1\) bins carry unique information (bins \(N/2+1\ldots N-1\) are conjugate-symmetric mirrors). Accordingly, phase_out must be pre-allocated to hold at least \(N/2+1\) doubles.

Interpretation:

• A pure cosine at an integer bin \(k_0\) (exact period in \(N\) samples) produces \(\phi(k_0) = 0\).
• A pure sine at the same bin produces \(\phi(k_0) = -\pi/2\).
• A time-delayed signal exhibits linear phase: \(\phi(k) = -2\pi k d/N\), where \(d\) is the delay in samples.

Note

Phase values at bins where the magnitude is near zero are numerically unreliable. Always examine MD_magnitude_spectrum() alongside the phase spectrum to identify significant bins.

Parameters

[in]	signal	Input signal of length N.
[in]	N	Signal length (must be >= 2).
[out]	phase_out	Pre-allocated array of at least N/2+1 doubles that receives the phase in radians.

Example

unsigned N = 1024;
double *sig = malloc(N * sizeof(double));
// ... fill sig ...
unsigned num_bins = N / 2 + 1;
double *phase = malloc(num_bins * sizeof(double));
MD_phase_spectrum(sig, N, phase);
// phase[k] in [-M_PI, M_PI]
free(sig);
free(phase);
MD_shutdown();

See also

MD_magnitude_spectrum(), MD_power_spectral_density(), MD_shutdown()

Compute the one-sided phase spectrum of a real signal.

The phase is the argument (angle) of each complex DFT coefficient: phi(k) = atan2(Im(X(k)), Re(X(k)))

This reuses the same FFT cache as MD_magnitude_spectrum() and MD_power_spectral_density() – no additional plan allocation.

Phase is scale-invariant (multiplying a signal by a positive constant does not change its phase), so no normalisation by N is needed.

Parameters

signal	Input signal of length N.
N	Number of samples (must be >= 2).
phase_out	Output: phase in radians for bins 0..N/2. Must be pre-allocated to N/2 + 1 doubles.

Definition at line 807 of file minidsp.c.

MD_power()

double MD_power	(	const double *	a,
		unsigned	N
	)

Compute signal power: energy divided by the number of samples.

Compute signal power: energy divided by the number of samples.

P = E / N

Power is the "average energy per sample". It lets you compare signals of different lengths on an equal footing.

Definition at line 333 of file minidsp.c.

MD_power_db()

double MD_power_db	(	const double *	a,
		unsigned	N
	)

Compute signal power in decibels: 10 * log10(power).

Compute signal power in decibels: 10 * log10(power).

P_dB = 10 * log10(P)

Decibels are a logarithmic scale commonly used in audio engineering. Every +10 dB means the power is 10x larger. A floor of 1e-10 is used to avoid log(0), which would be negative infinity.

Definition at line 349 of file minidsp.c.

MD_power_spectral_density()

void MD_power_spectral_density	(	const double *	signal,
		unsigned	N,
		double *	psd_out
	)

Compute the power spectral density (PSD) of a real-valued signal.

The PSD describes how a signal's power is distributed across frequencies. While the magnitude spectrum tells you the amplitude at each frequency, the PSD tells you the power – useful for noise analysis, SNR estimation, and comparing signals of different lengths.

This function computes the "periodogram" estimator:

PSD[k] = |X(k)|^2 / N = (Re(X(k))^2 + Im(X(k))^2) / N

where X(k) is the DFT of the input signal (unnormalised, as computed by FFTW).

Relationship to the magnitude spectrum: PSD[k] = |X(k)|^2 / N = (magnitude[k])^2 / N

Parseval's theorem (energy conservation): The one-sided PSD sums to the total signal energy: PSD[0] + 2 * sum(PSD[1..N/2-1]) + PSD[N/2] = sum(x[n]^2)

Parameters

signal	Input signal of length N.
N	Number of samples in the signal (must be >= 2).
psd_out	Output array, must be pre-allocated to at least N/2 + 1 doubles. On return, psd_out[k] = |X(k)|^2 / N.

Note

The caller must allocate psd_out. The required size is (N / 2 + 1) * sizeof(double).

Example:

double signal[1024];
// ... fill signal with audio samples ...
 
unsigned num_bins = 1024 / 2 + 1;  // = 513
double *psd = malloc(num_bins * sizeof(double));
MD_power_spectral_density(signal, 1024, psd);
 
// psd[0]   = DC power
// psd[k]   = power at frequency k * sample_rate / 1024
// psd[512] = Nyquist frequency power
 
free(psd);
MD_shutdown();  // free cached FFT plans when done

The PSD is the "periodogram" estimator: PSD[k] = |X(k)|^2 / N. It reuses the same FFT cache as MD_magnitude_spectrum() – both perform the same real-to-complex FFT, only the post-processing differs.

We compute |X(k)|^2 = re^2 + im^2 directly from the real and imaginary parts, rather than calling cabs() (which computes sqrt(re^2 + im^2)) and then squaring. This avoids a redundant sqrt and is both faster and more numerically precise.

Parameters

signal	Input signal of length N.
N	Number of samples (must be >= 2).
psd_out	Output: PSD for bins 0..N/2. Must be pre-allocated to N/2 + 1 doubles.

Definition at line 765 of file minidsp.c.

MD_sawtooth_wave()

void MD_sawtooth_wave	(	double *	output,
		unsigned	N,
		double	amplitude,
		double	freq,
		double	sample_rate
	)

Generate a sawtooth wave.

Fills the output buffer with a sawtooth wave that ramps linearly from −amplitude to +amplitude over each period:

\[ x[n] = A \left(\frac{\phi}{\pi} - 1\right) \]

where \(\phi = 2\pi f n / f_s \pmod{2\pi}\).

A sawtooth wave contains all integer harmonics (1f, 2f, 3f, …) whose amplitudes decay as 1/k — useful for demonstrating richer harmonic content compared to the square wave's odd-only series.

Parameters

output	Output buffer of length N (caller-allocated).
N	Number of samples to generate. Must be > 0.
amplitude	Peak amplitude of the sawtooth wave (e.g. 1.0).
freq	Frequency in Hz.
sample_rate	Sampling rate in Hz. Must be > 0.

double sig[1024];
MD_sawtooth_wave(sig, 1024, 1.0, 440.0, 44100.0);

Definition at line 685 of file minidsp.c.

MD_scale()

double MD_scale	(	double	in,
		double	oldmin,
		double	oldmax,
		double	newmin,
		double	newmax
	)

Map a single value from one range to another.

Example: MD_scale(5, 0, 10, 0, 100) returns 50.

Map a single value from one range to another.

Formula: out = (in - oldmin) * (newmax - newmin) / (oldmax - oldmin) + newmin

For example, mapping a value of 5 from the range [0, 10] into [0, 100] gives 50. This is sometimes called "lerp" (linear interpolation).

Definition at line 369 of file minidsp.c.

MD_scale_vec()

void MD_scale_vec	(	double *	in,
		double *	out,
		unsigned	N,
		double	oldmin,
		double	oldmax,
		double	newmin,
		double	newmax
	)

Map every element of a vector from one range to another.

Map every element of a vector from one range to another.

Definition at line 379 of file minidsp.c.

MD_shutdown()

void MD_shutdown

(

void

)

Free all internally cached FFT plans and buffers.

Call when done.

Definition at line 276 of file minidsp.c.

MD_sine_wave()

void MD_sine_wave	(	double *	output,
		unsigned	N,
		double	amplitude,
		double	freq,
		double	sample_rate
	)

Generate a sine wave.

Fills output[i] = amplitude * sin(2π * freq * i / sample_rate) for i in [0, N).

This is the simplest test signal in DSP — a pure tone at a single frequency. Use it to verify filter responses, check FFT bin alignment, or provide a clean input for any processing chain.

Parameters

output	Output buffer of length N (caller-allocated).
N	Number of samples to generate. Must be > 0.
amplitude	Peak amplitude of the sine wave (e.g. 1.0).
freq	Frequency in Hz.
sample_rate	Sampling rate in Hz. Must be > 0.

double sig[1024];
MD_sine_wave(sig, 1024, 1.0, 440.0, 44100.0);  // 440 Hz A4 tone

Definition at line 555 of file minidsp.c.

MD_square_wave()

void MD_square_wave	(	double *	output,
		unsigned	N,
		double	amplitude,
		double	freq,
		double	sample_rate
	)

Generate a square wave.

Fills the output buffer with a bipolar square wave that alternates between +amplitude and −amplitude at the given frequency:

\[ x[n] = \begin{cases} +A & 0 < \phi < \pi \\ -A & \pi < \phi < 2\pi \\ 0 & \phi = 0 \text{ or } \phi = \pi \end{cases} \]

where \(\phi = 2\pi f n / f_s \pmod{2\pi}\).

A square wave contains only odd harmonics (1f, 3f, 5f, …) whose amplitudes decay as 1/k — a classic demonstration of the Fourier series and the Gibbs phenomenon.

Parameters

output	Output buffer of length N (caller-allocated).
N	Number of samples to generate. Must be > 0.
amplitude	Peak amplitude of the square wave (e.g. 1.0).
freq	Frequency in Hz.
sample_rate	Sampling rate in Hz. Must be > 0.

double sig[1024];
MD_square_wave(sig, 1024, 1.0, 440.0, 44100.0);

Definition at line 666 of file minidsp.c.

MD_stft()

void MD_stft	(	const double *	signal,
		unsigned	signal_len,
		unsigned	N,
		unsigned	hop,
		double *	mag_out
	)

Compute the Short-Time Fourier Transform (STFT) of a real-valued signal.

The STFT slides a Hanning-windowed FFT over the signal in steps of hop samples, producing a time-frequency magnitude matrix.

For each frame f starting at sample f * hop, the function:

1. Multiplies the frame by a Hanning window.
2. Computes the real-to-complex FFT using FFTW.
3. Stores the magnitudes |X(k)| for bins 0..N/2 in the output.

The STFT formula for frame f and bin k is:

X_f(k) = SUM_{n=0}^{N-1} w[n] * x[f*hop + n] * e^{-j2pi*k*n/N}

mag_out[f * (N/2+1) + k] = |X_f(k)|

where w[n] is the Hanning window.

The output is not normalised by N – divide each value by N to get the "standard" DFT magnitude, consistent with MD_magnitude_spectrum().

The STFT reuses the same cached FFT plan as MD_magnitude_spectrum() and MD_power_spectral_density(). Only the Hanning window buffer is separate.

Parameters

signal	Input signal.
signal_len	Total number of samples in the signal.
N	FFT window size (must be >= 2).
hop	Hop size in samples (must be >= 1).
mag_out	Output array (row-major). Must be pre-allocated to at least MD_stft_num_frames(signal_len, N, hop) * (N/2+1) doubles. On return, mag_out[f*(N/2+1) + k] = |X_f(k)|.

Note

If signal_len < N, the function returns immediately without writing any output (zero frames fit). Use MD_stft_num_frames() to check in advance.

The caller must allocate mag_out. A typical pattern:

unsigned N   = 512;
unsigned hop = 128;
unsigned num_frames = MD_stft_num_frames(signal_len, N, hop);
unsigned num_bins   = N / 2 + 1;
 
double *mag_out = malloc(num_frames * num_bins * sizeof(double));
MD_stft(signal, signal_len, N, hop, mag_out);
 
// mag_out[f * num_bins + k] = |X_f(k)|
// Convert bin k to Hz: freq_hz = k * sample_rate / N
// Convert frame f to seconds: time_s = (double)(f * hop) / sample_rate
 
free(mag_out);
MD_shutdown();  // free cached FFT plans when done

Compute the Short-Time Fourier Transform (STFT) of a real-valued signal.

Slides a Hanning-windowed r2c FFT over the signal in steps of hop samples. For each frame, the window is applied by multiplying directly into the shared _spec_in buffer (no separate memcpy pass), then FFTW executes the plan and magnitudes |X(k)| are written to the output row.

The shared spec* plan is reused so that interleaving calls to MD_stft(), MD_magnitude_spectrum(), and MD_power_spectral_density() with the same N incurs no extra plan-rebuild overhead.

Parameters

signal	Input signal.
signal_len	Number of samples (0 frames if < N, see header).
N	FFT window size (>= 2).
hop	Hop size (>= 1).
mag_out	Row-major output: mag_out[f*(N/2+1) + k] = |X_f(k)|. Must be pre-allocated by caller.

Definition at line 860 of file minidsp.c.

MD_stft_num_frames()

unsigned MD_stft_num_frames	(	unsigned	signal_len,
		unsigned	N,
		unsigned	hop
	)

Compute the number of STFT frames for the given signal length and parameters.

The formula is: num_frames = (signal_len - N) / hop + 1 when signal_len >= N num_frames = 0 when signal_len < N

Use this function to size the output buffer before calling MD_stft().

Parameters

signal_len	Total number of samples in the signal.
N	FFT window size (samples per frame).
hop	Hop size (samples between successive frame starts).

Returns

Number of complete frames that fit in the signal.

Example:

// 1 s of audio at 16 kHz, 32 ms window, 8 ms hop
unsigned signal_len = 16000;
unsigned N   = 512;   // 32 ms at 16 kHz
unsigned hop = 128;   // 8 ms (75% overlap)
unsigned num_frames = MD_stft_num_frames(signal_len, N, hop);
// num_frames = (16000 - 512) / 128 + 1 = 121

Compute the number of STFT frames for the given signal length and parameters.

Returns (signal_len - N) / hop + 1 when signal_len >= N, else 0. Integer division truncates, so only complete frames are counted.

Definition at line 835 of file minidsp.c.

MD_white_noise()

void MD_white_noise	(	double *	output,
		unsigned	N,
		double	amplitude,
		unsigned	seed
	)

Generate Gaussian white noise.

Fills output with normally distributed random samples (mean 0, standard deviation amplitude) using the Box-Muller transform. White noise has equal energy at all frequencies – its power spectral density is approximately flat.

Use white noise to test filters, measure impulse responses, or as an additive noise source for SNR experiments.

Parameters

output	Output buffer of length N (caller-allocated).
N	Number of samples to generate. Must be > 0.
amplitude	Standard deviation of the noise (e.g. 1.0).
seed	Seed for the random number generator. Using the same seed produces the same output sequence, which is useful for reproducible tests.

double noise[4096];
MD_white_noise(noise, 4096, 1.0, 42);  // reproducible Gaussian noise

Definition at line 566 of file minidsp.c.