Let the camera be placed at the origin. Let its front direction be the z-axis, its upward direction be the y-axis and its right direction be the x-axis. This is a left-handed system. We work in homogeneous coordinates. Thus a point (x,y,z)\in\mathbb{R}^{3} becomes (x:y:z:1)\in\mathbb{P}^{3} in the homogeneous coordinates.
Let the near plane be z=z_{0} and the far plane be z=z_{1}. We assume that 0<z_{0}<z_{1}. We consider the perspective projection to the near plane: \begin{aligned} (x,y,z) \mapsto (\frac{xz_{0}}{z}, \frac{yz_{0}}{z}, z_{0}) \quad\text{for $z\neq 0$} \end{aligned} or in homogeneous coordinates \begin{aligned} (x:y:z:w) \mapsto (\frac{xz_{0}}{z}: \frac{yz_{0}}{z}: z_{0} :1) = (z_{0}x: z_{0}y: z_{0}z: z) \end{aligned} which then becomes a linear map. The disadvantage is that we lose the depth information. Thus we make the following adjustment in the z-coordinate \begin{aligned} (x:y:z:w) \mapsto (\frac{xz_{0}}{z}: \frac{yz_{0}}{z}: z_{0}- \frac{w}{z} :1) = (z_{0}x: z_{0}y: z_{0}z - w: z). \end{aligned} This is still a linear map and the matrix is given by \begin{aligned} P:= \begin{pmatrix} z_{0} & 0 & 0 & 0 \\ 0 & z_{0} & 0 & 0 \\ 0 & 0 & z_{0} & -1 \\ 0 & 0 & 1 & 0 \end{pmatrix}. \end{aligned} Note that points with coordinates (x,y,z) for z\in[z_{0},z_{1}] are mapped to \mathbb{R}\times \mathbb{R}\times [z_{0}- \frac{1}{z_{0}}, z_{0} - \frac{1}{z_{1}}]. We note that z_{0}- \frac{1}{z_{0}} < z_{0} - \frac{1}{z_{1}}.
Next we find a linear map that transforms [x_{0},x_{1}]\times [y_{0},y_{1}] \times [z_{0}- \frac{1}{z_{0}}, z_{0} - \frac{1}{z_{1}}] to [-1,1]\times[-1,1]\times[-1,1]. We want to preserve the orientation of each pair of corresponding line segments. It is easy to see that we should take \begin{aligned} &\rightarrow [-1,1] \\ x &\mapsto \frac{2x - x_{0} - x_{1}}{x_{1} - x_{0}} \end{aligned} and similarly for other coordinates. Thus in the homogeneous coordinates, the linear map is \begin{aligned} (x:y:z:w) \mapsto &(\frac{2x - (x_{0} + x_{1})w}{(x_{1} - x_{0})w} : \frac{2y - (y_{0} + y_{1})w}{(y_{1} - y_{0})w} : \frac{ 2z_{0}z_{1}z - (2z_{0}^{2}z_{1}-z_{0} - z_{1})w}{(z_{1} - z_{0})w} : 1) \\ &= (\frac{2x - (x_{0} + x_{1})w}{(x_{1} - x_{0})} : \frac{2y - (y_{0} + y_{1})w}{(y_{1} - y_{0})} : \frac{ 2z_{0}z_{1}z - (2z_{0}^{2}z_{1}-z_{0} - z_{1})w}{(z_{1} - z_{0})} : w) \end{aligned} and the matrix is \begin{aligned} Q:= \begin{pmatrix} \frac{2}{x_{1}-x_{0}} & 0 & 0 & - \frac{x_{0}+x_{1}}{x_{1}-x_{0}} \\ 0 & \frac{2}{y_{1}-y_{0}} & 0 & - \frac{y_{0}+y_{1}}{y_{1}-y_{0}} \\ 0 & 0 & \frac{2z_{0}z_{1}}{z_{1}-z_{0}} & -\frac{2z_{0}^{2}z_{1}-z_{0} - z_{1}}{z_{1} - z_{0}} \\ 0 & 0 & 0 & 1 \end{pmatrix}. \end{aligned}
Consider the frustum which is the part of the polygonal cone consisting of rays from the origin through points in [x_{0},x_{1}]\times [y_{0},y_{1}] \times \{z_{0}\} between the near plane z=z_{0} and the far plane z=z_{1}. Its image under the linear transformation given by QP is then [-1,1]\times[-1,1]\times[-1,1]. We have \begin{aligned} QP = \begin{pmatrix} \frac{2z_{0}}{x_{1}-x_{0}} & 0 & - \frac{x_{0}+x_{1}}{x_{1}-x_{0}} & 0 \\[.5em] 0 & \frac{2z_{0}}{y_{1}-y_{0}} & - \frac{y_{0}+y_{1}}{y_{1}-y_{0}} & 0 \\[.5em] 0 & 0 & \frac{z_{0}+z_{1}}{z_{1}-z_{0}} & -\frac{2z_{0}z_{1}}{z_{1} - z_{0}} \\[.5em] 0 & 0 & 1 & 0 \end{pmatrix}. \end{aligned} This is the projection matrix. This matrix is defined up to scalar multiplication.
For OpenGL, the z-axis is in the opposite direction to our setup. Let the z'-axis be the opposite direction to our z-axis. Let z'=-z_{0}' be the near plane and z=-z_{1}' be the far plane with 0 < z_{0}' < z_{1}'. Thus the projection matrix from the view space to the NDC (normalised device coordinates) in OpenGL is given by \begin{aligned} QP \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \end{aligned} where z_{0} = z_{0}' and z_{1} = z_{1}', namely \begin{aligned} \begin{pmatrix} \frac{2z_{0}}{x_{1}-x_{0}} & 0 & \frac{x_{0}+x_{1}}{x_{1}-x_{0}} & 0 \\[.5em] 0 & \frac{2z_{0}}{y_{1}-y_{0}} & \frac{y_{0}+y_{1}}{y_{1}-y_{0}} & 0 \\[.5em] 0 & 0 & -\frac{z_{0}+z_{1}}{z_{1}-z_{0}} & -\frac{2z_{0}z_{1}}{z_{1} - z_{0}} \\[.5em] 0 & 0 & -1 & 0 \end{pmatrix}. \end{aligned} This is the projection matrix in the OpenGL setting.