Wave Equation Boundary Conditions

I am trying to apply the boundary conditions to the wave equation given by: $$c^2\frac{\partial^2 u}{\partial x^2} = \frac{\partial^2 u}{\partial t^2}$$ What I am interested is the Neumann boundary conditions, but my question applies to any type of BCs.

It is straight forward to show that the solutions of the wave equation have the form: $$u_i = A_i e^{j(k_i x - \omega_i t)} + B_i e^{j(k_i x + \omega_i t)}$$ Where $$A_i$$, and $$B_i$$ are constant, $$j = \sqrt{-1}$$, $$c k_i = \omega_i$$, and $$i \in [-\infty, \infty]$$. I have seen a variety of sources where people apply the boundary conditions at this point, e.g. 28, and 29 Wolfram Mathworld. But surely the actual solution is: $$u = \sum_i u_i = \sum_i \left[ A_i e^{j(k_i x - \omega_i t)} + B_i e^{j(k_i x + \omega_i t)} \right]$$ Can you not have a case where $$A_i e^{j(k_i x - \omega_i t)} + B_i e^{j(k_i x + \omega_i t)}$$ on its own does does not satisfy the boundary conditions, but when you sum over all the frequencies you get a solution that satisfies them? Why do you assume that $$u_i$$ on its own must satisfy the boundary conditions?

In fact, one can take the Fourier transform of the sum in the time dimension, using the fact $$k_i = -k_{-i}$$ in order to flip the order of the sum over $$A_i e^{j(k_i x - \omega_i t)}$$: $$U = \sum_i \delta(\omega - \omega_i) \left( A_{-i} e^{-jk_i x} + B_{i} e^{j k_i x} \right)$$ So this is the frequency domain representation sampled at every frequency because the Dirac delta forms a continuous comb function. One can turn the sum into a continuous function: $$U = A(-\omega) e^{-jx \omega / c} + B(\omega) e^{jx \omega / c}$$ Applying the inverse Fourier transform gives the familiar solution in time domain: $$u = a(-t + x/c) + b(t + x/c)$$ I would then say that you need to apply the boundary conditions on $$U = A(-\omega) e^{-jx \omega / c} + B(\omega) e^{jx \omega / c}$$. I can certainly see that for $$u(0, t) = 0$$ this results in having each frequency component equal to zero. But that wouldn't be the case for $$u(0, t) = f(t)$$.

If $f(t)\in\bf{L}^{2}(\mathbb{R})$ you expand the solution in the uncontable basis of the infinite dimentional Hilbert space of solutions $$u(x, t)=\int_{\mathbb{R}^{2}}\tilde{u}_{-}(k, \omega)e^{i\omega{t}-ikx}d\omega{dk}+\int_{\mathbb{R}^{2}}\tilde{u}_{+}(k, \omega)e^{i\omega{t}+ikx}d\omega{dk}$$ then the coefficients $\tilde{u}_{\pm}(k, \omega)$ are given by the inner-product $$\tilde{u}_{\pm}(k, \omega)=\hat{u}_{\pm}(k)\frac{1}{(2\pi)}\int_{\mathbb{R}}u(0, t)e^{-i\omega{t}}dt$$ Where, say $$\hat{u}_{\pm}(k)=\frac{1}{2\pi}\int_{\mathbb{R}}u(x, 0)e^{\mp{ikx}}dx$$ This is for the Dirichlet boundary conditions. For neumann you will have the derrivatives of solutions specified at the boundary, but procedure is clearly similar. If it happends that the u(0, t) is defined on some interval, not the entire line (or it is periodic ofcourse), the basis for the Hilbert space becomes countable and you replace the integral by the sum, so is for the $u(x, 0)$.