FFT-based spectrum analysis: magnitude spectrum, PSD, phase spectrum, STFT, mel filterbanks, MFCCs, FFT-based F0 estimation, and MD_shutdown(). More...

#include "minidsp.h"
#include "minidsp_internal.h"

Include dependency graph for minidsp_spectrum.c:

Go to the source code of this file.

Macros
#define	MD_F0_FFT_PROMINENCE_RATIO 2.5
	F0-FFT peak must stand out from the average in-range spectrum.

#define	MD_F0_FFT_GLOBAL_RATIO 0.2
	Candidate must also be a meaningful fraction of the frame's global peak.

#define	MD_MFCC_LOG_FLOOR 1e-12
	Floor before logarithm in MFCC log-mel compression.

Functions
static double	md_parabolic_offset (double y_left, double y_mid, double y_right)
	Three-point parabolic refinement around a discrete peak index.

static double	md_hz_to_mel (double hz)
	HTK mel scale conversion: Hz -> mel.

static double	md_mel_to_hz (double mel)
	HTK mel scale conversion: mel -> Hz.

static void	_spec_free (void)
	Free spectrum analysis buffers.

static void	_spec_teardown (void)
	Tear down spectrum analysis plan and buffers.

static void	_spec_setup (unsigned N)
	Ensure the spectrum FFT cache matches the requested length.

static void	_stft_teardown (void)
	Free the cached STFT Hanning window.

static void	_mel_teardown (void)
	Free mel filterbank cache.

static void	_mel_setup (unsigned N, double sample_rate, unsigned num_mels, double min_freq_hz, double max_freq_hz)
	Build the cached mel filterbank for the requested parameters.

static void	_stft_setup (unsigned N)
	Ensure the shared r2c plan is ready for length N, and that the STFT Hanning window buffer matches N.

void	MD_shutdown (void)
	Free all internally cached FFT plans and buffers.

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-valued signal.

double	MD_f0_fft (const double *signal, unsigned N, double sample_rate, double min_freq_hz, double max_freq_hz)
	Estimate the fundamental frequency (F0) using FFT peak picking.

unsigned	MD_stft_num_frames (unsigned signal_len, unsigned N, unsigned hop)
	Compute the number of complete STFT frames for the given parameters.

void	MD_stft (const double signal, unsigned signal_len, unsigned N, unsigned hop, double mag_out)
	Compute the Short-Time Fourier Transform (STFT) magnitude matrix.

void	MD_mel_filterbank (unsigned N, double sample_rate, unsigned num_mels, double min_freq_hz, double max_freq_hz, double *filterbank_out)
	Build a mel-spaced triangular filterbank matrix.

void	MD_mel_energies (const double signal, unsigned N, double sample_rate, unsigned num_mels, double min_freq_hz, double max_freq_hz, double mel_out)
	Compute mel-band energies from a single frame.

void	MD_mfcc (const double signal, unsigned N, double sample_rate, unsigned num_mels, unsigned num_coeffs, double min_freq_hz, double max_freq_hz, double mfcc_out)
	Compute mel-frequency cepstral coefficients (MFCCs) from a single frame.

Variables
static unsigned	_spec_N = 0

static double *	_spec_in = nullptr

static fftw_complex *	_spec_out = nullptr

static fftw_plan	_spec_plan = nullptr

static double *	_stft_win = nullptr

static unsigned	_stft_win_N = 0

static unsigned	_mel_N = 0

static unsigned	_mel_num_mels = 0

static double	_mel_sample_rate = 0.0

static double	_mel_min_freq_hz = 0.0

static double	_mel_max_freq_hz = 0.0

static double *	_mel_fb = nullptr

Detailed Description

FFT-based spectrum analysis: magnitude spectrum, PSD, phase spectrum, STFT, mel filterbanks, MFCCs, FFT-based F0 estimation, and MD_shutdown().

Author: Chuck Wooters woote.nosp@m.rs@h.nosp@m.ey.co.nosp@m.m

Copyright: 2013 International Computer Science Institute

Definition in file minidsp_spectrum.c.

Macro Definition Documentation

◆ MD_F0_FFT_GLOBAL_RATIO

#define MD_F0_FFT_GLOBAL_RATIO 0.2

Candidate must also be a meaningful fraction of the frame's global peak.

Definition at line 53 of file minidsp_spectrum.c.

◆ MD_F0_FFT_PROMINENCE_RATIO

#define MD_F0_FFT_PROMINENCE_RATIO 2.5

F0-FFT peak must stand out from the average in-range spectrum.

Definition at line 51 of file minidsp_spectrum.c.

◆ MD_MFCC_LOG_FLOOR

#define MD_MFCC_LOG_FLOOR 1e-12

Floor before logarithm in MFCC log-mel compression.

Definition at line 55 of file minidsp_spectrum.c.

Function Documentation

◆ _mel_setup()

static void _mel_setup	(	unsigned	N,
		double	sample_rate,
		unsigned	num_mels,
		double	min_freq_hz,
		double	max_freq_hz
	)

static

Build the cached mel filterbank for the requested parameters.

Definition at line 139 of file minidsp_spectrum.c.

◆ _mel_teardown()

static void _mel_teardown ( void )

static

Free mel filterbank cache.

Definition at line 127 of file minidsp_spectrum.c.

◆ _spec_free()

static void _spec_free ( void )

static

Free spectrum analysis buffers.

Definition at line 83 of file minidsp_spectrum.c.

◆ _spec_setup()

static void _spec_setup ( unsigned N )

static

Ensure the spectrum FFT cache matches the requested length.

Rebuilds the plan and buffers only when N changes.

Definition at line 103 of file minidsp_spectrum.c.

◆ _spec_teardown()

static void _spec_teardown ( void )

static

Tear down spectrum analysis plan and buffers.

Definition at line 92 of file minidsp_spectrum.c.

◆ _stft_setup()

static void _stft_setup ( unsigned N )

static

Ensure the shared r2c plan is ready for length N, and that the STFT Hanning window buffer matches N.

Rebuilds the window only when N changes.

Definition at line 222 of file minidsp_spectrum.c.

◆ _stft_teardown()

static void _stft_teardown ( void )

static

Free the cached STFT Hanning window.

Called only from MD_shutdown().

Definition at line 119 of file minidsp_spectrum.c.

◆ MD_f0_fft()

double MD_f0_fft	(	const double *	signal,
		unsigned	N,
		double	sample_rate,
		double	min_freq_hz,
		double	max_freq_hz
	)

Estimate the fundamental frequency (F0) using FFT peak picking.

The signal is internally Hann-windowed, transformed with a one-sided FFT, then the dominant magnitude peak in the requested frequency range is mapped back to Hz:

\[ f_0 = \frac{k_\text{peak} f_s}{N} \]

where \(k_\text{peak}\) is the selected FFT bin.

This method is simple and fast but generally less robust than MD_f0_autocorrelation() for noisy or strongly harmonic signals.

Parameters

signal	Input signal of length N.
N	Number of samples (must be >= 2).
sample_rate	Sampling rate in Hz (must be > 0).
min_freq_hz	Minimum search frequency in Hz (must be > 0).
max_freq_hz	Maximum search frequency in Hz (must be > min_freq_hz).

Returns: Estimated F0 in Hz, or 0.0 if no reliable F0 is found.

Note: This function uses the internal FFT cache. Call MD_shutdown() when done with FFT-based APIs.; Invalid API usage is guarded by assertions. A return value of 0.0 means the call was valid but no usable pitch peak was found.

double f0 = MD_f0_fft(frame, frame_len, 16000.0, 80.0, 400.0);

MD_f0_fft

double MD_f0_fft(const double *signal, unsigned N, double sample_rate, double min_freq_hz, double max_freq_hz)

Estimate the fundamental frequency (F0) using FFT peak picking.

Definition minidsp_spectrum.c:378

Definition at line 378 of file minidsp_spectrum.c.

◆ md_hz_to_mel()

static double md_hz_to_mel ( double hz )

static

HTK mel scale conversion: Hz -> mel.

Definition at line 71 of file minidsp_spectrum.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.

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 276 of file minidsp_spectrum.c.

◆ MD_mel_energies()

void MD_mel_energies	(	const double *	signal,
		unsigned	N,
		double	sample_rate,
		unsigned	num_mels,
		double	min_freq_hz,
		double	max_freq_hz,
		double *	mel_out
	)

Compute mel-band energies from a single frame.

Processing steps: 1) Apply an internal Hann window. 2) Compute one-sided PSD bins via FFT: |X(k)|^2 / N. 3) Apply mel filterbank weights and sum per band.

This function uses the same internal FFT cache as MD_stft() and spectrum APIs.

Parameters

signal	Input frame of length N.
N	Frame length / FFT size (must be >= 2).
sample_rate	Sampling rate in Hz (must be > 0).
num_mels	Number of mel bands (must be > 0).
min_freq_hz	Requested lower frequency bound in Hz.
max_freq_hz	Requested upper frequency bound in Hz (must be > min_freq_hz).
mel_out	Output mel energies, caller-allocated length num_mels.

Definition at line 531 of file minidsp_spectrum.c.

◆ MD_mel_filterbank()

void MD_mel_filterbank	(	unsigned	N,
		double	sample_rate,
		unsigned	num_mels,
		double	min_freq_hz,
		double	max_freq_hz,
		double *	filterbank_out
	)

Build a mel-spaced triangular filterbank matrix.

This function creates num_mels triangular filters over one-sided FFT bins (0..N/2), using the HTK mel mapping:

mel(f) = 2595 * log10(1 + f / 700)

The requested frequency range [min_freq_hz, max_freq_hz] is clamped at runtime to [0, sample_rate/2]. If the clamped range is empty, all output weights are zero.

Output layout is row-major: filterbank_out[m * (N/2 + 1) + k] where m is the mel band index and k is the FFT bin index.

Parameters

N	FFT size (must be >= 2).
sample_rate	Sampling rate in Hz (must be > 0).
num_mels	Number of mel filters (must be > 0).
min_freq_hz	Requested lower frequency bound in Hz.
max_freq_hz	Requested upper frequency bound in Hz (must be > min_freq_hz).
filterbank_out	Output matrix, caller-allocated, length num_mels * (N/2 + 1) doubles.

Definition at line 513 of file minidsp_spectrum.c.

◆ md_mel_to_hz()

static double md_mel_to_hz ( double mel )

static

HTK mel scale conversion: mel -> Hz.

Definition at line 77 of file minidsp_spectrum.c.

◆ MD_mfcc()

void MD_mfcc	(	const double *	signal,
		unsigned	N,
		double	sample_rate,
		unsigned	num_mels,
		unsigned	num_coeffs,
		double	min_freq_hz,
		double	max_freq_hz,
		double *	mfcc_out
	)

Compute mel-frequency cepstral coefficients (MFCCs) from a single frame.

Conventions used by this API:

HTK mel mapping for filter placement.
Internal Hann windowing and one-sided PSD mel energies.
Natural-log compression: log(max(E_mel, 1e-12)).
DCT-II with HTK-C0 profile: c0 uses sqrt(1/M) normalization, c_n (n>0) uses sqrt(2/M), where M is num_mels.
Coefficient 0 (C0) is written to mfcc_out[0].

Parameters

signal	Input frame of length N.
N	Frame length / FFT size (must be >= 2).
sample_rate	Sampling rate in Hz (must be > 0).
num_mels	Number of mel bands (must be > 0).
num_coeffs	Number of cepstral coefficients to output (must be in [1, num_mels]).
min_freq_hz	Requested lower frequency bound in Hz.
max_freq_hz	Requested upper frequency bound in Hz (must be > min_freq_hz).
mfcc_out	Output array, caller-allocated length num_coeffs.

Definition at line 567 of file minidsp_spectrum.c.

◆ md_parabolic_offset()

static double md_parabolic_offset	(	double	y_left,
		double	y_mid,
		double	y_right
	)

static

Three-point parabolic refinement around a discrete peak index.

Returns a fractional offset in [-0.5, 0.5].

Definition at line 59 of file minidsp_spectrum.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-valued signal.

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 356 of file minidsp_spectrum.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 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 314 of file minidsp_spectrum.c.

◆ MD_shutdown()

void MD_shutdown ( void )

Free all internally cached FFT plans and buffers.

Call when done.

Definition at line 238 of file minidsp_spectrum.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) magnitude matrix.

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 478 of file minidsp_spectrum.c.

◆ MD_stft_num_frames()

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

Compute the number of complete STFT frames for the given parameters.

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 453 of file minidsp_spectrum.c.

Variable Documentation

◆ _mel_fb

double* _mel_fb = nullptr

static

Definition at line 48 of file minidsp_spectrum.c.

◆ _mel_max_freq_hz

double _mel_max_freq_hz = 0.0

static

Definition at line 47 of file minidsp_spectrum.c.

◆ _mel_min_freq_hz

double _mel_min_freq_hz = 0.0

static

Definition at line 46 of file minidsp_spectrum.c.

◆ _mel_N

unsigned _mel_N = 0

static

Definition at line 43 of file minidsp_spectrum.c.

◆ _mel_num_mels

unsigned _mel_num_mels = 0

static

Definition at line 44 of file minidsp_spectrum.c.

◆ _mel_sample_rate

double _mel_sample_rate = 0.0

static

Definition at line 45 of file minidsp_spectrum.c.

◆ _spec_in

double* _spec_in = nullptr

static

Definition at line 22 of file minidsp_spectrum.c.

◆ _spec_N

unsigned _spec_N = 0

static

Definition at line 21 of file minidsp_spectrum.c.

◆ _spec_out

fftw_complex* _spec_out = nullptr

static

Definition at line 23 of file minidsp_spectrum.c.

◆ _spec_plan

fftw_plan _spec_plan = nullptr

static

Definition at line 24 of file minidsp_spectrum.c.

◆ _stft_win

double* _stft_win = nullptr

static

Definition at line 33 of file minidsp_spectrum.c.

◆ _stft_win_N

unsigned _stft_win_N = 0

static

Definition at line 34 of file minidsp_spectrum.c.

Macros

Functions

Variables

Detailed Description

Macro Definition Documentation

◆ MD_F0_FFT_GLOBAL_RATIO

◆ MD_F0_FFT_PROMINENCE_RATIO

◆ MD_MFCC_LOG_FLOOR

Function Documentation

◆ _mel_setup()

◆ _mel_teardown()

◆ _spec_free()

◆ _spec_setup()

◆ _spec_teardown()

◆ _stft_setup()

◆ _stft_teardown()

◆ MD_f0_fft()

◆ md_hz_to_mel()

◆ MD_magnitude_spectrum()

◆ MD_mel_energies()

◆ MD_mel_filterbank()

◆ md_mel_to_hz()

◆ MD_mfcc()

◆ md_parabolic_offset()

◆ MD_phase_spectrum()

◆ MD_power_spectral_density()

◆ MD_shutdown()

◆ MD_stft()

◆ MD_stft_num_frames()

Variable Documentation

◆ _mel_fb

◆ _mel_max_freq_hz

◆ _mel_min_freq_hz

◆ _mel_N

◆ _mel_num_mels

◆ _mel_sample_rate

◆ _spec_in

◆ _spec_N

◆ _spec_out

◆ _spec_plan

◆ _stft_win

◆ _stft_win_N