Bessel type function? The matrix below
$\left(
\begin{array}{ccccc}
1 & 0 & 1 & 0 & 1 & 0 & 1 & \dots\\
\frac{1}{2} & \frac{1}{2} & 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 & \dots\\
\frac{1}{3} & \frac{1}{3} & \frac{1}{3} & 0 & 0 & 0 & \frac{1}{3} & \dots\\
\frac{1}{4} & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & 0 & 0 & 0 & \dots\\
\frac{1}{5} & \frac{1}{5} & \frac{1}{5} & \frac{1}{5} & \frac{1}{5} & 0 & 0 & \dots\\
\frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & 0 & \dots\\
\frac{1}{7} & \frac{1}{7} & \frac{1}{7} & \frac{1}{7} & \frac{1}{7} & \frac{1}{7} & \frac{1}{7} & \dots\\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots\\
\end{array}
\right)$
follows a pattern whereby the first line begins $\{1,0,1,0,1,0,1,0,1,0,1,0\dots\}$
the second: $\{\frac{1}{2},\frac{1}{2},0,0,\frac{1}{2},\frac{1}{2},0,0,\frac{1}{2},\frac{1}{2},0,0\dots\}$
the third: $\{\frac{1}{3},\frac{1}{3},\frac{1}{3},0,0,0,\frac{1}{3},\frac{1}{3},\frac{1}{3},0,0,0\dots\}$
and so on.
Summing the first $m$ rows together and subtracting $\frac{1}{2}\sum_{k=1}^{m}\frac{1}{k}$ from each element from the resultant array gives a sequence of numbers that fluctuates between positive and negative values. The resultant cummulative plot looks like this:

y = 3000; m = 500;
ListPlot[Accumulate[Total[Take[Flatten[ConstantArray
[#,Ceiling[(y)/Length@#]]], y] & /@ Table[Join[ConstantArray[1/row, row], 
ConstantArray[0, row]], {row, 1, m}]] - Sum[1/k, {k, 1, m}]/2]]

where $y$ is the number of columns, and $m$ the number of rows, and is apparently scalable by dividing by $\dfrac{m}{4}$:

Is it true that this is related to the family of Bessel functions (inclusive of trigonometric integrals, which are related to spherical Bessel functions)?
Below is a plot of a Bessel function of the second kind, for illustrative purposes:

given by
$$
Y_0(z)=\dfrac{2}{\pi}\bigg(\bigg(\log(z/2)+\gamma\bigg)\sum_{m=0}^{\infty}\dfrac{(-1)^m}{m!\ \Gamma(m+1)}(z/2)^{2m}+\sum_{k=1}^{\infty}(-1)^{k+1}\sum_{i=1}^{\infty}\dfrac{1}{i}\bigg(\dfrac{(z^2/4)^k}{(k!)^2}\bigg)\bigg)
$$
Whether it is the case or not, I would be interested to know why it follows this pattern. The periodicity up to $y$ follows a particular type of randomness, seen also in the differences between zeta zeros, for example. I am curious to find out whether there is any connection.
Update

plotted with this (ridiculously long) piecewise code, thanks to Ruslan's answer below.
 A: This is unlikely to be related to Bessel functions or Fresnel integral or sine/cosine integral. Let's take your code, assigning the ListPlot argument to variable data:
y = 3000; m = 500;
data = Accumulate[
   Total[Take[Flatten[ConstantArray[#, Ceiling[(y)/Length@#]]], y] & /@
       Table[Join[ConstantArray[1/row, row], 
        ConstantArray[0, row]], {row, 1, m}]] - 
    Sum[1/k, {k, 1, m}]/2];
ListPlot[data]


Now we'd need to smooth this somewhat. Let's use a moving average with e.g. 30 points:
smdata = ListConvolve[#/Total[#] &@Table[1., {i, 1, 30}], data];
ListPlot[smdata, PlotRange -> All]


Now this looks somewhat more tractable by differentiating, although it does have some noise, we can still see some structure in the derivative. For this, make an interpolation:
di = Interpolation[smdata];

And plot the derivative of the interpolant:
Plot[di'[x], {x, 1, First@Dimensions@data - 500}, PlotRange -> All]


We can see that the derivative may even have some cusps at its extrema, so it doesn't look like a smooth function, which would be if it were a Bessel function or Fresnel integral or anything supposed in the comments.
The plot, however, gives some ideas. What if it's just "broken" and (quasi)periodically reflected? Let's try adding -0.45-di'[x] to the plot:
Plot[{di'[x], -0.45 - di'[x]}, {x, 1, First@Dimensions@data - 500}, PlotRange -> All]


This seems to smoothly fit the guess. I think we could do the similar thing with other cusps and recover something like a logarithmic function. Here're some more added plots:
Plot[{di'[x], -0.45 - di'[x], di'[x] - 0.70, -0.85 - di'[x], di'[x] - 1}, {x, 1, 
  First@Dimensions@data - 500}, PlotRange -> All]


Defining a piecewise function, we can plot the single curve:
pw[x_] = Piecewise[{
    {di'[x], x < 500},
    {-0.45 - di'[x], x < 1000},
    {di'[x] - 0.70, x < 1500},
    {-0.85 - di'[x], x < 2000},
    {di'[x] - 1, True}}];
Plot[pw[x], {x, 1, First@Dimensions@data - 500}, PlotRange -> All]


Now, more interesting. Let's define analytically what we get. First, we can notice that our matrix can be generated by the following function:
$$f(x,y)=\frac1{2y}\left\{1+\text{sgn}\left[\sin\left(\frac{\pi x}y-\varepsilon\right)\right]\right\},$$
where to get the matrix we should take $\lim_{\varepsilon\to+0}$ and use $x,y=1,2,...$.
Now, we aren't interested in any pixellation noise, so we can switch from taking sums to taking integrals, and from matrix to the generator function itself. Then we don't have to bother with $\varepsilon$, and can set $\varepsilon=0$. Now the sum over column index $y$ and subtracting the sum over $\frac1k$ will give us a function of row index $x$:
$$s(x)=\int_1^m\frac1{2y}\text{sgn}\left(\sin\frac{\pi x}y\right)\text{d} y,$$
where $m$ is the same as m in your code.
Now, taking $y=300$, $m=50$, I get the resulting $s(x)$:

Here, at, e.g. $x=30$, we have:
$$s(30)=-\frac12\left(\log\frac35+\sum_{n=2}^{30}(-1)^n \log\frac n{n-1}\right).$$
Similar expressions in terms of logarithms are for other $x$ values.
So, the correct description of the function is in terms of simple logarithms. Then the function which you thought as Bessel-related, being an integral of this function, is also expressible in terms of logarithms.
Note however, that the function isn't really smooth. See it magnified:

