Test waveforms for Scott Wurcer's digital RIAA article in Volume 10

The following collection of test- and source files was submitted by Scott Wurcer to complement his article on Volume 10 on digital RIAA correction.

For those of you who wish to download the material, there is a pair of .zip files at the end which contain all the files and matrial posted here.

FIR Filters for RIAA    

There are several sets of files containing minimum phase FIR filter coefficients for RIAA equalization in both text and wav format. There are files for 48K, 96K and 192K sampling rates as well as duplicates employing the extra 50K time constant that some recommend to compensate for cutter head roll-off/phase. The text files are in the SoX .dat format as shown below. 

; Sample Rate 48000
; Channels 1
0.000000000e+00 5.966590663e-08
2.083333333e-05 -5.982466236e-08
4.166666667e-05 5.998405294e-08
6.250000000e-05 -6.014408199e-08
…

Here the first two lines contain the sample rate and the number of channels and each subsequent line contains the sample time and value at that time. Excel can be set to import just the second column skipping the first two lines from its import dialog so a simple text file of just the coefficients can be easily created. The files all have the peak value and hence the processing delay at different points, in general a slight delay which is equal in both channels does not matter for the purpose of simply listening. The text version of the files are all normalized to 1.0 so the peak can be found by a simple search in the second column if this is important in a particular application.

SoX Batch Files

It’s very easy to create a drag and drop .bat file to facilitate file conversions and other signal processing in SoX. It only requires an awareness of a small subset of DOS batch file commands and syntax. For instance, to convert a .wav file to a .dat (plain text) file one would do this…

cd %~dp0
FOR %%A IN (%*) DO sox -D %%A %%~nA.dat
pause

When you drag a set of files onto this batch file ~dp0 contains the directory that the batch file is in and you cd to that directory first. Then for each file in all the files dropped run sox with no dither (-D) and save the result in this directory as a .datfile with the same name. Whenever possible SoX uses the file extension to figure out the “right thing” to do so dropping a .wav file results in the simple translation of .wav to .dat with no dither. There are many web resources for translating the sometimes cryptic batch file syntax, and the SoX manual describes building a command line which can include both IIR and FIR filters as well as several options for adjusting gain or normalization of the output. The reverse process…

cd %~dp0
FOR %%A IN (%*) DO sox -D %%A -e signed -b32 %%~nA.wav gain -2
pause

Here we read a .dat file and set the output format to 32bit signed integers, the highest resolution 32 bit format. I have set the gain to 2dB below full scale for .wav files because I have run across some low cost sound cards that clip internally with files that go to full scale. For both of these SoX has to be in your $path or the same folder as the batch file.

min_phase_FIR_RIAA_48k.dat    min_phase_FIR_RIAA_96k.dat   min_phase_FIR_RIAA_192k.dat

min_phase_FIR_RIAA_48k.wav   min_phase_FIR_RIAA_96k.wav   min_phase_FIR_RIAA_192k.wav

min_phase_FIR_RIAAp50_192k.dat   min_phase_FIR_RIAAp50_96k.dat   min_phase_FIR_RIAAp50_48k.dat

min_phase_FIR_RIAAp50_48k.wav   min_phase_FIR_RIAAp50_96k.wav   min_phase_FIR_RIAAp50_192k.wav

The following two batch files facilitate conversion of one file type to the other:  dattowav.bat      wavtodat.bat

Test Waveforms 

These folders contain test waveforms (all in .wav format) with several possible uses. All files come in both a 48kHz and 96kHz versions. They are based on the concept of spectral contamination in the presence of a noise or music like stimulus. The first set consists of four multi-tone files where there are 30 tones at 1/3 octave spacing with random phases. In the time domain it sounds a little like one big 6 handed chord. There are both a general use, flat, version and a version with RIAA pre-emphasis applied which can be used as a simple test for RIAA conformance. The numerical noise filling the space between the tones is below the noise floor of a 24 bit DAC so any imperfections in the signal path will fill the space between the tones with IM/THD components. The files with a flat spectrum are generally useful for DAC/ADC testing as well as any amplifier/speaker chain. It may be difficult to find a DAC or ADC to reproduce these without excessive spectral contamination of their own, in that case a comparative analysis might be useful. These file are designed to be analyzed with a 65536 point FFT and with no (rectangular) windowing though some window functions have low enough leakage to still be useful. With any other length FFT the relationship of the tones to the FFT bins is lost. The next two folders contain pseudo-random noise waveforms containing equal energy at every FFT bin but again with random phase so the resulting time domain waveform looks and sounds like random noise. Each waveform comes in versions with all tones and with four logarithmically spaced gaps to again do the spectral contamination test, and with or without RIAA pre-emphasis for a total of eight files. Since the RIAA pre-emphasized files see the full 40dB of the RIAA curve and do not take into account that the spectral content of most music rolls off with frequency, they are an extreme test of your signal path.

30 tone multitone test files:                30multitone_IRIAA_96.wav   30multitone_flat_96.wav​             30multitone_IRIAA_48.wav   30multitone_flat_48.wav

pseudo-noise flat test files:                  psrnoise_flat_48.wav            psrnoise_flat_48_holes.wav        psrnoise_flat_96.wav            psrnoise_flat_96_holes.wav

pseudo-noise_inverse RIAA test files:  psrnoise_IRIAA_48.wav        psrnoise_IRIAA_48_holes.wav    psrnoise_IRIAA_96.wav       psrnoise_IRIAA_96_holes.wav

Coefficient tables

To close this section, here are the tables with filter coefficients as shown in the Article Appendix. If you want to play with them you can cut-'n-paste from these files:

v10 sw appendix a tables      v10 sw appendix b tables

Python functions and filter Optimizer

Optimizer

This is a very basic example of using scipy’s brute force optimizer to optimize the z domain poles and zeros to match an s domain prototype filter. Since the optimizer starts by doing an exhaustive search over a range of values for any number of variables, this becomes very slow rapidly as the number of variables or fineness of the search grid grows. This example is set up for a peaking equalizer where the gain , Q, and frequency of interest (fo) are specified. I suggest the reader familiarize themselves with the optimizer by reading the online description as well as The Audio EQ Cookbook. 
The function fun(x) is minimized by the optimizer via the array (x) of variables. Here x[0] is fo, x[1] is A (gain), and x[2] is Q. The poles and zeros are extracted as the roots of the characteristic polynomial and the peak to peak error between the z domain response and the s domain response is returned.
The values in the body of the routine will need to be edited for different frequencies, other parameters, and different types of filters. Start and stop give a range of frequencies over which to optimize. The ranges are the set of ranges for the array of variables (x) passed by the optimizer (brute()). I use these to multiply the filter parameters so (.8, 1.2) is an ~ +- 20% range while the variable Ns is the fineness of the search. There is lots of room for experimentation with the values that are left up to the user.
The optimizer returns the optimized values for the x array which are then used to compute the optimized poles and zeros and the coefficients for the IIR filter which are printed out at the end.
This example is for a 5dB, 10kHz peaking filter with a Q of .5 and a sampling frequency of 48kHz. I chose to optimize only to 12kHz because increasing this rapidly degrades the fit due to the frequency warping. In fact it is easy to get no better a fit than by simply using the suggested “tweek” to the Q mentioned in the main article. Running the script generates output which consists of the poles and zeros of the prototype followed by the optimum z domain poles and zeros with the computed IIR filter coefficients.

optimizer.py note: rename to optimizer.py to run!    

Python functions library

These are the basic functions that I use to go back and forth between the s and z domain, compute IIR coefficients and plot phase and frequency responses. As mentioned in the article I am not a programmer by trade and there are probably many improvements possible here. It is easy with Python to get a little confused around arrays, lists, and other ways of packing and unpacking data. Usually what seems obvious works but there are times when an explicit declaration is necessary. Also there is not thorough checking for bad inputs so they will simply fail or give a bad output. As noted also I use positive frequency for both the s and z domain in my computations, the built-in roots() function returns negative values as it should and they are inverted as necessary in my code.

dB(x) - Returns 20*log10(abs(x)) for any length argument.

f_warp(f, fs) – Returns the “warped” frequency of interest (f) for IIR filter coefficients. Fs is the sampling frequency.

f_from_z(z, fs) – Returns a frequency from the point z inside the unit circle.

z_from_f(f, fs) – Returns the point z inside the unit circle corresponding to frequency f.

Fs_at_f(Poles, Zeros, f, norm = 0) – Returns the complex value/s for s domain frequency response at list of frequencies f. The transfer function is described by its poles and zeros. If norm is true (1) the output is normalized to 1 at the maximum value.

Fz_at_f(Poles, Zeros , f, fs, norm = 0) – Returns the same for the z domain poles and zeros (here fs is needed also). The pole and zero frequencies are those that would be returned from f_from_z(z, fs).

z_coeff(Poles, Zeros, fs, g, fg, fo = 'none') – Returns the a and b coefficients of the IIR filter defined by the poles and zeros. G and fg are the gain in db and frequency at which that gain is set. Fo is the frequency of interest, that is, for a peaking filter fo would be the peak while for a low pass or high pass filter it would be the break point. If ‘none’ or not specified there is no frequency warping and all mapping is simply by the bilinear transform. 

biquad_filter(xin, z_coeff) – Returns xin filtered by the biquad filter defined by the z domain coefficients, [[a0, a1, a2] ,[ b0,b1,b2]]. By convention the b coefficients are the numerator and the a are the denominator with a0 = 1.

write_txt(filename, y) – Writes a simple text file of the values of the array y.

read_txt(filename) – Returns an array from the values in a text file.

read_dat(filename) – Returns a mono or stereo file in the SoX .dat format.

write_dat(filename, sr, left, right='none') – Writes a file in the SoX .dat format. Sr = the sample rate. Left and right are equal sized arrays of time samples.

plot_fft_log(fs, data, diff = 'none', start = 0,stop = 0, delay = 0) – Plots the FFT of a time domain waveform. Fs is the sampling frequency and data are the samples. If diff is specified the difference between data and diff is plotted. Start and stop are the start and stop frequencies for the plot, 20Hz and 20kHz are the defaults. A delay in samples is specified to compensate for phase wrapping, the default is 0. This function compensates for actual bin frequencies so, for example, a 65536 point FFT of a waveform sampled at 48kHz gives the right frequencies.

plot_fft_lin(fs, data,diff = 'none', start = 0, stop = 0, delay = 0) – Same function for a linear with frequency plot.

plot_dat(sound) – Plots a sound file read by read_dat() or in the .dat format.

plot_gain_phase(Poles, Zeros, zPoles = [], zZeros = [], fs = 48000.,fo = 1000. ,start = 20, stop = 20000,style = 0) – Plots the gain and phase vs frequency of a frequency response defined by a set of Poles and Zeros. If zPoles and zZeros are specified it becomes a comparison between the s and z domain with fs being the sampling frequency. Fo is the frequency where both responses are equalized. Start and stop are again the limits on the plot while style = 1 is a difference plot and style = 0 is a plot of both responses on the same axes.

min_phase_FIR(Poles, Zeros, fs, outfile, plot = 0) – Generates a minimum phase FIR filter in mono .dat format from the s domain poles and zeros at a sampling frequency fs. The output is written to ‘outfile’ and it expects a HOME directory to exist. You can change the code to output to whatever you want. Plot = 1 plots the result so you can examine it.

Linear_Audio_Lib.py note: rename to Linear_Audio_Lib.py to run! 

All test- and data waveform files                     All Python optimizer example&source files