Noise shaper (1st and 2nd order) simulation as writing practice of GNU Octave

noise_shaper.m

% 1st and 2nd order noise shaper time/frequency domain plot by GNU Octave
 
clear -all
 
%
% parameters
%
 
calc_length	= 2**16;	% number of samples
t_domain_start	= 280;
t_domain_stop	= 850;
 
signal_freq	= 10000;	% Hz
sampling_freq	= 44100;	% Hz
oversampling_ratio	= 256;	% Fs
attenuation	= 1;	% dB to full-scale
 
signal_duration 	= calc_length / sampling_freq;
oversampling_freq	= sampling_freq * oversampling_ratio;
amplitude	= 10 ** (-attenuation/20);
 
%
% input signal
%
 
in	= sinetone( signal_freq, oversampling_freq, signal_duration, amplitude );
 
%
% 1st order noise shaper
%
 
noise	= 0;
for i = 1 : calc_length
qin	= in(i) - noise;
 
if ( qin < 0 )
out1(i)	= -1;
else
out1(i)	= 1;
endif
 
noise	= out1(i) - qin;
end
 
%
% 2nd order noise shaper (SAA7320)
%
 
s2	= 0;
s1	= 0;
for i = 1 : calc_length
if ( s2 < 0 )
out2(i)	= -1;
else
out2(i)	= 1;
endif
 
noise	= out2(i) - s2;
 
s2	= s1 - noise * 2;
s1	= in(i) + noise;
end
 
%
% window function (blackman-harris)
%
 
N	= calc_length-1;
a0	= 0.35875;
a1	= 0.48829;
a2	= 0.14128;
a3	= 0.01168;
 
for i = 0 : N
wndw(i+1)	=0.35875-a1*cos(2*pi*i/N)+a2*cos(4*pi*i/N)-a3*cos(6*pi*i/N);
end
 
%
% applying window function to 1 bit output data
%
 
w_mean	= mean( wndw );
wndw_in1	= (out1.') .* (wndw.');
wndw_in2	= (out2.') .* (wndw.');
 
%
% get power spectrum by FFT
%
 
spectrum1	= 20 * log10( abs( fft( wndw_in1,calc_length ) ) ./ ((calc_length / 2) * w_mean) );
spectrum2	= 20 * log10( abs( fft( wndw_in2,calc_length ) ) ./ ((calc_length / 2) * w_mean) );
spectrum1	= spectrum1( 1:calc_length/2 );
spectrum2	= spectrum2( 1:calc_length/2 );
 
%
% plot timedomain (waveform)
%
 
subplot( 2, 1, 1 );
t	= t_domain_start:1:t_domain_stop;
ph_t0	= plot( t, in(t_domain_start:t_domain_stop) );
hold on
ph_t1	= plot( t, out1(t_domain_start:t_domain_stop) );
ph_t2	= plot( t, out2(t_domain_start:t_domain_stop) );
hold off
set( ph_t1, "color", "blue" );
set( ph_t2, "color", "red" );
legend( ph_t1, "1st order" );
legend( ph_t2, "2nd order" );
axis( [ t_domain_start, t_domain_stop, -1.1, 1.1 ] );
xlabel( "time [sample]" );
ylabel( "amplitude [fullscale = -1..+1]" );
 
%
% plot frequency domain (spectrum)
%
 
subplot( 2, 1, 2 );
f	= 0 : oversampling_freq / (calc_length - 1) : oversampling_freq / 2;
ph_f1	= semilogx( (f.'), spectrum1 );
hold on;
ph_f2	= semilogx( (f.'), spectrum2 );
hold off
set( ph_f1, "color", "blue" );
set( ph_f2, "color", "red" );
legend( ph_f1, "1st order" );
legend( ph_f2, "2nd order" );
axis([ 100, oversampling_freq / 2, -150, 0 ] );
xlabel( "frequency [Hz]" );
ylabel( "power [dBFS]" );
grid on;

Info about this code (Japanese): http://goo.gl/Tb7Pud