% speakercable()

function speakercable

% frequencies to evaluate, Hz
points = 100;
f = logspace(1.3, 4.3, points).';
% resistance of AWG 34, ~30, 18, 12, 9 solid wires, half of an MI cable in ohms
awg34 = roundskineffect(f, 80E-6);
awg30p5 = roundskineffect(f, 119E-6);
awg18 = roundskineffect(f, 501E-6);
awg12 = roundskineffect(f, 1.025E-3);
awg9 = roundskineffect(f, 1.45E-3);
mi1 = flatskineffect(f, 0.0254 * 0.375, 0.0254 * 0.010);
mi2 = flatskineffect(f, 0.0254 * 0.750, 0.0254 * 0.010);
mi3 = flatskineffect(f, 0.0254 * 1.250, 0.0254 * 0.012);

% find total cable loop impedance (hence the multiplies by 2 here too)
% capacitance is in parallel so it counters the phase rotation from
% inductance, same as the lumped RLCG approximation of a transmission line
% section
rDielectric = 1E12; % ohm/m
z4VS = 2 * (awg18 / 4 + 1/2.5 * j * 2 * pi * f * 596E-9 + (1 + 1/2.5 * j * 2 * pi * f * 340E-12 * rDielectric) / rDielectric);
z8VS = 2 * (awg18 / 8 + 1/2.5 * j * 2 * pi * f * 378E-9 + (1 + 1/2.5 * j * 2 * pi * f * 744E-12 * rDielectric) / rDielectric);
z8TC = 2 * (awg18 / 8 + 1/2.5 * j * 2 * pi * f * 345E-9 + (1 + 1/2.5 * j * 2 * pi * f * 812E-12 * rDielectric) / rDielectric);
zCrosslink = 2 * (awg30p5 / 108 + 3.28 * j * 2 * pi * f * 110E-9 + (1 + 3.28 * j * 2 * pi * f * 55E-12 * rDielectric) / rDielectric);
zGoldenCross = 2 * (awg34 / 752 + 3.28 * j * 2 * pi * f * 36.8E-9 + (1 + 3.28 * j * 2 * pi * f * 154E-12 * rDielectric) / rDielectric);
zNeutralReference = 2 * (awg34 / 450 + 3.28 * j * 2 * pi * f * 34E-9 + (1 + 3.28 * j * 2 * pi * f * 117E-12 * rDielectric) / rDielectric);
zMI1 = 2 * (mi1 + 3.28 * j * 2 * pi * f * 10E-9 + (1 + 3.28 * j * 2 * pi * f * 500E-12 * rDielectric) / rDielectric);
zMI2 = 2 * (mi2 + 3.28 * j * 2 * pi * f * 6E-9 + (1 + 3.28 * j * 2 * pi * f * 950E-12 * rDielectric) / rDielectric);
zMI3 = 2 * (mi3 + 3.28 * j * 2 * pi * f * 4E-9 + (1 + 3.28 * j * 2 * pi * f * 1.5E-9 * rDielectric) / rDielectric);
zZip = 2 * (awg12 + j * 2 * pi * f * 190E-9 + (1 + j * 2 * pi * f * 18E-12 * rDielectric) / rDielectric);

% plot skin effect
% loop resistance is resistance per foot of cable which includes both the
% outgoing and return paths; hence the multiplication by two
figure(1);
lines = semilogx(f, 1E3 * real(zZip), 'y', ...
                 f, 1E3 * real(zCrosslink), 'g', ... 
                 f, 1E3 * real(zNeutralReference), 'g', ...
                 f, 1E3 * real(zGoldenCross), 'g', ...
                 f, 1E3 * real(zMI1), 'r', ...
                 f, 1E3 * real(zMI2), 'r', ...
                 f, 1E3 * real(zMI3), 'r', ...
                 f, 1E3 * real(z4VS), 'b', ... 
                 f, 1E3 * real(z8VS), 'b', ...
                 f, 1E3 * real(z8TC), 'b');
ylabel('R_{loop}, m\Omega/m');
configurefigure(lines);

% plot magnitude of cable impedance
figure(2);
lines = loglog(f, 1E3 * abs(zZip), 'y', ...
               f, 1E3 * abs(zCrosslink), 'g', ...
               f, 1E3 * abs(zNeutralReference), 'g', ...
               f, 1E3 * abs(zGoldenCross), 'g', ...
               f, 1E3 * abs(zMI1), 'r', ...
               f, 1E3 * abs(zMI2), 'r', ...
               f, 1E3 * abs(zMI3), 'r', ...
               f, 1E3 * abs(z4VS), 'b', ...
               f, 1E3 * abs(z8VS), 'b', ...
               f, 1E3 * abs(z8TC), 'b');
grid('on');
ylabel('|Z_{loop}|, m\Omega/m');
configurefigure(lines);

% find fraction of load seen by power amp which is the cable
% use zero order approximation of speaker impedance
zSpeaker = 4; % ohms
cableLength = 2.43; % meters = 8 feet
cableLength = 1.8;  % meters = 6 feet
cableLength = 0.66; % meters = 26 inches
cableLength = 0.15; % meters = 6 inches

% assuming speaker is in linear drive range, then the voltage developed
% over the cable instead of the speaker is just the fraction of the total
% resistance comprised by the cable
% cable impedances are loop impedances, so no need to multiply cable length
% by two
rTermination = 0.0E-3;
lTermination = 0E-9;
zTermination = 2 * (rTermination + j * 2 * pi * f * lTermination); % ohms

lossZip = abs((cableLength * zZip + zTermination) ./ (zSpeaker + cableLength * zZip + zTermination));
loss4VS = abs((cableLength * z4VS + zTermination) ./ (zSpeaker + cableLength * z4VS + zTermination));
loss8VS = abs((cableLength * z8VS + zTermination) ./ (zSpeaker + cableLength * z8VS + zTermination));
loss8TC = abs((cableLength * z8TC + zTermination) ./ (zSpeaker + cableLength * z8TC + zTermination));
lossCrosslink = abs((cableLength * zCrosslink + zTermination) ./ (zSpeaker + cableLength * zCrosslink + zTermination));
lossGoldenCross = abs((cableLength * zGoldenCross + zTermination) ./ (zSpeaker + cableLength * zGoldenCross + zTermination));
lossNeutralReference = abs((cableLength * zNeutralReference + zTermination) ./ (zSpeaker + cableLength * zNeutralReference + zTermination));
lossMI1 = abs((cableLength * zMI1 + zTermination) ./ (zSpeaker + cableLength * zMI1 + zTermination));
lossMI2 = abs((cableLength * zMI2 + zTermination) ./ (zSpeaker + cableLength * zMI2 + zTermination));
lossMI3 = abs((cableLength * zMI3 + zTermination) ./ (zSpeaker + cableLength * zMI3 + zTermination));

% plot cable losses
figure(3);
lines = semilogx(f, 20*log10(lossZip), 'y', ...
                 f, 20*log10(lossCrosslink), 'g', ...
                 f, 20*log10(lossNeutralReference), 'g', ...
                 f, 20*log10(lossGoldenCross), 'g', ...
                 f, 20*log10(lossMI1), 'r', ...
                 f, 20*log10(lossMI2), 'r', ...
                 f, 20*log10(lossMI3), 'r', ...
                 f, 20*log10(loss4VS), 'b', ...
                 f, 20*log10(loss8VS), 'b', ...
                 f, 20*log10(loss8TC), 'b');
ylabel(sprintf('cable loss @ %.2fm, dB, Z_{T}=%.1fm\\Omega + %.1fnH', cableLength, 1E3 * rTermination, 1E9 * lTermination));
configurefigure(lines);


% set up pretty figure colors, markers, and other things
function configurefigure(lines)
hold('on');
f = subsample(lines(1), 'XData');
markerLines = plot(f, subsample(lines(1), 'YData'), 'ys', ...
                   f, subsample(lines(2), 'YData'), 'go', ...
                   f, subsample(lines(3), 'YData'), 'gd', ...
                   f, subsample(lines(4), 'YData'), 'g^', ...
                   f, subsample(lines(5), 'YData'), 'ro', ...
                   f, subsample(lines(6), 'YData'), 'rd', ...
                   f, subsample(lines(7), 'YData'), 'r^', ...
                   f, subsample(lines(8), 'YData'), 'bo', ...
                   f, subsample(lines(9), 'YData'), 'bd', ...
                   f, subsample(lines(10), 'YData'), 'b^');
hold('off');
set(markerLines, 'MarkerSize', 5);
grid('on');
xlabel('frequency, Hz');
legend(markerLines, 'AWG 12 zip', 'Cardas Crosslink', 'Cardas Neutral Reference', 'Cardas Golden Cross', 'Goertz MI1', 'Goertz MI2', 'Goertz MI3', 'Kimber 4VS', 'Kimber 8VS', 'Kimber 8TC', 2);

set(gcf, 'Color', [ 0.1, 0.1, 0.1 ]);
set(gca, 'Color', [ 0.1, 0.1, 0.1 ]);
set(legend, 'Color', [ 0.4, 0.4, 0.4 ]);

set(gca, 'XColor', [ 0.8, 0.8, 0.8 ]);
set(gca, 'YColor', [ 0.8, 0.8, 0.8 ]);

set(gca, 'XLim', [10, 30E3]);


function subsampledData = subsample(line, data)
skip = 15;
values = get(line, data);
subsampledData = values(round(0.3*skip):skip:length(values));


% do numerical integration of current falloff with increasing depth due to
% skin effect for a round, solid wire
% integration here is a zero order approximation with sampling at the edges
% of the bands; quite crude, but skin effect is smooth and a lot of points
% are thrown at the problem so the error here is not huge and the results
% are a decent first order approximation of the resistance of a stranded
% cable
function r = roundskineffect(f, wireRadius)

% conductivity of copper
sigma = 58.6E6; % Siemens

% number of points along radius of wire
integrationPoints = 100;
depth = linspace(0, wireRadius, integrationPoints).';

ringRadius = linspace(wireRadius, 0, integrationPoints).';
ringSurfaceArea = pi * (ringRadius.^2 - [ ringRadius(2:length(ringRadius)); 0 ].^2);

% copper, gold, etc, all have permeability of free space
mu0 = 4*pi*1E-7;
% skin depth as a function of frequency
skinDepth = 1 ./ sqrt(pi*f*mu0*sigma);

% run numerical integration
r = ones(size(f));
for frequency = 1:length(f)
    currentDensity = exp(-depth / skinDepth(frequency));
    weightedSurfaceArea = sum(currentDensity .* ringSurfaceArea);
    r(frequency) = 1 / (weightedSurfaceArea * sigma);
end


% do numerical integration of current falloff for Goertz planar cables
% essentially the same as the code above, but edge effects at the edges of
% the metal ribbons are neglected and the assumption is all current
% clusters against the inside faces of the two traces---this quite
% reasonable as it's the lowest impedance position and the back face
% current is almost negligible for any microstrip type structure at 
% operating at a reasonably high frequency
function r = flatskineffect(f, width, thickness)

sigma = 58.6E6;

integrationPoints = 100;
depth = linspace(0, thickness, integrationPoints).';
stripThickness = thickness / (integrationPoints - 1);

mu0 = 4*pi*1E-7;
skinDepth = 1 ./ sqrt(pi*f*mu0*sigma);

r = ones(size(f));
for frequency = 1:length(f)
    currentDensity = exp(-depth / skinDepth(frequency));
    weightedSurfaceArea = sum(currentDensity .* stripThickness) * width;
    r(frequency) = 1 / (weightedSurfaceArea * sigma);
end