Autocorrelation function of a triangular wave I am interested to find the analytical expression for the autocorrelation function of a signal that comprises triangular pulses. I have followed the derivation in &quot;Statistical Theory of Communication&quot; by Lee (1960, p.16) 
However, if I replicate this in MATLAB, then I do not obtain the autocorrelation function as derived here. The general shape is the same, but it looks as follows:

with the autocorrelation function dipping below zero. In this case, the triangular pulses are spaced such that each triangle is half of the period, with the remainder of each period being zero.
Why is this different from Lee's Derivation, and how can I go about adapting his technique to find the analytical form that matches the empirical result?
The code I used in MATLAB was as follows:
function [x, acf1, acf2] = ACFTriangle(leftbase,rightbase,height,gap,repeat)


x = [zeros(1,gap) (1:leftbase)*(height/leftbase) (rightbase-1:-1:1)*(height/rightbase) 0];

x = repmat(x,[1,repeat]);

t = 0:length(x)-1;

acf = xcov(x,'unbiased');

n = length(acf);

s = (n-1)/2 + 1;

acf1 = acf(s:n)/var(x);

acf2 = autocorr(x,length(x)-1);

x = [t; x]';

end

The reason for my question is that I have a very large number of time series in which there may be a hidden triangular pulse wave. If I have an analytical expression for the autocorrelation sequence, then I can identify their characteristics from the shape of the autocorrelation function.
New edit
Notwithstanding the excellent answer below, how would one go about modifying Lee's mathematical derivation to allow for removal of the means? I can attempt this doing the following:
\begin{equation}
\varphi_{11} = \frac{1}{T_1}\int_0^{b-\tau} \left(\frac{E_m}{b}t - m\right)\left(\frac{E_m}{b}\left(t+\tau\right) - m\right)\,dt
\end{equation}
where $m$ is the mean of the series.  However, following this through gives:
\begin{equation}
\varphi_{11} = \frac{E_m^2}{6b^2T_1}\left(\tau^3-3b^2\tau + 2b^3\right) + m^2\left(b-\tau\right) - \frac{mE_m}{bT_1}\left(b-\tau\right)^2 - \frac{mE_m}{bT_1}\tau\left(b-\tau\right)
\end{equation}
which clearly gives $\varphi_{11}=0$ for $\tau=b$.  Where is the error in this working?
I couldn't quite follow your code (you calculate two different versions of the ACF?), but I believe the problem with the plot being shifted toward zero is that xcov calculates the cross covariance not correlation. Remember that cross covariance is subtracting the mean, so there should be a shift! You should be calling xcorr instead.
The following code snippet reproduces the Figure 9 from your screenshot.
tau = 1;
b = 100;
Em = sqrt(6);
T1 = 2*b;
signal = zeros(T1, 1);
step = Em/(b+1);
signal(tau:tau+b) = linspace(0, Em, b+1);
signal = repmat(signal, 1000, 1);
[xcorrOut, lag] = xcorr(signal, 'unbiased');

figure
hold on
grid on
plot(lag, xcorrOut)
xlim([-5*b, 5*b])

Edit
In response to the new edit in the question, the OP wants a way to convert the cross correlation result to a cross covariance result. The cross covariance can be obtained from the cross correlation by subtracting the product of the mean and conjugate mean, https://en.wikipedia.org/wiki/Cross-covariance:
$$ \varphi_{11}(\tau) - |\mu_{f_1}|^2$$
where $\mu_{f_1}$ is the mean of $f_1(t)$, as shown in Figure 7 of the screen capture. This additional code snippet builds upon the above to show how adjusting the correlation can get you the same result as computing the covariance.
mu = mean(signal);
xcovFromXcorr = xcorrOut - mu*conj(mu);
[xcovOut, lag] = xcov(signal, 'unbiased');

figure
hold on
grid on
plot(lag, xcovFromXcorr, 'LineWidth', 2)
plot(lag, xcovOut, '--', 'LineWidth', 2)
xlim([-5*b, 5*b])