Let $H$ be a $\mathbb C$-Hilbert space, $A\in\mathfrak L(H)$ be self-adjoint and $B\in\mathfrak L(H)$.
If $\xi_A$ is the spectral measure associated with $A$, I would like to show that $A$ and $B$ commute iff $\xi_A(I)$ and $B$ commute for all $I\in\mathcal B(\mathbb R)$.
For simplicity, I would like to assume $\dim H\in\mathbb N$. Letting $\pi_A(\lambda)$ denote the orthogonal projection onto $E_A(\lambda):=\mathcal N(\lambda-A)$, the desired claim should become:
- $AB=BA$;
- $\pi_A(\lambda)B=B\pi_A(\lambda)$ for all $\lambda\in\sigma(A)\setminus\{0\}$.
Since we can write $$A=\sum_{\lambda\in\sigma(A)\setminus\{0\}}\lambda\pi_A(\lambda)\tag1,$$ it is clear that (2.) implies (1.).
But I really struggle to show the converse.
Assume $(1.)$. Let $x\in H$ and $\lambda\in\sigma(A)\setminus\{0\}$. Maybe it's helpful to consider the cases $x\in E_A(\lambda)$ and $x\not\in E_A(\lambda)$ separately. In the first case, $$B\pi_A(\lambda)x=BAx=ABx=\sum_{\mu\in\sigma(A)\setminus\{0\}}\pi_A(\mu)Bx\tag2.$$ However, I've don't know how we can show that the right-hand side is equal to $\pi_A(\lambda)Bx$ ...
Since you are staying in the finite-dimensional case, you can do the following.
From $AB=BA$ you can easily get $A^nB=BA^n$ for all $n\in\mathbb N$. It then follows that $p(A)B=Bp(A)$ for all $p\in\mathbb C[x]$.
Now with a fixed $\lambda\in\sigma(A)$ (doesn't matter if it is nonzero), since $\sigma(A)$ is finite you can choose a polynomial $p$ such that $p(\lambda)=1$ and $p(\eta)=0$ for all $\eta\in\sigma(A)\setminus\{\lambda\}$. Then $$ p(A)=\sum_{\alpha\in\sigma(A)}p(\alpha)\,\pi_A(\alpha)=\pi_A(\lambda) $$ and $$ \pi_A(\lambda)B=p(A)B=Bp(A)=B\pi_A(\lambda). $$ The same idea can be used in infinite dimension but it requires the Borel functional calculus and to restate your problem a bit since $A$ may not have any eigenvalues.