The following sections contain the three Appendices which are part of the article by Malcolm Hawksford in Linear Audio Volume 7 on Audio Equalisation.
Theses Appendices are listed in text format and can be used as-is in Matlab or it's equivalent applications.

Appendix 1	Matlab routine for series-feedback phono preamplifier

% RIAA_series_feedback.m
% 04.08.13
% series feedback RIAA preamp with servo
% preferred values can be entered in lines 49 to 60
home; clc; clear; colordef white; close
fprintf('SERIES FEEDBACK RIAA PREAMP WITH SERVO\n')
 
% set ww = 1 to write Excel file of frequency responses
ww=0;
 
% RIAA data
T1=3180*10^-6; % RIAA EQ time constant, second
T2=318*10^-6; % RIAA EQ time constant, second
T3=75*10^-6; % RIAA EQ time constant, second
f0=20; % 3-dB break frequency, Hz
 
% frequency vector (1 Hz to 100 kHz, N points/dec)
N=1000;
f=10.^(0:1/N:6);
 
% component values
R1=110000; % RIAA dc resistor, ohm
R0=100; % feedback resistor to ground, ohm
Rs1=10^6; % servo resistor, ohm
Rs3=10000; % series output resistor of servo
Rx=1000; % resistance of RC filter associate with T4
 
% response correction for series feedback topology
Gdc=1+R1/(R0*Rs3/(R0+Rs3)); % dc gain neglecting servo
Tm2=(Gdc*T2+T1*T3/T2-T1-T3)/(Gdc-1);
T4=T1*T3/(Gdc*T2);
 
% series feedback RIAA preamp
C1=(T1*T3/Tm2)/R1;
C2=(T1+T3-Tm2-T1*T3/Tm2)/R1;
R2=R1*Tm2/(T1+T3-Tm2-T1*T3/Tm2);
 
% servo amplifier components
T0=1/(2*pi*f0);
Rs2=Rs1*Tm2/(T1-Tm2);
Cs0=(R1/Rs3)/(2*pi*f0*(Rs1+Rs2));
Cs1=T1/Rs1;
rs=2*pi*f0*T3*(Rs1+Rs2)*Rs3/R1;
 
% low-pass correction filter, T4
Cx=T4/Rx;
 
% enter preferred values
R0=R0; %#ok<ASGSL>
R1=R1; %#ok<ASGSL>
R2=R2; %#ok<ASGSL>
C1=C1; %#ok<ASGSL>
C2=C2; %#ok<ASGSL>
Rs1=Rs1; %#ok<ASGSL>
Rs2=Rs2; %#ok<ASGSL>
Rs3=Rs3; %#ok<ASGSL>
Cs0=Cs0; %#ok<ASGSL>
Cs1=Cs1; %#ok<ASGSL>
Rx=Rx; %#ok<ASGSL>
Cx=Cx; %#ok<ASGSL>
 
% computed (or preferred value) data for series feedback preamplifier
fprintf('\n\nRIAA EQ COMPONENTS FOR SERIES FEEDBACK MC PRE-AMPLIFIER\n')
R0 %#ok<NOPTS>
R1 %#ok<NOPTS>
R2 %#ok<NOPTS>
C1 %#ok<NOPTS>
C2 %#ok<NOPTS>
Rs1 %#ok<NOPTS>
Rs2 %#ok<NOPTS>
Rs3 %#ok<NOPTS>
Cs0 %#ok<NOPTS>
Cs1 %#ok<NOPTS>
Rx %#ok<NOPTS>
Cx %#ok<NOPTS>
 
% plot reference RIAA frequency response
Griaa=abs(Gdc*(1+i*2*pi*f*T2)./((1+i*2*pi*f*T1).*(1+i*2*pi*f*T3)).*(i*2*pi*f*T0./(1+i*2*pi*f*T0)));
semilogx(f,20*log10(Griaa),'k')
title('Reference RIAA frequency response')
xlabel('Frequency, Hz')
ylabel('Amplitude, dB')
grid; pause(3); close
 
% theoretical frequency response for RIAA TF with and without T0 correction
G0=abs((1+(Gdc-1)*(1+i*2*pi*f*Tm2).*(i*2*pi*f*T0./(1+i*2*pi*f*T0))./((1+i*2*pi*f*T1).*(1+i*2*pi*f*T3)))./(1+i*2*pi*f*T4));
G1=abs(1+(Gdc-1)*(1+i*2*pi*f*Tm2).*(i*2*pi*f*T0./(1+i*2*pi*f*T0))./((1+i*2*pi*f*T1).*(1+i*2*pi*f*T3)));
semilogx(f,20*log10(G0),'b')
hold
semilogx(f,20*log10(G1),'r')
title('Theoretical series-feedback RIAA frequency response with (blue) and without (red) T4')
xlabel('Frequency, Hz')
ylabel('Amplitude, dB')
grid; pause(3); close
 
% frequency response calculations using component values
Zriaa=R1*(1+i*2*pi*f*R2*C2)./(1+i*2*pi*f*(R1*C2+R1*C1+R2*C2)-(2*pi*f).^2*R1*C1*R2*C2);
Zs=(Rs1+Rs2).*(1+i*2*pi*f*Rs1*Rs2/(Rs1+Rs2)*Cs1)./(1+i*2*pi*f*Rs1*Cs1);
servo=(rs+(i*2*pi*f*Cs0).^-1)./Zs;
Gamp=(1+Zriaa/(R0*Rs3/(R0+Rs3)))./(1+i*2*pi*f*Rx*Cx);
Gamp=Gamp./(1+servo.*Zriaa/Rs3);
 
% frequency response derived from series feedback amp
semilogx(f,20*log10(abs(Gamp)),'k')
title('Series feedback amplifier RIAA frequency response')
xlabel('Frequency, Hz')
ylabel('Amplitude, dB')
grid; pause(3); close
 
% error frequency response derived from series feedback amp
semilogx(f,20*log10(abs(Gamp./Griaa)),'k')
title('Series feedback amplifier RIAA error frequency response')
xlabel('Frequency, Hz')
ylabel('Amplitude, dB')
grid; pause(.5);
 
% write Excel file for frequency responses
if ww==1 % write Excel file
xlswrite('RIAA_frequency_response',[f;20*log10(G0);20*log10(G1)]');
xlswrite('error_frequency_response',[f;20*log10(abs(Gamp./Griaa))]');
end
pause; close; % close graphs
 
return
% *****************************************************************

Appendix 2	Matlab routine for “tone control” evaluation

% Matlab routine for tone-control evaluation
% Tone control frequency response
% 04.08.13
close; home; clear; colordef white
fprintf('Tone-control frequency response\n')
 
% set ww = 1 to write Excel file of frequency responses
ww=0;
 
% frequency vector (1 Hz to 100 kHz, N points/dec)
N=1000; % frequency bins per decade
f=10.^(1:1/N:5); % frequency range 1 Hz to 100 kHz
dum=size(1:1/N:5);
NS=dum(2); clear dum; % number of frequency bins
 
% set data
fprintf('Load circuit and computation data\n')
maxg=5; % max gain/atten in dB
R1=5000;
RL=10000;
RH=10000;
C1=1*10^-6;
C2=.001*10^-6;
 
% display component values
mg=10^(round(maxg)/20);
R0=2*R1*(mg-1);
R2=2*R1; % identical lf and hf gain range
R0 %#ok<NOPTS>
RL %#ok<NOPTS>
RH %#ok<NOPTS>
R1 %#ok<NOPTS>
R2 %#ok<NOPTS>
C1 %#ok<NOPTS>
C2 %#ok<NOPTS>
 
% main program
fprintf('Evaluate and plot tone-control frequency response\n')
a=R0/(2*R1);
b=R0/R2;
c=RL/(8*R1);
p=0;
T0=i*2*pi*R0*C2;
T1=i*2*pi*R1*C1;
T2=i*2*pi*R2*C2;
T3=i*2*pi*RH*C2;
NP=(2*round(maxg)+1)^2;
tfdB=zeros(NP,NS);
for m1=-1:.2:1
for m2=-1:.2:1
p=p+1; % plot number
lambda1=a./((1+c*(1+m1)*(1-m1)+(c*(1+m1)*(1-m1)+.5)*f*T1));
lambda2=b*f*T2./(1+f*T2+0.25*T3*(1+m2)*(1-m2));
num=1+lambda1*(1-m1)/2+lambda2*(1-m2)/2;
den=1+lambda1*(1+m1)/2+lambda2*(1+m2)/2;
tfdB(p,1:NS)=20*log10(abs(num./den));
 
% plot tone-control frequency response
semilogx(f,tfdB(p,1:NS),'k')
hold on
end
end
title('Family of tone-control frequency responses')
xlabel('Frequency, Hz (logarithmic scale)')
ylabel('Response level, dB')
grid; pause(.5);
 
% write Excel file for frequency responses
if ww==1 % write Excel file
xlswrite('tone_control_response',[f;tfdB]');
end
pause; close; % close graphs
 
% end of program
fprintf('Program terminated\n')
return
% ******************************

Appendix 3	Matlab routine for parametric equalizer evaluation

% Parametric eq frequency response
% 04.08.13
close; home; clear; colordef white
fprintf('Parametric equalizer frequency response\n')
 
% set ww = 1 to write Excel file of frequency responses
ww=0;
 
% frequency vector (1 Hz to 100 kHz, N points/dec)
N=1000; % frequency bins per decade
f=10.^(1:1/N:5); % frequency range 10 Hz to 100 kHz
dum=size(1:1/N:5);
NS=dum(2); clear dum; % number of frequency bins
 
% set data
fprintf('Load circuit and computation data\n')
gmax(1)=15; % band 1 max boost/cut, dB
gmax(2)=15; % band 2 max boost/cut, dB
gmax(3)=15; % band 3 max boost/cut, dB
gmax(4)=15; % band 4 max boost/cut, dB
gmax(5)=15; % band 5 max boost/cut, dB
Q(1)=10; % band 1 Q-factor
Q(2)=1; % band 2 Q-factor
Q(3)=5; % band 3 Q-factor
Q(4)=3; % band 4 Q-factor
Q(5)=10; % band 5 Q-factor
fc(1)=15000; % band 1 centre frequency
fc(2)=5000; % band 2 centre frequency
fc(3)=1000; % band 3 centre frequency
fc(4)=250; % band 4 centre frequency
fc(5)=50; % band 5 centre frequency
 
% main program
fprintf('Evaluate and plot parametric equalizer frequency response\n')
q=1./Q; % 1/Q vector
gm=10.^(gmax(1:5)/20)-1; % gain vector
tfdB=zeros(31,NS); p=0;
for m=-1:1/15:1 % gain parameter
p=p+1;
num=zeros(1,NS); den=num;
for x=1:5;
num=num+.5*gm(x)*(1-m)*(i*q(x)*(f/fc(x)))./(1+i*q(x)*(f/fc(x))-(f/fc(x)).^2);
den=den+.5*gm(x)*(1+m)*(i*q(x)*(f/fc(x)))./(1+i*q(x)*(f/fc(x))-(f/fc(x)).^2);
end; 
tfdB(p,1:NS)=20*log10(abs((1+num)./(1+den)));
semilogx(f,tfdB(p,1:NS),'k')
hold on
end
title('Parametric equalizer frequency response')
xlabel('Frequency, Hz')
ylabel('Response level, dB')
grid; pause(.5);
 
% write Excel file for frequency responses
if ww==1 % write Excel file
xlswrite('parametric_response',[f;tfdB]');
end
pause; close; % close graphs

% end of program
fprintf('Program terminated\n')
return
% ***********************************