f_samplerate1 = 44100; f_samplerate2 = 96000; f_lcm = 14112000; %! holy crap this would be rediculous to calculate! better use 44100 and 96000...

f_samplerate1 = 44100; f_samplerate2 = 88200; f_exact = 176400; % will compare to the exact waveform values at the upsampled rate of 176.4kHz

T = 1; %period in seconds

t1 = [0:f_samplerate1*T-1]/ f_samplerate1; % 1 second audio, hence FFT frequencies will be spaced by 1/period = 1Hz

t2 = [0:f_samplerate2*T-1]/ f_samplerate2; % 1 second audio, hence FFT frequencies will be spaced by 1/period = 1Hz

te = [0:f_exact*T-1]/ f_exact; %

% this is our mathematically exact definition of the signal

flim = 20000;

fs = 1:flim; % the frequencies present

a = -1+2*rand([1, flim]) ; % cosine coefficients for a 20kHz bandwidth limited signal---white noise

b = -1+2*rand([1, flim]); % sine coefficients for a 20kHz bandwidth limited signal---white noise

A_Squared = sum( a.^2 + b.^2 ); % parseval says sum( abs(y_j)^2 ) = sum ( abs(fft(y_j))^2 ), so we can normalize the signal to a nominal +/- 1 amplitude

a = a/sqrt(A_Squared);

b = b/sqrt(A_Squared); % -> these coefficients tend to give rms amplitude of 0.7071 (ie 1/sqrt(2)), but samples can get up to +/- 3 or so (ie. +/-4 std)

% to be safe, lets cut the samples down by a factor of 4 to avoid clipping

a = a/4;

b = b/4;

Fy1 = zeros(size(t1));

Fy2 = zeros(size(t2));

Fe = zeros(size(te));

Fy1(2:flim+1) = a + 1i*b; % fill in DFT coefficients

Fy1 = 0.5* (Fy1 + fliplr(circshift(Fy1, [0,-1])));

Fy2(2:flim+1) = a + 1i*b; % fill in DFT coefficients

Fy2 = 0.5* (Fy2 + fliplr(circshift(Fy2, [0,-1])));

Fe(2:flim+1) = a + 1i*b; % fill in DFT coefficients

Fe= 0.5* (Fe+ fliplr(circshift(Fe, [0,-1])));

y1 = ifft(Fy1*f_samplerate1, 'symmetric');

y2 = ifft(Fy2*f_samplerate2, 'symmetric');

ye = ifft(Fe*f_exact, 'symmetric');

dither_function = @(y) randn(size(y)); % i dont' know how dither is actually implemented, can somebody point me towards a easy noise-shaping algorithm? I could calculate white noise (like here), and in fft-space, apply a normalization curve, invert the fft, and rescale for a nominal amplitude of +/- 1 bit?

% now we want to cut the samples down to 16bits and 24bits

y1_24 = round( 2^23 * y1 + dither_function(y1)) ; % now, integers in the +/- 2^23 range

y2_16 = round( 2^15 * y2 + dither_function(y2)) ; % now, integers in the +/- 2^15 range

% "playback" on an upsampling DAC or with upsampling DSP

y1_24_up176 = 4*ifft(ifftshift( [zeros([1, length(y1)*1.5]), fftshift(fft(y1_24/(2^23))), zeros([1, length(y1)*1.5])] ),'symmetric');

y2_16_up176 = 2*ifft(ifftshift( [zeros([1, length(y2)*0.5]), fftshift(fft(y2_16/(2^15))), zeros([1, length(y2)*0.5])] ),'symmetric');

% compare to the original mathematically defined waveform

err1 = ye-y1_24_up176;

err2 = ye-y2_16_up176;

% compare the square error of the two? How would you quantify the error, or the quality, etc...?

[sum(err1.^2), sum(err2.^2)]

10*log10([sum(err1.^2), sum(err2.^2)])

figure(1), plot(t1, y1, te, ye-y1_24_up176); % shows the 1st sample of the signal and a line which indicates really good reconstruction (ie small error)

figure(2), plot(t2, y2, te, ye-y2_16_up176); % shows the 2st sample of the signal and a line which indicates really good reconstruction (ie small error)