diff options
| -rw-r--r-- | src/common.h | 2 | ||||
| -rw-r--r-- | src/dft.c | 17 | ||||
| -rw-r--r-- | src/dft.h | 10 | ||||
| -rw-r--r-- | src/string.h | 2 | ||||
| -rw-r--r-- | src/wav.c | 6 | ||||
| -rw-r--r-- | src/wav.h | 14 | ||||
| -rw-r--r-- | tests/dft.c | 53 | ||||
| -rw-r--r-- | tests/wav.c | 123 |
8 files changed, 210 insertions, 17 deletions
diff --git a/src/common.h b/src/common.h index 1154bdf..f1e74f8 100644 --- a/src/common.h +++ b/src/common.h @@ -13,7 +13,7 @@ #endif // !PI #ifndef F64_EPSILON -#define F64_EPSILON = 1e-9 +#define F64_EPSILON 1e-9 #endif // !F64_EPSILON #define F64_EQ(x, y, eps) (fabs((x) - (y)) <= (eps)) @@ -139,6 +139,7 @@ wav_data_t dft_data_to_wav(smrt_arena_t *arena, dft_data_t dft, u64 sample_rate, f64 p = dft.phases[freq_index]; amp += cos((t * f * 2.0 * PI) + p) * a; } + amp = amp < -1.0 ? -1.0 : amp > 1.0 ? 1.0 : amp; data[s_num] = (u8)((amp+1.0) * (UINT8_MAX/2)); } @@ -148,9 +149,12 @@ wav_data_t dft_data_to_wav(smrt_arena_t *arena, dft_data_t dft, u64 sample_rate, stft_data_t short_time_fourier_transform(smrt_arena_t *arena, u64 window_size, u64 hop_size, f64 *samples, u64 sample_count, u64 sample_rate) { assert(F64_EQ(round(log2(window_size)), log2(window_size), 1e-9) && "STFT input window_size must be a power of 2"); + assert(hop_size != 0 && "STFT input hop_size must be nonzero"); sample_count = ALIGN_UP_POW2(sample_count * sizeof(f64), STFT_SAMPLE_ALIGN_BYTES(window_size)) / sizeof(f64); + assert(sample_count >= window_size && "STFT input sample_count must be at least window_size"); + u64 segment_count = ((sample_count - window_size) / hop_size) + 1; stft_segment_t *segments = SMRTA_ALLOC_ARRAY(arena, stft_segment_t, segment_count); @@ -179,8 +183,11 @@ stft_data_t short_time_fourier_transform(smrt_arena_t *arena, u64 window_size, u return (stft_data_t){ .sample_rate=sample_rate, + .window_size=window_size, + .hop_size=hop_size, .segment_count=segment_count, .segments=segments, + .total_samples=sample_count, }; } @@ -233,10 +240,13 @@ void inverse_fast_fourier_transform(smrt_arena_t *arena, f64 const *real, f64 co smrta_scratch_end(scratch); } -f64 *inverse_short_time_fourier_transform(smrt_arena_t *arena, stft_data_t stft, u64 window_size, u64 hop_size, smrt_arena_t **conflicts, u64 num_conflicts) { +f64 *inverse_short_time_fourier_transform(smrt_arena_t *arena, stft_data_t stft, smrt_arena_t **conflicts, u64 num_conflicts) { + u64 window_size = stft.window_size; + u64 hop_size = stft.hop_size; + smrta_temp_t scratch = smrta_scratch_start(conflicts, num_conflicts); - u64 output_len = (stft.segment_count - 1) * hop_size + window_size; + u64 output_len = stft.total_samples; dft_data_t *frames = SMRTA_ALLOC_ARRAY(scratch.arena, dft_data_t, stft.segment_count); @@ -277,6 +287,9 @@ f64 *inverse_short_time_fourier_transform(smrt_arena_t *arena, stft_data_t stft, } void reconstruct_spectrum(smrt_arena_t *arena, dft_data_t const *data, u64 sample_count, f64 **real_o, f64 **imag_o) { + assert(data->freq_count == (sample_count / 2) + 1 && + "reconstruct_spectrum: data->freq_count must equal (sample_count / 2) + 1"); + f64 *real = *real_o = SMRTA_ALLOC_ARRAY(arena, f64, sample_count); f64 *imag = *imag_o = SMRTA_ALLOC_ARRAY(arena, f64, sample_count); @@ -39,7 +39,10 @@ typedef struct { typedef struct { stft_segment_t * segments; u64 segment_count; + u64 window_size; + u64 hop_size; u64 sample_rate; + u64 total_samples; } stft_data_t; /// Byte alignment required for a samples buffer to safely back @@ -59,6 +62,9 @@ stft_data_t short_time_fourier_transform(smrt_arena_t * arena , u64 sample_rate); /// Run inverse FFT on a complex spectrum (real/imag, length sample_count). +/// No 1/N normalization is applied here: real/imag must already be a correctly-scaled +/// full N-point spectrum (e.g. as produced by reconstruct_spectrum), not a raw +/// forward-FFT output, or the result will be off by a factor of N. /// Pass NULL for either out-param to skip allocating/returning that component void inverse_fast_fourier_transform(smrt_arena_t * arena , f64 const * real , @@ -68,11 +74,9 @@ void inverse_fast_fourier_transform(smrt_arena_t * arena , f64 ** imag_o , smrt_arena_t **conflicts, u64 num_conflicts); -/// Reconstruct samples from STFT data via overlap-add inverse FFT synthesis +/// Reconstruct samples from STFT data via overlap-add inverse FFT synthesis. f64 *inverse_short_time_fourier_transform( smrt_arena_t * arena , stft_data_t stft , - u64 window_size , - u64 hop_size , smrt_arena_t **conflicts, u64 num_conflicts); /// Reconstruct a full N-point complex spectrum from one-sided DFT amplitude/phase data diff --git a/src/string.h b/src/string.h index 7e3cf66..ef6d966 100644 --- a/src/string.h +++ b/src/string.h @@ -17,4 +17,4 @@ strng_t * strng_dup(smrt_arena_t *arena, strng_t const *src); b32 strng_set(strng_t *string, char const *c); void strng_clear(strng_t *string); -#define STRNG_FMT(s) (int)s->len, (char *)((u8*)s+STRNG_BASE_POS) +#define STRNG_FMT(s) (i32)s->len, (char *)((u8*)s+STRNG_BASE_POS) @@ -254,9 +254,11 @@ f64 *read_32bps_float_data(smrt_arena_t *arena, wav_data_t data, u16 num_channel u8 two = data.samples[(i*4*num_channels)+1+channel*4]; u8 three = data.samples[(i*4*num_channels)+2+channel*4]; u8 four = data.samples[(i*4*num_channels)+3+channel*4]; - u32 intermediate = one | two << 8 | three << 16 | four << 24; + u32 intermediate = one | two << 8 | three << 16 | (u32)four << 24; - memcpy(&vs[i], &intermediate, sizeof(u32)); + f32 f; + memcpy(&f, &intermediate, sizeof(f32)); + vs[i] = (f64)f; } return vs; @@ -44,11 +44,11 @@ typedef struct { /// Load WAV data from wav_file. Data will be allocated on arena, with allocation /// size aligned up to align_up_memoryn. -wav_data_t load_wav_file(smrt_arena_t * arena , - FILE *wav_file , - wav_master_chunk_t *master_o , - wav_fmt_chunk_t *format_o , - u64 align_up_memoryn); +wav_data_t load_wav_file(smrt_arena_t * arena , + FILE * wav_file , + wav_master_chunk_t * master_o , + wav_fmt_chunk_t * format_o , + u64 align_up_memoryn); /// Write WAV headers and data to wav_file. b32 write_wav_file(FILE * wav_file , @@ -61,14 +61,14 @@ wav_fmt_chunk_t make_wav_fmt_chunk(u32 num_channels , u16 bits_per_sample); /// Load WAV data into [-1.0, 1.0] f64 amplitudes. Populates -o inputs with information. -/// samples will be amplitude data for each channel from the WAV file, of length sample_count. +/// samples will be amplitude data for each channel from the WAV file, of length sample_count /// align_up_memoryn will align allocation up to nearest multiple (For FFT) void wav_load(smrt_arena_t *arena, FILE *wav, f64 *** samples_o , u16 *channel_count_o , u64 * sample_count_o , u32 * sample_rate_o , u64 align_up_memoryn , - smrt_arena_t **conflicts, u64 num_conflicts ); + smrt_arena_t **conflicts, u64 num_conflicts); f64 * read_8bps_data(smrt_arena_t *arena, wav_data_t data, u16 num_channels, u16 channel, u64 align_up_memoryn); f64 * read_16bps_data(smrt_arena_t *arena, wav_data_t data, u16 num_channels, u16 channel, u64 align_up_memoryn); diff --git a/tests/dft.c b/tests/dft.c index 9c406f5..3f3d7f7 100644 --- a/tests/dft.c +++ b/tests/dft.c @@ -1,6 +1,7 @@ #include "../src/dft.h" #include <math.h> +#include <signal.h> #include <stddef.h> #include <stdint.h> @@ -351,6 +352,56 @@ Test(dft, reconstruct_spectrum_mirrors_conjugate) { smrt_arena_destroy(arena); } +Test(dft, data_to_wav_clamps_out_of_range_amplitude) { + smrt_arena_t *arena = smrt_arena_create(KiB(16), KiB(4), false); + + // two in-phase DC-like bins summing to 1.6, well outside [-1, 1] + dft_data_t dft = { .freq_count = 2 }; + dft.frequencies = SMRTA_ALLOC_ARRAY(arena, f64, 2); + dft.amplitudes = SMRTA_ALLOC_ARRAY(arena, f64, 2); + dft.phases = SMRTA_ALLOC_ARRAY(arena, f64, 2); + dft.frequencies[0] = 0.0; dft.amplitudes[0] = 0.8; dft.phases[0] = 0.0; + dft.frequencies[1] = 0.0; dft.amplitudes[1] = 0.8; dft.phases[1] = 0.0; + + wav_data_t wav = dft_data_to_wav(arena, dft, 8, 1.0); + + cr_assert_eq(wav.sample_count, 8); + for (u64 i = 0; i < wav.sample_count; i++) { + // matches the clamped-amp=1.0 case in data_to_wav_varies_with_time: + // (1.0 + 1.0) * (UINT8_MAX/2 truncated to 127) = 254 + cr_expect_eq(wav.samples[i], 254); + } + + smrt_arena_destroy(arena); +} + +Test(dft, stft_rejects_zero_hop_size, .signal = SIGABRT) { + smrt_arena_t *arena = smrt_arena_create(KiB(16), KiB(4), false); + + f64 samples[8] = { 0 }; + short_time_fourier_transform(arena, 8, 0, samples, 8, 8); +} + +Test(dft, stft_rejects_zero_sample_count, .signal = SIGABRT) { + smrt_arena_t *arena = smrt_arena_create(KiB(16), KiB(4), false); + + f64 samples[1] = { 0 }; + short_time_fourier_transform(arena, 8, 4, samples, 0, 8); +} + +Test(dft, reconstruct_spectrum_rejects_mismatched_freq_count, .signal = SIGABRT) { + smrt_arena_t *arena = smrt_arena_create(KiB(16), KiB(4), false); + + // sample_count=8 expects freq_count == 5; give it 3 instead + dft_data_t data = { .freq_count = 3 }; + data.frequencies = SMRTA_ALLOC_ARRAY(arena, f64, 3); + data.amplitudes = SMRTA_ALLOC_ARRAY(arena, f64, 3); + data.phases = SMRTA_ALLOC_ARRAY(arena, f64, 3); + + f64 *real_o, *imag_o; + reconstruct_spectrum(arena, &data, 8, &real_o, &imag_o); +} + Test(dft, istft_round_trip_recovers_samples) { smrt_arena_t *arena = smrt_arena_create(MiB(1), KiB(4), false); @@ -362,7 +413,7 @@ Test(dft, istft_round_trip_recovers_samples) { stft_data_t stft = short_time_fourier_transform(arena, window_size, hop_size, samples, sample_count, sample_rate); cr_assert_eq(stft.segment_count, 7); - f64 *output = inverse_short_time_fourier_transform(arena, stft, window_size, hop_size, NULL, 0); + f64 *output = inverse_short_time_fourier_transform(arena, stft, NULL, 0); for (u64 i = 0; i < sample_count; i++) { cr_expect(F64_EQ(output[i], samples[i], 1e-9)); diff --git a/tests/wav.c b/tests/wav.c index dae2ef1..7635223 100644 --- a/tests/wav.c +++ b/tests/wav.c @@ -340,3 +340,126 @@ Test(wav, write_sine, .init = write_sine_setup, .fini = write_sine_teardown) { smrt_arena_destroy(arena); } + +static char float32_wav_path[] = "/tmp/steez_testwav_f32_XXXXXX"; +static int float32_wav_fd; + +void float32_wav_setup(void) { + float32_wav_fd = mkstemp(float32_wav_path); + cr_assert_geq(float32_wav_fd, 0, "mkstemp failed"); +} + +void float32_wav_teardown(void) { + unlink(float32_wav_path); + close(float32_wav_fd); +} + +Test(wav, read_32bps_float_round_trip, .init = float32_wav_setup, .fini = float32_wav_teardown) { + FILE *wav_file = fdopen(float32_wav_fd, "wb"); + cr_assert_not_null(wav_file); + + smrt_arena_t *arena = smrt_arena_create(KiB(64), KiB(4), false); + + f32 expected[5] = { 0.0f, 0.5f, -0.5f, 1.0f, -1.0f }; + u64 sample_count = 5; + + wav_data_t data; + data.samples = (u8*)expected; + data.sample_count = sample_count; + + wav_fmt_chunk_t fmtchunk = make_wav_fmt_chunk(1, 44100, 32); + fmtchunk.audio_format = 3; // IEEE 754 float + + cr_assert(write_wav_file(wav_file, &fmtchunk, data)); + + fclose(wav_file); + + FILE *readback_file = fopen(float32_wav_path, "rb"); + cr_assert_not_null(readback_file); + + wav_master_chunk_t mchunk; + wav_fmt_chunk_t read_fmtchunk; + + wav_data_t read_data = load_wav_file(arena, readback_file, &mchunk, &read_fmtchunk, sizeof(f64)); + + fclose(readback_file); + + cr_expect_eq(read_fmtchunk.audio_format, 3); + cr_expect_eq(read_fmtchunk.bits_per_sample, 32); + cr_assert_not_null(read_data.samples); + cr_expect_eq(read_data.sample_count, sample_count); + + f64 *decoded = read_32bps_float_data(arena, read_data, 1, 0, sizeof(f64)); + cr_assert_not_null(decoded); + + for (u64 i = 0; i < sample_count; i++) { + cr_expect(F64_EQ(decoded[i], (f64)expected[i], 1e-6)); + } + + smrt_arena_destroy(arena); +} + +static char pcm24_wav_path[] = "/tmp/steez_testwav_24bit_XXXXXX"; +static int pcm24_wav_fd; + +void pcm24_wav_setup(void) { + pcm24_wav_fd = mkstemp(pcm24_wav_path); + cr_assert_geq(pcm24_wav_fd, 0, "mkstemp failed"); +} + +void pcm24_wav_teardown(void) { + unlink(pcm24_wav_path); + close(pcm24_wav_fd); +} + +Test(wav, read_24bps_round_trip, .init = pcm24_wav_setup, .fini = pcm24_wav_teardown) { + FILE *wav_file = fdopen(pcm24_wav_fd, "wb"); + cr_assert_not_null(wav_file); + + smrt_arena_t *arena = smrt_arena_create(KiB(64), KiB(4), false); + + i32 expected[4] = { 0, 4194304 /* +0.5 */, -4194304 /* -0.5 */, 8388607 /* max */ }; + u64 sample_count = 4; + + u8 raw[4*3]; + for (u64 i = 0; i < sample_count; i++) { + i32 v = expected[i]; + raw[i*3+0] = (u8)( v & 0xFF); + raw[i*3+1] = (u8)((v >> 8) & 0xFF); + raw[i*3+2] = (u8)((v >> 16) & 0xFF); + } + + wav_data_t data; + data.samples = raw; + data.sample_count = sample_count; + + wav_fmt_chunk_t fmtchunk = make_wav_fmt_chunk(1, 44100, 24); + + cr_assert(write_wav_file(wav_file, &fmtchunk, data)); + + fclose(wav_file); + + FILE *readback_file = fopen(pcm24_wav_path, "rb"); + cr_assert_not_null(readback_file); + + wav_master_chunk_t mchunk; + wav_fmt_chunk_t read_fmtchunk; + + wav_data_t read_data = load_wav_file(arena, readback_file, &mchunk, &read_fmtchunk, sizeof(f64)); + + fclose(readback_file); + + cr_expect_eq(read_fmtchunk.bits_per_sample, 24); + cr_assert_not_null(read_data.samples); + cr_expect_eq(read_data.sample_count, sample_count); + + f64 *decoded = read_24bps_data(arena, read_data, 1, 0, sizeof(f64)); + cr_assert_not_null(decoded); + + for (u64 i = 0; i < sample_count; i++) { + f64 expected_amplitude = (f64)expected[i] / (f64)(0b1 << 23); + cr_expect(F64_EQ(decoded[i], expected_amplitude, 1e-9)); + } + + smrt_arena_destroy(arena); +} |