|
Packit |
c948fe |
/*
|
|
Packit |
c948fe |
* Copyright (C) 2004 Steve Harris
|
|
Packit |
c948fe |
*
|
|
Packit |
c948fe |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
c948fe |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
c948fe |
* the Free Software Foundation; either version 2 of the License, or
|
|
Packit |
c948fe |
* (at your option) any later version.
|
|
Packit |
c948fe |
*
|
|
Packit |
c948fe |
* This program is distributed in the hope that it will be useful,
|
|
Packit |
c948fe |
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
c948fe |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
Packit |
c948fe |
* GNU General Public License for more details.
|
|
Packit |
c948fe |
*
|
|
Packit |
c948fe |
* $Id: distortion-test.c $
|
|
Packit |
c948fe |
*/
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
#include <unistd.h>
|
|
Packit |
c948fe |
#include <string.h>
|
|
Packit |
c948fe |
#include <math.h>
|
|
Packit |
c948fe |
#include <stdio.h>
|
|
Packit |
c948fe |
#include <stdint.h>
|
|
Packit |
c948fe |
#include <stdlib.h>
|
|
Packit |
c948fe |
#include <complex.h>
|
|
Packit |
c948fe |
#include <fftw3.h>
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
#include "gdither.h"
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
#define SIZE 16384
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
#ifndef M_PI
|
|
Packit |
c948fe |
#define M_PI 3.1415926
|
|
Packit |
c948fe |
#endif
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
float in[SIZE];
|
|
Packit |
c948fe |
int32_t out[SIZE];
|
|
Packit |
c948fe |
double sin_tbl[SIZE];
|
|
Packit |
c948fe |
double sinless_tbl[SIZE];
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
void find_snr(GDitherType dt, const char *desc);
|
|
Packit |
c948fe |
void find_res(GDitherType dt, const char *desc);
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
int main(int argc, char *argv[])
|
|
Packit |
c948fe |
{
|
|
Packit |
c948fe |
int i;
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
for (i=0; i
|
|
Packit |
c948fe |
sin_tbl[i] = sin((double)i * M_PI * 0.05) * 0.01;
|
|
Packit |
c948fe |
in[i] = (float)sin_tbl[i];
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
find_snr(GDitherNone, "no dithering");
|
|
Packit |
c948fe |
find_snr(GDitherRect, "rectangular dithering");
|
|
Packit |
c948fe |
find_snr(GDitherTri, "triangular dithering");
|
|
Packit |
c948fe |
find_snr(GDitherShaped, "noise-shaped dithering");
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
printf("ok\n\n");
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
find_res(GDitherNone, "no dithering");
|
|
Packit |
c948fe |
find_res(GDitherRect, "rectangular dithering");
|
|
Packit |
c948fe |
find_res(GDitherTri, "triangular dithering");
|
|
Packit |
c948fe |
find_res(GDitherShaped, "noise-shaped dithering");
|
|
Packit |
c948fe |
printf("ok\n");
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
return 0;
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
void find_res(GDitherType dt, const char *desc)
|
|
Packit |
c948fe |
{
|
|
Packit |
c948fe |
int i, j;
|
|
Packit |
c948fe |
GDither *ds = gdither_new(dt, 1, 16, 0);
|
|
Packit |
c948fe |
int16_t out16[SIZE];
|
|
Packit |
c948fe |
fftw_plan plan_rc;
|
|
Packit |
c948fe |
fftw_complex *freq = fftw_malloc(sizeof(fftw_complex) * (SIZE + 1));
|
|
Packit |
c948fe |
double amp[SIZE];
|
|
Packit |
c948fe |
double floor_peak = -1000.0;
|
|
Packit |
c948fe |
double floor_sum = 0.0;
|
|
Packit |
c948fe |
double floor_cnt = 0.0;
|
|
Packit |
c948fe |
double harm_peak = -1000.0;
|
|
Packit |
c948fe |
double harm_sum = 0.0;
|
|
Packit |
c948fe |
double harm_cnt = 0.0;
|
|
Packit |
c948fe |
int harmonic[SIZE];
|
|
Packit |
c948fe |
//char tmp[256];
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
plan_rc = fftw_plan_dft_r2c_1d(SIZE, sinless_tbl, freq, FFTW_PRESERVE_INPUT);
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
gdither_runf(ds, 0, SIZE, in, out16);
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
for (i=0; i
|
|
Packit |
c948fe |
sinless_tbl[i] = ((double)out16[i] / 32768.0) - sin_tbl[i];
|
|
Packit |
c948fe |
sinless_tbl[i] *= -0.5 * cos(2.0f * M_PI * (float) i /
|
|
Packit |
c948fe |
(float) SIZE) + 0.5;
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
fftw_execute(plan_rc);
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
/* divide bins into ones that are where harmonic peaks could be
|
|
Packit |
c948fe |
and ones that aren't */
|
|
Packit |
c948fe |
for (i=0; i
|
|
Packit |
c948fe |
harmonic[i] = 0;
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
for (i=1; i<20; i++) {
|
|
Packit |
c948fe |
for (j=-4; j<5; j++) {
|
|
Packit |
c948fe |
harmonic[i*SIZE/40 + j] = 1;
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
for (i=2; i
|
|
Packit |
c948fe |
amp[i] = 20.0 * log10(cabs(freq[i]));
|
|
Packit |
c948fe |
if (harmonic[i]) {
|
|
Packit |
c948fe |
if (amp[i] > harm_peak) harm_peak = amp[i];
|
|
Packit |
c948fe |
if (1 || amp[i] > -100.0) {
|
|
Packit |
c948fe |
harm_sum += amp[i];
|
|
Packit |
c948fe |
harm_cnt++;
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
} else {
|
|
Packit |
c948fe |
if (amp[i] > floor_peak) floor_peak = amp[i];
|
|
Packit |
c948fe |
floor_sum += amp[i];
|
|
Packit |
c948fe |
floor_cnt++;
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
//sprintf(tmp, "%f\t%g\n", 48000.0 / (double)SIZE * (double)i, amp[i]);
|
|
Packit |
c948fe |
//write(5, tmp, strlen(tmp));
|
|
Packit |
c948fe |
//printf("%d\t%g\n", i, amp_db);
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
printf("%s mean harmonic resonace = %.1fdB\n", desc,
|
|
Packit |
c948fe |
harm_sum / harm_cnt - floor_sum / floor_cnt);
|
|
Packit |
c948fe |
printf("%s peak harmonic resonance = %.1fdB\n", desc,
|
|
Packit |
c948fe |
harm_peak - floor_peak);
|
|
Packit |
c948fe |
if (dt != GDitherNone && harm_sum / harm_cnt - floor_sum / floor_cnt > 6.0) {
|
|
Packit |
c948fe |
printf("excessive resonance peaks, failing\n");
|
|
Packit |
c948fe |
exit(1);
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
if (dt != GDitherNone && harm_peak - floor_peak > 6.0) {
|
|
Packit |
c948fe |
printf("excessive resonance peaks, failing\n");
|
|
Packit |
c948fe |
exit(1);
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
void find_snr(GDitherType dt, const char *desc)
|
|
Packit |
c948fe |
{
|
|
Packit |
c948fe |
GDither *ds = gdither_new(dt, 1, GDither32bit, 24);
|
|
Packit |
c948fe |
int i;
|
|
Packit |
c948fe |
double sum_sq = 0.0;
|
|
Packit |
c948fe |
double val;
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
gdither_runf(ds, 0, SIZE, in, out);
|
|
Packit |
c948fe |
for (i=0; i
|
|
Packit |
c948fe |
sinless_tbl[i] = ((double)out[i] / 2147483648.0) - sin_tbl[i];
|
|
Packit |
c948fe |
sum_sq += sinless_tbl[i] * sinless_tbl[i];
|
|
Packit |
c948fe |
//printf("%g\n", sinless_tbl[i]);
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
val = 20.0 * log10(sqrt(sum_sq / (double)SIZE));
|
|
Packit |
c948fe |
printf("RMS SNR for %s = %.0fdB\n", desc, val);
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
switch (dt) {
|
|
Packit |
c948fe |
case GDitherNone:
|
|
Packit |
c948fe |
break;
|
|
Packit |
c948fe |
case GDitherRect:
|
|
Packit |
c948fe |
if (val > -141.0) {
|
|
Packit |
c948fe |
printf("excessive noise, failed\n");
|
|
Packit |
c948fe |
exit(1);
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
break;
|
|
Packit |
c948fe |
case GDitherTri:
|
|
Packit |
c948fe |
if (val > -144.0) {
|
|
Packit |
c948fe |
printf("excessive noise, failed\n");
|
|
Packit |
c948fe |
exit(1);
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
break;
|
|
Packit |
c948fe |
case GDitherShaped:
|
|
Packit |
c948fe |
if (val > -130.0) {
|
|
Packit |
c948fe |
printf("excessive noise, failed\n");
|
|
Packit |
c948fe |
exit(1);
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
break;
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
}
|
|
Packit |
c948fe |
|
|
Packit |
c948fe |
/* vi:set ts=8 sts=4 sw=4: */
|