<!-- Archive by https://little-link.co/archive -->
% Except for the creation of the sampled signal, this file was made by someone
% else, and is not in any way the creation of Kurt Sanner
% recent modification made by Kurt Sanner
% December 16, 1999
% Elec 301 Group Project
% this file performs the sigma - delta processing of the sampled signal after
%it has been 16-times oversampled.


% Demonstrate simulateDSM and simulateSNR


echo off;
clear;


clc
fprintf(1,'\t\t\tDiscrete-Time Simulation\n');
if exist('LiveDemo','var') == 0
    LiveDemo=0;
end
R=32;


  % first, we form a matrix for all the frequencies we want from 0 to 32000 Hz
  % We assume they are normalized to 1 Volt maximum for loudest signal
  % recordable.

  w=0:20*2*pi:32000*2*pi;
  numfreqs = length (w);

  % Next, we build the inverse Chebyshev filter
  Gp = 2;
  Gs = 30;
  Wp = 2*pi*18000;
  Ws = 2*pi*22000;
  [n, Wc] = cheb2ord (Wp, Ws, Gp, Gs, 's');
  [num, den] = cheby2 (n, Gs, Wc, 's');

  % finally, we create the magnitude and phase shift vectors for each value of w
  [MAG, PHASE] = bode(num, den, w);

  %x3 is the filtered input

  % there should be 1601 individual frequencies w for which the magnitude and
  % phase shift have been calculated.  We need the particular magnitude and
  % phase shift for the frequency of the sinusoid that is part of the input.

  % sin3 is a 30kHz filtered sine wave.  It yields 30 thousand cylces every
  % second. It repeats a cycle every 4.41/3 sample points, when sampling at 44.1
  % kHz.  Notice that this will result in aliasing, since we are not obeying the
  % Nyquist sampling rate requirement.  Therfore, the filter was designed to
  % eliminate this signal as much as possible.

  % shift3 is for the 30kHz signal
  % shift2 is for the 10kHz signal
  % shift1 is for the 1 kHz signal

  shift3 = 2*pi/360 * PHASE (round(numfreqs * 30000/32000));
  shift2 = 2*pi/360 * PHASE (round(numfreqs * 10000/32000));
  shift1 = 2*pi/360 * PHASE (round(numfreqs * 1000/32000));

shiftsin3 = .3*sin(2*pi*[0:44100]/(4.41/3) + shift3);
shiftsin2 = .1*sin(2*pi*[0:44100]/4.41 + shift2);
shiftsin1 = .1*sin(2*pi*[0:44100]/44.1 + shift1);

  filt3 = MAG (round(numfreqs * 30000/32000)) * shiftsin3;
  filt2 = MAG (round(numfreqs * 10000/32000)) * shiftsin2;
  filt1 = MAG (round(numfreqs * 1000/32000)) * shiftsin1;

% We assume the signal has already been normalized, so that 1 Volt is the peak
% possible for the amplitude.

qfilt3 = round(filt3*65536)/65536;
qfilt2 = round(filt2*65536)/65536;
qfilt1 = round(filt1*65536)/65536;

  x3 = qfilt1 + qfilt2 + qfilt3;


x4 = x3(1:441);

for k = 1:length(x4)-1,
	val1 = x4(k);
	val2 = x4(k+1);
	step = (val2 - val1)/16;
	bigmatrix (16*(k-1)+1) = val1;
	for n = 1:15,
		inval = (round((val1+step*n)*65536)) / 65536;
		bigmatrix(16*(k-1)+n+1) = inval;
	end
end
bigmatrix(16*k + 1) = val2;


stem (bigmatrix, 'r');
ylabel ('amplitude, normalized to 1 volt')
title ('stem plot of oversampled signal after filtering and sampling');
xlabel ('sample points, time equal to (sample point)#/705,600');


pause;


echo on;

H = synthesizeNTF(5,R,1);
N = 7041;
fB = ceil(N/(2*R)); f=floor(2/3*fB);
u = bigmatrix;
v = simulateDSM(u,H); 
echo off;

figure(1); clf;
t = 0:200;
stairs(t, u(t+1),'r');
hold on;
stairs(t,v(t+1),'g');
axis([0 200 -1.2 1.2]);
xlabel('Sample Number');
ylabel('u, v');
title('Modulator Input & Output');
if LiveDemo
    set(1,'position',[4 205 777 355]);
    changeFig(18,2,12);
else
	fprintf(1,'paused\n');
end
pause
if LiveDemo
    set(1,'position',[4 360 300 200]);
    changeFig;
end

echo on
spec = fft(v.*hann(N))/(N/4);
echo off;
figure(2); clf;
plot(linspace(0,0.5,N/2), dbv(spec(1:N/2)))
axis([0 0.5 -120 0]);
text(0.25, -90, sprintf('BW=%7.5f',1.5/N));
grid on;
xlabel('Normalized Frequency')
ylabel('dBFS')
title('Output Spectrum');
if LiveDemo
    set(2,'position',[80 57 662 505]);
    changeFig(18,2,12);
else
	fprintf(1,'paused\n');
end
pause

echo on
snr = calculateSNR(spec(3:fB),f-2);
echo off;
text_handle = text(0.25,-110, sprintf('SNR = %4.1fdB @ OSR=%d\n',snr,R));
if LiveDemo
    set(text_handle,'fontsize',18);
    pause;
	set(2,'position',[328 360 450 200]);
    changeFig;
	pause;	% This ensures that the graphs are updated during the next calc.
else
	fprintf(1,'paused\n'); pause
end

echo on
[snr_pred,amp_pred] = predictSNR(H,R);
[snr,amp] = simulateSNR(H,R);
echo off
figure(3); clf;
plot(amp_pred,snr_pred,'-',amp,snr,'og');
grid on;
figureMagic([-100 0], 10, 1, [0 100], 10, 1);
xlabel('Input Level (dBFS)');
ylabel('SNR (dB)');
title('SNR curve- theory and simulation');
text(-49,15,sprintf('peak SNR = %4.1fdB\n(OSR=%d)',max(snr),R));
if LiveDemo
    set(3,'position',[80 57 662 505]);
    changeFig(18,2,12);
else
	fprintf(1,'paused\n');
end
pause
if LiveDemo
    set(3,'position',[450 60 330 260]);
    changeFig;
end

clc; 
fprintf(1,'\t\t\t Bandpass Modulator\n');
echo on
f0 = 1/8; R=64;
H = synthesizeNTF(8,R,1,[],f0);
fB = ceil(N/(2*R)); f=round(f0*N + 1/3*fB);
u = 0.5*sin(2*pi*f/N*[0:N-1]);	% half-scale sine-wave input
v = simulateDSM(u,H); 
echo off;

figure(1); clf;
t = 0:100;
stairs(t, u(t+1),'r');
hold on;
stairs(t,v(t+1),'g');
axis([0 100 -1.2 1.2]);
xlabel('Sample Number');
ylabel('u, v');
title('Modulator Input & Output');
if LiveDemo
    set(1,'position',[4 205 777 355]);
    changeFig(18,2,12);
else
	fprintf(1,'paused\n');
end
pause
if LiveDemo
    set(1,'position',[4 360 300 200]);
    changeFig;
end

echo on
spec = fft(v.*hann(N))/(N/4);
echo off;
figure(2); clf;
plot(linspace(0,0.5,N/2), dbv(spec(1:N/2)))
axis([0 0.5 -140 0]);
text(0.25, -90, sprintf('BW=%7.5f f_s',1.5/N));
f1 = round((f0-0.25/R)*N);
f2 = round((f0+0.25/R)*N);
snr = calculateSNR(spec(f1:f2),f-f1+1);
text(0.25,-110, sprintf('SNR = %4.1f dB\n(OSR=%d)',snr,R));
grid on;
xlabel('Normalized Frequency')
ylabel('dBFS')
title('Output Spectrum');
if LiveDemo
    set(2,'position',[80 57 662 505]);
    changeFig(18,2,12);
    pause;
	set(2,'position',[328 360 450 200]);
    changeFig;
	pause;	% This ensures that the graphs are updated during the next calc.
else
	fprintf(1,'paused\n'); pause
end

echo on
[snr_pred,amp_pred] = predictSNR(H,R,[],f0);
[snr,amp] = simulateSNR(H,R,[],f0);
echo off
figure(3); clf;
plot(amp_pred,snr_pred,'-b',amp,snr,'og');
grid on;
figureMagic([-110 0], 10, 1, [0 110], 10, 1);
xlabel('Input Level (dBFS)');
ylabel('SNR (dB)');
title('SNR curve- theory and simulation');
text(-69,15,sprintf('peak SNR = %4.1fdB\n(OSR=%d)',max(snr),R));
if LiveDemo
    set(3,'position',[80 57 662 505]);
    changeFig(18,2,12);
else
	fprintf(1,'paused\n');
end
pause
if LiveDemo
    set(3,'position',[450 60 330 260]);
    changeFig;
end

clc;
fprintf(1,'\t\t\t 16-level Lowpass\n');
fprintf(1,'\t\t\t(no SNR prediction)\n\n');
echo on
R=16; M = 16;
H = synthesizeNTF(5,R,1,5);
N = 8192;
fB = ceil(N/(2*R)); f=floor(2/3*fB);
u = 0.5*M*sin(2*pi*f/N*[0:N-1]);	% half-scale sine-wave input
v = simulateDSM(u,H,M); 
echo off;

figure(1); clf;
t = 0:100;
stairs(t, u(t+1),'r');
hold on;
stairs(t,v(t+1),'g');
figureMagic([0 100],10,2, [-M+1 M-1],2,3);
xlabel('Sample Number');
ylabel('u, v');
title('Modulator Input & Output');
if LiveDemo
    set(1,'position',[4 205 777 355]);
    changeFig(18,2,12);
else
	fprintf(1,'paused\n');
end
pause
if LiveDemo
    set(1,'position',[4 360 300 200]);
    changeFig;
end

echo on
spec = fft(v.*hann(N))/(N/4);
echo off;
figure(2); clf;
plot(linspace(0,0.5,N/2), dbv(spec(1:N/2)))
snr = calculateSNR(spec(3:fB),f-2);
text(0.25,-110, sprintf('SNR = %4.1fdB @ OSR=%d\n',snr,R));
text(0.25, -90, sprintf('BW=%7.5f',1.5/N));
figureMagic([0 0.5], 0.1, 2, [-160 20], 20, 2);
grid on;
xlabel('Normalized Frequency')
ylabel('dBFS')
title('Output Spectrum');
if LiveDemo
    set(2,'position',[80 57 662 505]);
    changeFig(18,2,12);
    pause;
	set(2,'position',[328 360 450 200]);
    changeFig;
	pause;	% This ensures that the graphs are updated during the next calc.
else
	fprintf(1,'paused\n'); pause
end

echo on
[snr,amp] = simulateSNR(H,R,[],[],M);
echo off;
figure(3); clf;
plot(amp,snr,'og',amp,snr,'--b');
grid on;
figureMagic([-120 0], 10, 2, [0 120], 10, 2);
xlabel('Input Level (dBFS)');
ylabel('SNR (dB)');
title('SNR curve');
text(-79,15,sprintf('peak SNR = %4.1fdB\n(OSR=%d)\n',max(snr),R));
if LiveDemo
    set(3,'position',[80 57 662 505]);
    changeFig(18,2,12);
else
	fprintf(1,'paused\n');
end
pause
if LiveDemo
    set(3,'position',[450 60 330 260]);
    changeFig;
end

echo off